Programming Question Calculate TLE to Orbiter Elements

SolarLiner · Nov 30, 2013

Hi,
Since my (somewhat buried) call from my blog post didn't get any responses, I write my problem here, in hope for answers.

In the process of writing AOSP, I decided to implement the initialization of OrbElements (the class that handle the "elements" line of the scenario) using TLE. The thing is, the conversion does not properly convert into Orbiter-readable elements.

Here's an example, using This TLE file:

Code:

ISS (ZARYA)
    ELEMENTS 13531.1087724418 0.0001222 1.6527 3.9354 87.3022 17.9391 56642.5793151855

TIANGONG 1
    ELEMENTS 13246.368845725 0.0006051 2.7752 72.5622 1.8501 177.2008 56642.5152893066

SOYUZ-TMA 11M
    ELEMENTS 13542.8059613007 9.89E-05 1.6486 112.8496 126.4609 354.0773 56620.6011352539

ISS OBJECT A
    ELEMENTS 13505.3219610729 0.0005079 1.6524 4.4834 59.9147 79.1278 56642.4480285645

ISS OBJECT B
    ELEMENTS 13494.7038836642 0.0005227 1.6522 4.1416 61.3932 84.1956 56642.5113220215

ISS OBJECT C
    ELEMENTS 13504.5505380456 0.0005437 1.6518 4.4822 61.9227 76.8127 56642.4476318359

TECHEDSAT-3P
    ELEMENTS 13483.3133946145 0.000245 1.6485 14.8294 347.8341 12.31 56640.3591003418

Try one of these, and you'll see that the resulting is not quite what I expect ...

Can anyone help me on this one?

Blog post: http://www.orbiter-forum.com/blog.php?b=1342
OrbElements class: http://sourceforge.net/p/orbaosp/code/HEAD/tree/trunk/OrbElements.cs

blixel · Nov 30, 2013

Check the description for links.

SolarLiner · Nov 30, 2013

Too bad there is no source code for it ... I don't want to put extra programs to use with the AOSP dll, so I want to port the code into AOSP directly.

EDIT: Oh, in the rush I didn't looked at the image ... got an Instantanious Headache trying to understand it ! :lol:

Ajaja · Dec 1, 2013

SolarLiner said:
Too bad there is no source code for it ...

Source code of eq_to_ecl?

Code:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

double gm = 3.9860044e14;
double ob_ecl=.40927970966248346; //23.45 deg;

double gr(double x)
{
  x=x*57.295779513082320875;
  while ((x<0)||(x>=360.0))
  {
    if (x<0) {x+=360; continue;}
    if (x>=360.0) {x-=360; continue;}
  }
  return x;
}

double rg(double x)
{
  return x*0.01745329251994329577;
}


void calc(double tm,double e,double eq_i,double eq_RAAN,double eq_ap,double eq_ma,double r)
{
  double pi=3.1415926535897932;
  int tm1=int(tm/1000.0);
  double tm2=tm-tm1*1000.0;
  double mjd=tm2+51544+365*tm1+(tm1-1)/4;
  // double tm3=(tm1+tm2/365)/100.0;
  // double ob_ecl=rg(23.4392911111-0.0130041666*tm3-0.00000163888*tm3*tm3+0.0000005036111*tm3*tm3*tm3);
  // printf("%.8f\n",gr(ob_ecl));
  double a = 0.720e3 * pow(0.5e1, 0.1e1 / 0.3e1) * pow(pow(r, -0.2e1) * gm * pow(pi, -0.2e1), 0.1e1 / 0.3e1); //24*60*60
  //double a = pow(0.2e1, 0.2e1 / 0.3e1) * pow(0.21541e5, 0.2e1 / 0.3e1) * pow(pow(r, -0.2e1) * gm * pow(pi, -0.2e1), 0.1e1 / 0.3e1); //86164
  //double a = 0.1228940547021246e4 * pow(pow(r, -0.2e1) * gm * pow(pi, -0.2e1), 0.1e1 / 0.3e1); //86164.09
  //double a = 0.1228940547021246e4 * pow(pow(r, -0.2e1) * gm * pow(pi, -0.2e1), 0.1e1 / 0.3e1);
  double eq_ml = eq_ma + eq_RAAN + eq_ap;

  double t1 = cos(eq_i);
  double t2 = cos(ob_ecl);
  double t3 = t1*t2;
  double t4 = sin(eq_i);
  double t5 = cos(eq_RAAN);
  double t7 = sin(ob_ecl);
  double t8 = t4*t5*t7;
  double t11 = t5*t5;
  double t12 = t2*t2;
  double t14 = t1*t1;
  double t16 = t14*t12;
  double t19 = sqrt(-2.0*t3*t8+1.0-t11+t12*t11+t14*t11-t16-t16*t11);
  double ecl_i = atan2(t19,t3+t8);
  double t21 = sin(eq_RAAN);
  double t23 = 1/t19;
  double t25 = t2*t4;
  double ecl_RAAN = atan2(t4*t21*t23,(t25*t5-t1*t7)*t23);
  double t30 = sin(eq_ap);
  double t33 = cos(eq_ap);
  double t36 = t1*t5;
  double t37 = t30*t7;
  double ecl_ap = atan2((t4*t30*t2-t21*t33*t7-t36*t37)*t23,(t25*t33-t33*t7*t36+t37*t21)*t23);

  double ecl_ma = eq_ma;
  double ecl_lp = ecl_ap+ecl_RAAN;
  double ecl_ml = ecl_lp+ecl_ma;

  printf("ELEMENTS %.8f ",a);
  printf("%.8f ",e);
  printf("%.8f ",gr(ecl_i));
  printf("%.8f ",gr(ecl_RAAN));
  printf("%.8f ",gr(ecl_lp));
  printf("%.8f ",gr(ecl_ml));
  printf("%.8f\n",mjd);
}

int main(int argc, char *argv[])
{
  FILE *f;
  char c1[256];
  char c2[256];
  char c3[256];

  if (argc!=2)
  {
    puts("Use eq_to_ecl.exe <filename.tle>");
    return -1;
  }

  f=fopen(argv[1],"rt");
  if (!f)
  {
    puts("File error!");
    return -1;
  }

  while (fgets(c1,255,f)&&fgets(c2,255,f))
  {
    char temp[255];
    char *stopstring;
    if(c1[0]!='1')
    {
      if ((c2[0]!='1')||(!fgets(c3,255,f)||(c3[0]!='2')))
      {
        puts("Error in the TLE file!");
        return -1;
      }
    }
    else
    {
      strcpy(c3,c2);
      strcpy(c2,c1);
      strncpy(c1,c2+3,6);
      c1[6]='\n';
      c1[7]='\0';
    }
    printf("%s",c1);

    strncpy(temp,c2+18,14);
    temp[14]='\0';
    double tm=strtod(temp,&stopstring);

    strncpy(temp,c3+26,7);
    temp[7]='\0';
    double e=strtod(temp,&stopstring);
    e=e/10000000.0;

    strncpy(temp,c3+8,8);
    temp[8]='\0';
    double eq_i=strtod(temp,&stopstring);

    strncpy(temp,c3+17,8);
    temp[8]='\0';
    double eq_RAAN=strtod(temp,&stopstring);

    strncpy(temp,c3+34,8);
    temp[8]='\0';
    double eq_ap=strtod(temp,&stopstring);

    strncpy(temp,c3+43,8);
    temp[8]='\0';
    double eq_ma=strtod(temp,&stopstring);

    strncpy(temp,c3+52,10);
    temp[10]='\0';
    double r=strtod(temp,&stopstring);

    calc(tm,e,rg(eq_i),rg(eq_RAAN),rg(eq_ap),rg(eq_ma),r);
  }
  fclose(f);

}

Quick_Nick · Dec 1, 2013

Shouldn't be hard to convert from TLE.
I posted this last little section to the wiki five years ago: http://www.orbiterwiki.org/wiki/Orbital_elements#Orbiter_Elements

Ajaja · Dec 1, 2013

SolarLiner said:
EDIT: Oh, in the rush I didn't looked at the image ... got an Instantanious Headache trying to understand it !

Maple (http://www.maplesoft.com/) solution of equations:

Code:

> restart; with(codegen):
> eq1:=cos(ecl_i)=cos(eq_i)*cos(ob_ecl)+sin(eq_i)*cos(eq_RAAN)*sin(ob_ecl):
> eq2:=cos(eq_i)=cos(ecl_i)*cos(ob_ecl)-sin(ecl_i)*cos(ecl_RAAN)*sin(ob_ecl):
> eq3:=sin(ecl_i)*sin(ecl_RAAN)=sin(eq_i)*sin(eq_RAAN):
> eq4:=sin(ecl_i)*sin(ecl_ap)=sin(eq_i)*sin(eq_ap)*cos(ob_ecl)-sin(eq_RAAN)*cos(eq_ap)*sin(ob_ecl)-cos(eq_i)*cos(eq_RAAN)*sin(eq_ap)*sin(ob_ecl):
> eq5:=sin(eq_i)*sin(eq_ap)=sin(ecl_i)*sin(ecl_ap)*cos(ob_ecl)+sin(ecl_RAAN)*cos(ecl_ap)*sin(ob_ecl)+cos(ecl_i)*cos(ecl_RAAN)*sin(ecl_ap)*sin(ob_ecl):
> res:=solve({eq1,eq2,eq3,eq4,eq5},[ecl_i,ecl_RAAN,ecl_ap]):
> res_all:=allvalues(res):
> res_all_1:=simplify(res_all[1]):
> res_all_2:=simplify(res_all[2]):
> C(res_all_1[1],optimized):
      t1 = cos(eq_i);
      t2 = cos(ob_ecl);
      t3 = t1*t2;
      t4 = sin(eq_i);
      t5 = cos(eq_RAAN);
      t7 = sin(ob_ecl);
      t8 = t4*t5*t7;
      t11 = t5*t5;
      t12 = t2*t2;
      t14 = t1*t1;
      t16 = t14*t12;
      t19 = sqrt(1.0-2.0*t3*t8-t11+t12*t11+t14*t11-t16-t16*t11);
      ecl_i = atan2(t19,t3+t8);
      t21 = sin(eq_RAAN);
      t23 = 1/t19;
      ecl_RAAN = atan2(t4*t21*t23,(t2*t4*t5-t1*t7)*t23);
      t30 = sin(eq_ap);
      t33 = cos(eq_ap);
      t36 = t1*t5;
      ecl_ap = atan2(-(-t4*t30*t2+t21*t33*t7+t36*t30*t7)*t23,-(-t4*t33*t2-t30*
t21*t7+t7*t33*t36)*t23);

Quick_Nick · Dec 1, 2013

@SolarLiner: your first and third numbers for each ELEMENTS are wrong; others may be correct.
The first number should be semi-major axis, but you're giving date (and in TLE format). The third element should be inclination...and I'm not sure what you're giving.

Ajaja · Dec 1, 2013

@Quick_Nick
Don't forget about equatorial->ecliptical plane conversion. Orbiter uses ecliptical plane, TLE - equatorial.

Quick_Nick · Dec 1, 2013

Ajaja said:
@Quick_Nick
Don't forget about equatorial->ecliptical plane conversion. Orbiter uses ecliptical plane, TLE - equatorial.

I see. I was really unsure. I may been quite wrong five years ago.

Ignore that then. :lol: I guess that would make the longitude of ascending node and other parameters different as well.

SolarLiner · Dec 1, 2013

Quick_Nick said:
@SolarLiner: your first and third numbers for each ELEMENTS are wrong; others may be correct.
The first number should be semi-major axis, but you're giving date (and in TLE format). The third element should be inclination...and I'm not sure what you're giving.

So the actual calculations are wrong, the SMA variable is really put in the first element of the list.
The inclination was given in radians I believe, so I made the conversion, but after the post.

Ajaja said:
Source code of eq_to_ecl?
{snip}

Ohoh, thanks Ajaja !

---------- Post added at 13:28 ---------- Previous post was at 12:38 ----------

Hmmm ... Eccentricity and inclination is good now, thanks to Ajaja and the eq_to_ecl, but now the SMA fails:

Yep, 28 000 km !

---------- Post added at 16:20 ---------- Previous post was at 13:28 ----------

Hmm ... I don't even know what's wrong ...
I looked at Ajaja solutions (the source code and the equation solution), and my code is all good, but it still gives me this weird apoapsis ...

Code:

private static OrbElements tlecalc(double tm, double e, double eq_i, double eq_RAAN, double eq_ap, double eq_ma, double r)
{
    double pi = 3.1315926535897932;
    double gm = 3.9860044e14;
    double ob_ecl = 0.40927970966248346;

    int tm1 = (int)Math.Truncate(tm / 1000.0);
    double tm2 = tm - tm1 * 1000.0;
    double mjd = tm2 + 51544 + 365 * tm1 + (tm1 - 1) / 4;

    double a = 0.720e3 * Math.Pow(0.5, 0.1 / 0.3) * Math.Pow(Math.Pow(r, -0.2) * gm * Math.Pow(pi, -0.2), 0.1 / 0.3);

    double eq_ml = eq_ma + eq_RAAN + eq_ap;

    double t1 = Math.Cos(eq_i);
    double t2 = Math.Cos(ob_ecl);
    double t3 = t1 * t2;
    double t4 = Math.Sin(eq_i);
    double t5 = Math.Cos(eq_RAAN);
    double t7 = Math.Sin(ob_ecl);
    double t8 = t4 * t5 * t7;
    double t11 = t5 * t5;
    double t12 = t2 * t2;
    double t14 = t1 * t1;
    double t16 = t14 * t12;
    double t19 = Math.Sqrt(-2.0 * t3 * t8 + 1.0 - t11 + t12 * t11 + t14 * t11 - t16 - t16 * t11);
    double ecl_i = Math.Atan2(t19, t3 + t8);
    double t21 = Math.Sin(eq_RAAN);
    double t23 = 1 / t19;
    double t25 = t2 * t4;
    double ecl_RAAN = Math.Atan2(t4 * t21 * t23, (t25 * t5 - t1 * t7) * t23);
    double t30 = Math.Sin(eq_ap);
    double t33 = Math.Cos(eq_ap);
    double t36 = t1 * t5;
    double t37 = t30 * t7;
    double ecl_ap = Math.Atan2((t4 * t30 * t2 - t21 * t33 * t7 - t36 * t37) * t23, (t25 * t33 - t33 * t7 * t36 + t37 * t21) * t23);

    double ecl_ma = eq_ma;
    double ecl_lp = ecl_ap + ecl_RAAN;
    double ecl_ml = ecl_lp + ecl_ma;

    return new OrbElements(a, e, gr(ecl_i), gr(ecl_RAAN), gr(ecl_lp), gr(ecl_ml), mjd);
}

The scenario that results:

Code:

BEGIN_DESC
This scenario was generated using TLE data the 12/01/2013 at 16:19:36 UTC.

Source data: http://celestrak.com/NORAD/elements/stations.txt
END_DESC

BEGIN_ENVIRONMENT
  System Sol
END_ENVIRONMENT

BEGIN_FOCUS
  Ship ISS_(ZARYA)
END_FOCUS

BEGIN_SHIPS
ISS_(ZARYA):ProjectAlpha_ISS
  STATUS Orbiting Earth
  ELEMENTS 34784825.25756 0.0001315 28.20286 358.52304 94.29286 135.12346 56627.55148
  AROT 0 0 0
  VROT 0 0 0
END
TIANGONG_1:Mir
  STATUS Orbiting Earth
  ELEMENTS 34710788.89436 0.0005963 38.43018 93.55864 68.04069 279.11119 56627.60525
  AROT 0 0 0
  VROT 0 0 0
END
SOYUZ-TMA_11M:ShuttlePB
  STATUS Orbiting Earth
  ELEMENTS 34788044.14496 0.0000989 63.38109 126.06352 228.30697 222.38427 56604.60114
  AROT 0 0 0
  VROT 0 0 0
END
ISS_OBJECT_A:Carina
  STATUS Orbiting Earth
  ELEMENTS 34777948.15994 0.0005397 28.20566 358.91475 67.62799 159.71019 56627.48188
  AROT 0 0 0
  VROT 0 0 0
END
ISS_OBJECT_B:Carina
  STATUS Orbiting Earth
  ELEMENTS 34774928.86567 0.0005167 28.20443 358.90145 69.81155 148.18575 56627.47712
  AROT 0 0 0
  VROT 0 0 0
END
ISS_OBJECT_C:Carina
  STATUS Orbiting Earth
  ELEMENTS 34777720.82043 0.0005417 28.20578 358.93300 67.30829 143.93279 56627.47868
  AROT 0 0 0
  VROT 0 0 0
END
TECHEDSAT-3P:Carina
  STATUS Orbiting Earth
  ELEMENTS 34769347.33907 0.0002825 28.21068 2.10562 8.47195 101.17735 56627.07795
  AROT 0 0 0
  VROT 0 0 0
END
END_SHIPS

N_Molson · Dec 1, 2013

Yep, 28 000 km !

Hehe someone went a bit crazy on the last orbit boost burn :lol:

Ajaja · Dec 1, 2013

@SolarLiner

Code:

double a = 0.720e3 * pow(0.5e1, 0.1e1 / 0.3e1) * pow(pow(r, -0.2e1) * gm * pow(pi, -0.2e1), 0.1e1 / 0.3e1);

Code:

double a = 0.720e3 * Math.Pow(0.5, 0.1 / 0.3) * Math.Pow(Math.Pow(r, -0.2) * gm * Math.Pow(pi, -0.2), 0.1 / 0.3);

0.1e1 = 1 not 0.1
-0.2e1 = -2 not -0.2
...

This code was autogenerated in Maple too, so it's look little weird.

SolarLiner · Dec 1, 2013

Ohh right! Sorry for the burnt eyes there :lol:
Yeah, Maple looks weird, but is fine when you finally see clear in it ^^

SolarLiner · Dec 3, 2013

Ajaja said:
@SolarLiner
0.1e1 = 1 not 0.1
-0.2e1 = -2 not -0.2

Well, the change brought the ISS back ... to 6M km SMA.

Code:

ELEMENTS 13573947.27966 0.0001292 28.89719 341.83065 87.79923 88.25263 56629.60682

Code:

double a = 0.720e3 * Math.Pow(0.5e1, 0.1e1 / 0.3e1) * Math.Pow(Math.Pow(r, -0.2e1) * gm * Math.Pow(pi, -0.2e1), 0.1e1 / 0.3e1);

Code:

double pi = 3.1315926535897932;
double gm = 3.9860044e14;

Ajaja · Dec 3, 2013

SolarLiner said:
Code:

double pi = 3.1315926535897932;

???

kuddel · Dec 3, 2013

we do all know that [ame="http://www.youtube.com/watch?v=O-Y-ua3WBi4"]this[/ame] is the correct value of PI :lol:

NOMAD · Dec 4, 2013

A small program I had done long time ago. It converts the TLE format to a simple format that contains the classical orbital elements of a spacecraft.
Hers's the source code. I built it in VS2010, if you want to build it in other IDE, you need to modify some functions like "strcpy_s","strcat_s", etc.

To use it, just download the data that you want to convert and save it in a txt file. Open the text file in the program and choose the convert option, it will convert all the data automatically.

SolarLiner · Dec 4, 2013

So, according to kuddel, I'll change my pi value to 3.000000000000
:lol: so much mistakes :facepalm:

---------- Post added at 13:39 ---------- Previous post was at 08:24 ----------

Well, it changed nothing ...
Rahh, what is wrong with my code? brawawawawawawawaw (<- sound of fake crying with real mainly tears)

Ajaja · Dec 4, 2013

@SolarLiner

Code:

            tle2[5] =line2.Substring(53, 11);
            double r = double.Parse(tle2[5]);

53->52

SolarLiner · Dec 6, 2013

Ajaja said:

And all of that was due to mistyping some number on the keyboard. :facepalm:

Anyway, it works wonder now ! Thanks to you and your "clear vision" over my code. Sometimes you're so much used to see it that you can't even see the big mistakes you've made.

Again big thanks Ajaja ! :cheers:

---------- Post added at 23:05 ---------- Previous post was at 22:55 ----------

And accurate !

On my way to release the testing tool, for those who may find any interests in it.

Programming Question Calculate TLE to Orbiter Elements

It's necessary, TARS.

Donator

It's necessary, TARS.

Active member

Passed the Turing Test

Active member

Passed the Turing Test

Active member

Passed the Turing Test

It's necessary, TARS.

Addon Developer

Active member

It's necessary, TARS.

It's necessary, TARS.

Active member

Donator

Member

Attachments

It's necessary, TARS.

Active member

It's necessary, TARS.

Similar threads