x y and z velocity coordinates

ncc1701d

Member
Joined
Aug 9, 2009
Messages
204
Reaction score
6
Points
18
Can someone tell me how I can find the X velocity coordinate, Y velocity coordinate, and Z velocity coordinates if I am given the velocity?
I believe I am looking for the Vx Vy Vz shown in this picture. https://en.wikipedia.org/wiki/Orbital_state_vectors#/media/File:Orbital_state_vectors.png
The video here
does a a good job explaining way to find velocity .
But I cant find much detail on web or examples on finding the X Y Z velocity coordinate if I am given the velocity.
I think you can add them up to get the velocity but that doesnt help me much. Any help apprecitated. Thank you.
 

asbjos

tuanibrO
Addon Developer
Joined
Jun 22, 2011
Messages
696
Reaction score
259
Points
78
Location
This place called "home".
It is not possible to find the xyz velocity vector components from the magnitude (speed) only. You need some other information in addition to determine the components.

Without further information, it's like saying "I'm driving at 80 km/h, which direction am I going?" You may be driving southward, westward, uphill to the right, or any other direction. Without any extra information, it's impossible to determine.

The orbital parameters are heavily interwoven, so even the most obscure information may help finding the components. Are you sure that you have no additional knowledge of the orbit at hand?
 

ncc1701d

Member
Joined
Aug 9, 2009
Messages
204
Reaction score
6
Points
18
It is not possible to find the xyz velocity vector components from the magnitude (speed) only. You need some other information in addition to determine the components.

Without further information, it's like saying "I'm driving at 80 km/h, which direction am I going?" You may be driving southward, westward, uphill to the right, or any other direction. Without any extra information, it's impossible to determine.

The orbital parameters are heavily interwoven, so even the most obscure information may help finding the components. Are you sure that you have no additional knowledge of the orbit at hand?
Ok Here is the information I have. I have spacecraft above orbiting Venus. The following info is all the info I have. The reference frame is going to be Venus. I have a 2 points. They are of where it was and where it is going. The 1st point has a known altitude above Venus and I know the X Y Z position cartesian coordinates of the spacecraft using Venus as the reference frame. I also know its latitude and longitude angles referenced from Venus. I also know the same information about the second point as I do about the first. I also know the time it takes to travel from one point to the next. My challenge is I need the velocity vector components of X Y and Z. Do I have enough information to find what I need?
 

francisdrake

Addon Developer
Addon Developer
Joined
Mar 23, 2008
Messages
1,060
Reaction score
864
Points
128
Website
francisdrakex.deviantart.com
You could try playing around with a scenario (.scn) file.
If a vessel is in orbit, the line RPOS gives the relativ position to the reference body in cartesian coordinates im [m],
while RVEL gives its relative velocity in [m/s].

If the position is given in polar coordinates (long, lat, rad) it needs to be transferred to cartesian coordinates.
In this case keep in mind Orbiter uses a left-handed coordinate system.

BEGIN_SHIPS GL-01:DeltaGlider STATUS Orbiting Earth RPOS 3626158.96 4307928.18 -3325004.36 RVEL 6623.108 -3432.497 2656.884 AROT -52.67 -56.93 90.32 ...
 

asbjos

tuanibrO
Addon Developer
Joined
Jun 22, 2011
Messages
696
Reaction score
259
Points
78
Location
This place called "home".
Ok Here is the information I have. I have spacecraft above orbiting Venus. The following info is all the info I have. The reference frame is going to be Venus. I have a 2 points. They are of where it was and where it is going. The 1st point has a known altitude above Venus and I know the X Y Z position cartesian coordinates of the spacecraft using Venus as the reference frame. I also know its latitude and longitude angles referenced from Venus. I also know the same information about the second point as I do about the first. I also know the time it takes to travel from one point to the next. My challenge is I need the velocity vector components of X Y and Z. Do I have enough information to find what I need?
Now we're talking! Yes, from an initial position, final position and time from initial to final position, you can determine the velocity vector (or any other parameter) of the orbit. This challenge is called Lambert's problem.
Wiki: https://en.wikipedia.org/wiki/Lambert's_problem
This page may also be of interest: https://space.stackexchange.com/que...ments-from-two-position-vectors-and-a-time-di

I have not done this before myself, though.
As I have done lots of conversion from Kepler elements to state vectors and back, I would have tried to first find the orbital elements by solving Lambert's problem. Then you can find the state vectors (and therefore the velocity vector) at any time by varying the true anomaly (i.e. time). (here are the methods for state vector to Kepler and Kepler to state vector)
But be wary that this may not be the most efficient solution. There may be some more direct way that I'm not aware of.

If you need actual help calculating for actual position vectors, then let me know. (no guarantees though, as I haven't done this before!)
 

ncc1701d

Member
Joined
Aug 9, 2009
Messages
204
Reaction score
6
Points
18
It is not possible to find the xyz velocity vector components from the magnitude (speed) only. You need some other information in addition to determine the components.

Without further information, it's like saying "I'm driving at 80 km/h, which direction am I going?" You may be driving southward, westward, uphill to the right, or any other direction. Without any extra information, it's impossible to determine.

The orbital parameters are heavily interwoven, so even the most obscure information may help finding the components. Are you sure that you have no additional knowledge of the orbit at hand?

Now we're talking! Yes, from an initial position, final position and time from initial to final position, you can determine the velocity vector (or any other parameter) of the orbit. This challenge is called Lambert's problem.
Wiki: https://en.wikipedia.org/wiki/Lambert's_problem
This page may also be of interest: https://space.stackexchange.com/que...ments-from-two-position-vectors-and-a-time-di

I have not done this before myself, though.
As I have done lots of conversion from Kepler elements to state vectors and back, I would have tried to first find the orbital elements by solving Lambert's problem. Then you can find the state vectors (and therefore the velocity vector) at any time by varying the true anomaly (i.e. time). (here are the methods for state vector to Kepler and Kepler to state vector)
But be wary that this may not be the most efficient solution. There may be some more direct way that I'm not aware of.

If you need actual help calculating for actual position vectors, then let me know. (no guarantees though, as I haven't done this before!)
thank you. Lots of good information to check out in your post. I was suspecting I needed more information about my orbit ellipse to get my answer. Looks like Lambert is the answer.
 

ncc1701d

Member
Joined
Aug 9, 2009
Messages
204
Reaction score
6
Points
18
Now we're talking! Yes, from an initial position, final position and time from initial to final position, you can determine the velocity vector (or any other parameter) of the orbit. This challenge is called Lambert's problem.
Wiki: https://en.wikipedia.org/wiki/Lambert's_problem
This page may also be of interest: https://space.stackexchange.com/que...ments-from-two-position-vectors-and-a-time-di

I have not done this before myself, though.
As I have done lots of conversion from Kepler elements to state vectors and back, I would have tried to first find the orbital elements by solving Lambert's problem. Then you can find the state vectors (and therefore the velocity vector) at any time by varying the true anomaly (i.e. time). (here are the methods for state vector to Kepler and Kepler to state vector)
But be wary that this may not be the most efficient solution. There may be some more direct way that I'm not aware of.

If you need actual help calculating for actual position vectors, then let me know. (no guarantees though, as I haven't done this before!)
If I look at your link https://en.wikipedia.org/wiki/Lambert's_problem Do you know where s1 and s2 come from? They dont illustrate it in their diagrams. Iam a visual learner so kinda need to see it. I dont have enough info to find F2 from this do I ?
 
Last edited:

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,651
Reaction score
785
Points
128
I got a source codes for a multiple different Lambert solvers, so, I can copy paste one if you like. But if you want to understand what you are doing then I can write a description to the very first I used which was my own design and comes from very simple basic equations. This is probably the only one I can understand. The source code for this been long gone.

Orbital Radius Is:
\[ r = \frac {p}{1 + e\cos(\upsilon)} \]

From here we can solve orbital parameter "p"
\[ p = r(1+e\cos(\upsilon)) \]

Next step is to "mark" two of them equal
\[ r_1(1+e\cos(\upsilon_1)) = r_2(1+e\cos(\upsilon_2)) \]

From here we can solve eccentricity "e"
\[ e = \frac {r_1 - r_2}{r_2\cos(\upsilon_2) - r_1\cos(\upsilon_1)} \]

If \( \beta \) is an angle between the two radius vectors and \( \upsilon_1 \) is you "true anomaly" which you have to guess to give you the required time of flight. I'd say that the eccentricity equation is the main equation for the solution.

\[ \upsilon_2 = \beta + \upsilon_1 \]

And to compute the Time of Flight:

Semi-major axis:
\[ a = \frac {p}{1-e^2} \]

Eccentric anomaly is:
\[ E_x = \cos^{-1}\frac{\cos(\upsilon_x) + e}{1+e\cos(\upsilon_x)} \]

Mean Anomaly is:
\[ M_x = E_x - e\sin(E_x) \]

And finally Time of Flight:
\[ T = |M_1-M_2| \sqrt{\frac {|a^3|}{\mu}} \]

C:
double tra2eca(double tra, double e)
{
    double eca;
    double cos_tra = cos(tra);
       
    if (e>1.0) {
        eca = acosh( (e+cos_tra) / (1.0+e*cos_tra) );  
        if (tra>PI) eca=-eca;
        return(eca);
    }

    eca = acos( (e+cos_tra) / (1.0+e*cos_tra) );
    if (tra>PI) eca = PI2 - eca;
    return(eca);
}

double eca2mna(double eca, double ecc)
{
    if (ecc<1.0) return eca - ( ecc * sin( eca ) );
    return ( ( ecc * sinh(eca) ) - eca );  
}
 

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,651
Reaction score
785
Points
128
This is also a classic book "BMW" among anyone interested about orbital mechanics over here. Mine has been read so much that covers have almost worn-off. There are three different lambert solvers discussed.

Fundamentals of Astrodynamics
 

ncc1701d

Member
Joined
Aug 9, 2009
Messages
204
Reaction score
6
Points
18
I got a source codes for a multiple different Lambert solvers, so, I can copy paste one if you like. But if you want to understand what you are doing then I can write a description to the very first I used which was my own design and comes from very simple basic equations. This is probably the only one I can understand. The source code for this been long gone.

Orbital Radius Is:
\[ r = \frac {p}{1 + e\cos(\upsilon)} \]

From here we can solve orbital parameter "p"
\[ p = r(1+e\cos(\upsilon)) \]

Next step is to "mark" two of them equal
\[ r_1(1+e\cos(\upsilon_1)) = r_2(1+e\cos(\upsilon_2)) \]

From here we can solve eccentricity "e"
\[ e = \frac {r_1 - r_2}{r_2\cos(\upsilon_2) - r_1\cos(\upsilon_1)} \]

If \( \beta \) is an angle between the two radius vectors and \( \upsilon_1 \) is you "true anomaly" which you have to guess to give you the required time of flight. I'd say that the eccentricity equation is the main equation for the solution.

\[ \upsilon_2 = \beta + \upsilon_1 \]

And to compute the Time of Flight:

Semi-major axis:
\[ a = \frac {p}{1-e^2} \]

Eccentric anomaly is:
\[ E_x = \cos^{-1}\frac{\cos(\upsilon_x) + e}{1+e\cos(\upsilon_x)} \]

Mean Anomaly is:
\[ M_x = E_x - e\sin(E_x) \]

And finally Time of Flight:
\[ T = |M_1-M_2| \sqrt{\frac {|a^3|}{\mu}} \]

C:
double tra2eca(double tra, double e)
{
    double eca;
    double cos_tra = cos(tra);
     
    if (e>1.0) {
        eca = acosh( (e+cos_tra) / (1.0+e*cos_tra) );
        if (tra>PI) eca=-eca;
        return(eca);
    }

    eca = acos( (e+cos_tra) / (1.0+e*cos_tra) );
    if (tra>PI) eca = PI2 - eca;
    return(eca);
}

double eca2mna(double eca, double ecc)
{
    if (ecc<1.0) return eca - ( ecc * sin( eca ) );
    return ( ( ecc * sinh(eca) ) - eca );
}
Thank you but I am not sure we are on the same page. I already have the time of flight between between the two points and the two points position in space. My understanding is I need at least the semi-major axis to reconstruct the orbital ellipse which fits between the two points. Then when you have that you know how find the velocity. I dont know. Are you answering my question? I am currently looking at this pdf as possible answer for me. Do you know of something similar but more simple? The download was free. https://www.researchgate.net/publication/236012521_Lambert_Universal_Variable_Algorithm
 
Last edited:

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,651
Reaction score
785
Points
128
I believe we are almost on a same page. I use the universal variable solution to solve a Lambert's problem in IMFD. It's a good method (easy to implement), But I don't understand any of it. (Why it works the way it does). The solution I posted is NOT better just easier to understand.

The Time of Flight is easy to compute. But if you already have a Time of Flight and you want to convert it to something else. Then the problem is that mankind does not posses a math to do that directly.

As an example the Mean Anomaly (i.e. Time) can be computed from Eccentric Anomaly:

\[ M = E - e\sin(E) \]

But if you already have the Mean Anomaly "M" (i.e. Time) and you want to solve eccentric anomaly "E". Then the math is unknown to a mankind. So, to solve it, one must guess the value of "E" that gives the "M" you already got. It's not blind guessing, there are methods like Newtons method or Bracket method that can be used:

This code uses a newtons method to solve "E" from "M"

C:
double mna2eca(double mna, double ecc)
{
    // iterative calculation of eccentric anomaly from mean anomaly
    register double eca, m, x;
    int i;
    if (ecc<1.0) {
        eca = mna + ecc*sin(mna);
          m = mna - eca+ecc*sin(eca);
        for (i=0; (fabs(m)>1e-14 && i<22); i++) {
            x = m / (1.0-ecc*cos(eca));
            if (x>1.0) x=1.0; else if (x<-1.0) x=-1.0;
            eca += x;
            m = mna - eca+ecc*sin(eca);
        }
        return eca;
    }
    else {
        eca=0;
        m = mna - ecc * sinh(eca) + eca;
        for (i = 0; fabs(m)>1e-14 && i<22; i++) {   
            x = m / (ecc * cosh(eca) - 1.0);
            if (x>1.0) x=1.0; else if (x<-1.0) x=-1.0;
            eca += x;
            m = mna - ecc * sinh(eca) + eca;
        }
    }
    return eca;
}

So, the equations in my previous post requires you to guess the value of true anomaly so that it matches the Time of Flight you got. Also the univesal variable solution has a "For" loop where the "guessing" work is done. In the "BMW" the universal variable solution is on page 231.

.
 
Last edited:

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,651
Reaction score
785
Points
128
I resurrected the code for the Lambert solver I mentioned earlier with some improvements. It has some stability issues in bad cases but should be easier to understand than some other solvers.

C:
VECTOR3 LambertXP(VECTOR3 p1, VECTOR3 p2, double time, double mu, double dm)
{
    double r1 = length(p1);
    double r2 = length(p2);
    double b = angle(p1, p2);
    double dT = 1.0;
    double tra = 0;
    double e, p, a, B1, B2;

    if (dm < 0) b = PI2 - b; // It's the long way around. dm = direction of motion

    // Compute brackets
    double w = r2 * sin(b) / sqrt(r1*r1 + r2*r2 - 2.0 * r1*r2*cos(b));
    double k = r1*r2*(1.0 - cos(b));
    double L = r1 + r2;
    double M = r1*r2*(1.0 + cos(b));

    double b1 = acos(w);
    double b2 = acos(-w);
    double b3 = limit(-acos(w));
    double b4 = limit(-acos(-w));
    double n1 = k / (L + sqrt(2.0*M));
    double n2 = k / (L - sqrt(2.0*M));

    if (dm > 0) {
        B2 = acos((n1 - r1) / r1);
        B1 = b3 - PI2;
    }
    else {
        B1 = -acos((n1 - r1) / r1);
        B2 =  b2;
    }
    
    // Stop if maximum iteration count is reached or error gets withing acceptable limit
    int i;
    
    for (i = 0; i < 42; i++)
    {
        tra = (B1 + B2) * 0.5; // Take a sample from between the brackets

        double tra2 = limit(tra + b);
        
        // Compute eccentricity from two radius, angle in between and true anomaly
        e = (r2 - r1) / (r1*cos(tra) - r2*cos(tra2));

        // Parameter
        p = (r1 * (1.0 + e * cos(tra)));   

        // Semi-major axis
        a = p / (1.0 - e*e);

        double E1 = tra2eca(limit(tra), e);
        double E2 = tra2eca(tra2, e);
        double M1 = eca2mna(E1, e);
        double M2 = eca2mna(E2, e);

        if (M2 < M1) M2 += PI2;

        // Mean motion 's/rad'           
        double n = sqrt(abs(a*a*a) / mu);

        // Time Error
        dT = time - (M2 - M1) * n;
        
        if (abs(dT) < 1e-5) break; // We are done

        // Move a bracket...
        if (dT > 0) B1 = tra;
        else        B2 = tra;
    }

    // Angular momentum
    double h = sqrt(mu*p);

    // "f" & "g" functions
    double f = 1.0 - r2*(1.0 - cos(b)) / p;
    double g = r1*r2*sin(b) / h;

    // Return Velocity
    VECTOR3 V = p2 / g - p1 * (f / g);

    return V;
}
 

ncc1701d

Member
Joined
Aug 9, 2009
Messages
204
Reaction score
6
Points
18
I resurrected the code for the Lambert solver I mentioned earlier with some improvements. It has some stability issues in bad cases but should be easier to understand than some other solvers.

C:
VECTOR3 LambertXP(VECTOR3 p1, VECTOR3 p2, double time, double mu, double dm)
{
    double r1 = length(p1);
    double r2 = length(p2);
    double b = angle(p1, p2);
    double dT = 1.0;
    double tra = 0;
    double e, p, a, B1, B2;

    if (dm < 0) b = PI2 - b; // It's the long way around. dm = direction of motion

    // Compute brackets
    double w = r2 * sin(b) / sqrt(r1*r1 + r2*r2 - 2.0 * r1*r2*cos(b));
    double k = r1*r2*(1.0 - cos(b));
    double L = r1 + r2;
    double M = r1*r2*(1.0 + cos(b));

    double b1 = acos(w);
    double b2 = acos(-w);
    double b3 = limit(-acos(w));
    double b4 = limit(-acos(-w));
    double n1 = k / (L + sqrt(2.0*M));
    double n2 = k / (L - sqrt(2.0*M));

    if (dm > 0) {
        B2 = acos((n1 - r1) / r1);
        B1 = b3 - PI2;
    }
    else {
        B1 = -acos((n1 - r1) / r1);
        B2 =  b2;
    }
   
    // Stop if maximum iteration count is reached or error gets withing acceptable limit
    int i;
   
    for (i = 0; i < 42; i++)
    {
        tra = (B1 + B2) * 0.5; // Take a sample from between the brackets

        double tra2 = limit(tra + b);
       
        // Compute eccentricity from two radius, angle in between and true anomaly
        e = (r2 - r1) / (r1*cos(tra) - r2*cos(tra2));

        // Parameter
        p = (r1 * (1.0 + e * cos(tra)));  

        // Semi-major axis
        a = p / (1.0 - e*e);

        double E1 = tra2eca(limit(tra), e);
        double E2 = tra2eca(tra2, e);
        double M1 = eca2mna(E1, e);
        double M2 = eca2mna(E2, e);

        if (M2 < M1) M2 += PI2;

        // Mean motion 's/rad'          
        double n = sqrt(abs(a*a*a) / mu);

        // Time Error
        dT = time - (M2 - M1) * n;
       
        if (abs(dT) < 1e-5) break; // We are done

        // Move a bracket...
        if (dT > 0) B1 = tra;
        else        B2 = tra;
    }

    // Angular momentum
    double h = sqrt(mu*p);

    // "f" & "g" functions
    double f = 1.0 - r2*(1.0 - cos(b)) / p;
    double g = r1*r2*sin(b) / h;

    // Return Velocity
    VECTOR3 V = p2 / g - p1 * (f / g);

    return V;
}
Do you by chance have an example that tested ok for you using this? I like to know I wrote my program correctly. I have to convert this to a different programming language. What value did you use for mu? Thanks this is helpful.
 

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,651
Reaction score
785
Points
128
I tested it by comparing the results to two other solvers using p-iteration and universal variable solution. All three are exact matches.

C:
    if (Program == 6)
    {
        VECTOR3 rpos, tpos;
        OBJHANDLE hBody = hVessel->GetGravityRef();
        OBJHANDLE hTgt = oapiGetVesselByName("Tgt");
        VESSEL *vTgt = oapiGetVesselInterface(hTgt);

        hVessel->GetRelativePos(hBody, rpos);
        vTgt->GetRelativePos(hBody, tpos);

        double mu = oapiGetMass(hBody) * GC;       // GC = 6.67259e-11
        double T = 700.0;
   
        VECTOR3 v_xp = LambertXP(rpos, tpos, T, mu, 1.0);

        VECTOR3 v_UB, v_PN;
        LambertUB(rpos, tpos, T, 1.0, mu, NULL, 0, &v_UB);
        LambertPN(rpos, tpos, T, 1.0, mu, NULL, &v_PN);
   
        TextA(5, pos, "Angle ", angle(v_xp, v_UB)*DEG); pos += ld;
        Text(5, pos, "Length ", length(v_xp) - length(v_UB)); pos += ld;

        TextA(5, pos, "Angle ", angle(v_PN, v_UB)*DEG); pos += ld;
        Text(5, pos, "Length ", length(v_PN) - length(v_UB)); pos += ld;
    }

In the following scenario when targeting the "Tgt" named vessel with 700s transfer time the result should be: (remember to start "paused")

x = 5603.860
y = -0.0669872
z = 6965.300
 

Attachments

  • LAMBERT.scn.txt
    2.9 KB · Views: 1

ncc1701d

Member
Joined
Aug 9, 2009
Messages
204
Reaction score
6
Points
18
I tested it by comparing the results to two other solvers using p-iteration and universal variable solution. All three are exact matches.

C:
    if (Program == 6)
    {
        VECTOR3 rpos, tpos;
        OBJHANDLE hBody = hVessel->GetGravityRef();
        OBJHANDLE hTgt = oapiGetVesselByName("Tgt");
        VESSEL *vTgt = oapiGetVesselInterface(hTgt);

        hVessel->GetRelativePos(hBody, rpos);
        vTgt->GetRelativePos(hBody, tpos);

        double mu = oapiGetMass(hBody) * GC;       // GC = 6.67259e-11
        double T = 700.0;
  
        VECTOR3 v_xp = LambertXP(rpos, tpos, T, mu, 1.0);

        VECTOR3 v_UB, v_PN;
        LambertUB(rpos, tpos, T, 1.0, mu, NULL, 0, &v_UB);
        LambertPN(rpos, tpos, T, 1.0, mu, NULL, &v_PN);
  
        TextA(5, pos, "Angle ", angle(v_xp, v_UB)*DEG); pos += ld;
        Text(5, pos, "Length ", length(v_xp) - length(v_UB)); pos += ld;

        TextA(5, pos, "Angle ", angle(v_PN, v_UB)*DEG); pos += ld;
        Text(5, pos, "Length ", length(v_PN) - length(v_UB)); pos += ld;
    }

In the following scenario when targeting the "Tgt" named vessel with 700s transfer time the result should be: (remember to start "paused")

x = 5603.860
y = -0.0669872
z = 6965.300
Thanks I will check out your info and let you know how it goes.
 

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,651
Reaction score
785
Points
128
I noticed a couple of problems with the code so here's an updated one:

C:
VECTOR3 LambertXP(VECTOR3 p1, VECTOR3 p2, double time, double mu, double dm)
{
    double r1 = length(p1);
    double r2 = length(p2);
    double b = angle(p1, p2);
    double dT = 1.0;
    double tra = 0;
    double e, p, a, B1, B2;

    if (dm < 0) b = PI2 - b;

    // Compute brackets
    double w = r2 * sin(b) / sqrt(r1*r1 + r2*r2 - 2.0 * r1*r2*cos(b));
    double k = r1*r2*(1.0 - cos(b));
    double L = r1 + r2;
    double M = r1*r2*(1.0 + cos(b));

    double b3 = limit(-acos(w));
    double n1 = k / (L + sqrt(2.0*M));
    double n2 = k / (L - sqrt(2.0*M));

    if (r2 > r1) {
        if (dm > 0) {
            B1 = b3 - PI2;
            B2 = acos((n1 - r1) / r1);
        }
        else {
            B1 = -acos((n1 - r1) / r1);
            B2 =  acos((n2 - r1) / r1);
        }
    }
    else { 
        if (dm > 0) {
            B1 = acos((n1 - r1) / r1);
            B2 = b3;
        }
        else {
            B1 = acos((n2 - r1) / r1);
            B2 = PI2 - acos((n1 - r1) / r1);
        }
    }

    // Stop if maximum iteration count is reached or error gets withing acceptable limit
    //
    for (int i = 0; i < 42; i++)
    {
        tra = (B1 + B2) * 0.5; // Take a sample from between the brackets

        double tra2 = limit(tra + b);
      
        // Compute eccentricity from two radius, angle in between and true anomaly
        e = (r2 - r1) / (r1*cos(tra) - r2*cos(tra2));

        // Parameter
        p = (r1 * (1.0 + e * cos(tra))); 

        // Semi-major axis for the "guess"
        a = p / (1.0 - e*e);

        double E1 = tra2eca(limit(tra), e);
        double E2 = tra2eca(tra2, e);
        double M1 = eca2mna(E1, e);
        double M2 = eca2mna(E2, e);

        if (M2 < M1) M2 += PI2;

        // Mean motion 's/rad'         
        double n = sqrt(abs(a*a*a) / mu);

        // Time Error
        dT = time - (M2 - M1) * n;
      
        if (abs(dT)/time < 1e-9) break; // We are done

        // Move brackets...
        if (r2>r1) {
            if (dT > 0) B1 = tra;
            else        B2 = tra;
        }
        else {
            if (dT < 0) B1 = tra;
            else        B2 = tra;
        }
    }

    // Angular momentum
    double h = sqrt(mu*p);

    // "f" & "g" functions
    double f = 1.0 - r2*(1.0 - cos(b)) / p;
    double g = r1*r2*sin(b) / h;

    // Return Velocity
    VECTOR3 V = p2 / g - p1 * (f / g);

    return V;
}


Here's the universal variable code for cross reference:
C:
void CompCS(double *C, double *S, double z)
{
    if (z == 0) {
        *C = 1.0 / 2.0;
        *S = 1.0 / 6.0;
    }
    else
    if (z>10e-9) {
        double sqrt_z = sqrt(z);
        *C = (1.0 - cos(sqrt_z)) / z;
        *S = (sqrt_z - sin(sqrt_z)) / (z*sqrt_z);
    }
    else
    if (z<-10e-9) {
        double sqrt_z = sqrt(-z);
        *C = (1.0 - cosh(sqrt_z)) / z;
        *S = (sinh(sqrt_z) - sqrt_z) / (-z*sqrt_z);
    }
    else {
        *C = 1.0 / 2.0 - z / 24.0 + z*z / 720.0 - z*z*z / 40320.0 + z*z*z*z / 3628800.0;
        *S = 1.0 / 6.0 - z / 120.0 + z*z / 5040.0 - z*z*z / 362880.0 + z*z*z*z / 39916800.0;
    }
}


VECTOR3 LambertUB(VECTOR3 _Pos0, VECTOR3 _Pos1, double time, double mu, double dm)
{
    double y = 0.0, z = 0.0, x, dt, zl, zu, C, S;

    double th = angle(_Pos0, _Pos1);
    double r0 = length(_Pos0);
    double r1 = length(_Pos1);

    double sqrt_mu = sqrt(mu);
    double r01 = r0 + r1;

    double A = dm * sqrt(r0*r1*(1.0 + cos(th)));

    if (th == 0 || A == 0) return _V(0,0,0);

    zu = 4.0*PI*PI;
    zl = -4.0*PI*PI;

    for (int i = 0; i<42; i++)
    {
        z = (zl + zu) * 0.5;

        CompCS(&C, &S, z);

        double sqrt_c = sqrt(C);

        y = r01 - A * (1.0 - z*S) / sqrt_c;

        if (y<0) zl = z;
        else {
            double sqrt_y = sqrt(y);

            x = sqrt_y / sqrt_c;

            dt = (x*x*x * S + A * sqrt_y) / sqrt_mu;

            if ((fabs(dt - time) / time) < 10e-7) break;

            if (dt <= time) zl = z;
            else            zu = z;
        }
    }

    double f = 1.0 - y / r0;
    double g = 1.0 / (A*sqrt(y / mu));

    return _Pos1*g - _Pos0*(f*g);
}

The test run was a vessel/vessel intercept trajectory with constant 700s or 2000s flight time depending if it was the short or the long way around. I didn't detect any error during a couple of orbits.
 

Attachments

  • Testing.jpg
    Testing.jpg
    259.4 KB · Views: 5
Last edited:

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,651
Reaction score
785
Points
128
I noticed a problem that the TrA iteration method fails if r1 = r2, first I thought there must be easy solution to that special case. But there appears to be none, so, the only workaround I can think is not to allow r2 to be equal with r1. On Low Earth Orbit it appears that 1.0 meter difference between r1, r2 is sufficient for stable operation. Therefore it's recommended to add the following line to a beginning of the code:

Code:
if (abs((r1 / r2) - 1) < 5e-7) r2 = r1 * (1.0 + 5e-7);

Here's also a plot that shows how the Ecc and Time of Flight behaves as function of TrA. Extreme cases leading to high eccentricity transfer orbits exists near the brackets and NaN-space.
 

Attachments

  • TrA-Plot.png
    TrA-Plot.png
    237.8 KB · Views: 5

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,651
Reaction score
785
Points
128
what programming language are you using?
I use C and C++ with Microsoft Visual Studio. The "community" edition is free for non-commercial use.
To what language you are trying to translate it to ?
 
Top