Skip to content

Commit

Permalink
Fixed bug in orbit initialization from TLE
Browse files Browse the repository at this point in the history
  • Loading branch information
@ericstoneking
ericstoneking committed Sep 19, 2024
1 parent 10a1b3b commit 4be580d
Show file tree
Hide file tree
Showing 8 changed files with 59 additions and 24 deletions.
2 changes: 1 addition & 1 deletion Demo/Inp_Region.txt
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
******************** Regions for 42 *******************
2 ! Number of Regions
1 ! Number of Regions
---------------------------------------------------------
TRUE ! Exists
"TAG" ! Name
Expand Down
4 changes: 2 additions & 2 deletions InOut/Inp_Sim.txt
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
<<<<<<<<<<<<<<<<< 42: The Mostly Harmless Simulator >>>>>>>>>>>>>>>>>
************************** Simulation Control **************************
FAST ! Time Mode (FAST, REAL, EXTERNAL, or NOS3)
10000.0 0.1 ! Sim Duration, Step Size [sec]
1000.0 0.1 ! Sim Duration, Step Size [sec]
10.0 ! File Output Interval [sec]
0 ! RNG Seed
TRUE ! Graphics Front End?
Expand All @@ -14,7 +14,7 @@ TRUE Orb_LEO.txt ! Input file name for Orb 0
TRUE 0 SC_Simple.txt ! Existence, RefOrb, Input file for SC 0
***************************** Environment *****************************
04 08 2024 ! Date (UTC) (Month, Day, Year)
12 00 0.00 ! Time (UTC) (Hr,Min,Sec)
00 00 0.00 ! Time (UTC) (Hr,Min,Sec)
37.0 ! Leap Seconds (sec)
USER ! F10.7, Ap (USER, NOMINAL or TWOSIGMA)
230.0 ! USER-provided F10.7
Expand Down
17 changes: 17 additions & 0 deletions InOut/Optics_Simple.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
<<<<<<<<<<<<<<<<<< Optics Inputs for 42 >>>>>>>>>>>>>>>>>>>>>>
Minimal Optical System ! Description
2 ! Number of Elements
******************* Element 0 *********************
0 0 0 ! S/C, Body, Node
APERTURE ! Element Type (APERTURE, DETECTOR, PLANAR, PARABOLIC, SPHERICAL, THINLENS)
0.0 0.0 1.0 ! Axis
0.1 ! Aperture Diameter, m
1.0 ! Focal Length (ignored for DETECTOR, PLANAR)
2.0 CONCAVE ! Radius of Curvature, CONCAVE or CONVEX (SPHERICAL only)
******************* Element 1 *********************
0 0 1 ! S/C, Body, Node
DETECTOR ! Element Type (APERTURE, DETECTOR, PLANAR, PARABOLIC, SPHERICAL, THINLENS)
0.0 0.0 1.0 ! Axis
0.01 ! Aperture Diameter, m
0.01 ! Focal Length (ignored for DETECTOR, PLANAR)
2.0 CONCAVE ! Radius of Curvature, CONCAVE or CONVEX (SPHERICAL only)
2 changes: 1 addition & 1 deletion Kit/Source/glkit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3248,7 +3248,7 @@ GLuint BuildShaderProgram(GLuint VtxShader, GLuint FragShader, const char *Name)
glLinkProgram(ShaderProgram);

glGetProgramiv(ShaderProgram,GL_LINK_STATUS,&Success);

if (!Success) printf("Error in %s ShaderProgram link\n",Name);
glGetProgramiv(ShaderProgram,GL_INFO_LOG_LENGTH,&LogSize);
if (LogSize > 0) {
Expand Down
6 changes: 6 additions & 0 deletions Kit/Source/iokit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,12 @@ int FileToString(const char *file_name, char **result_string,
printf("Error reading from file %s\n", file_name);
return -1;
}
if (ret > file_len) {
printf("Error: Number of characters read (%d) exceeds expected file size (%d) for file %s\n",
ret,(int) file_len,file_name);
return -1;
}
(*result_string)[ret] = '\0';

close(fd);

Expand Down
42 changes: 23 additions & 19 deletions Kit/Source/orbkit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,12 +406,13 @@ void RV2Eph(double time, double mu, double xr[3], double xv[3],
#undef TWOPI
}
/**********************************************************************/
#define TWOPI (6.283185307179586)
#define PI (3.141592653589793)
#define D2R (1.74532925199E-2)
void TLE2MeanEph(const char Line1[80], const char Line2[80], double JD,
double LeapSec, struct OrbitType *O)
{
#define EPS (1.0E-12)
#define TWOPI (6.283185307179586)
#define D2R (1.74532925199E-2)

char YearString[3];
char DOYstring[13];
Expand Down Expand Up @@ -454,36 +455,37 @@ void TLE2MeanEph(const char Line1[80], const char Line2[80], double JD,

strncpy(RAANstring,&Line2[17],9);
RAANstring[9] = 0;
O->RAAN = ((double) atof(RAANstring))*D2R;
O->RAAN0 = ((double) atof(RAANstring))*D2R;

strncpy(EccString,&Line2[26],7);
EccString[7] = 0;
O->ecc = ((double) atof(EccString))*1.0E-7;

strncpy(omgstring,&Line2[34],8);
omgstring[8] = 0;
O->ArgP = ((double) atof(omgstring))*D2R;
O->ArgP0 = ((double) atof(omgstring))*D2R;

strncpy(MeanAnomString,&Line2[43],8);
MeanAnomString[8] = 0;
O->MeanAnom = ((double) atof(MeanAnomString))*D2R;
O->MeanAnom0 = ((double) atof(MeanAnomString))*D2R;

strncpy(MeanMotionString,&Line2[52],11);
MeanMotionString[11] = 0;
O->MeanMotion = ((double) atof(MeanMotionString))*TWOPI/86400.0;
O->Period = TWOPI/(O->MeanMotion);

/* Time of Periapsis passage given in seconds since J2000 */
O->tp = O->Epoch - O->MeanAnom/(O->MeanMotion);
while ((O->tp-DynTime) < -(O->Period)) O->tp += O->Period;
while ((O->tp-DynTime) > O->Period ) O->tp -= O->Period;
O->tp = O->Epoch - O->MeanAnom0/(O->MeanMotion);
while ((DynTime - O->tp) > O->Period) O->tp += O->Period;
while ((DynTime - O->tp) < -(O->Period)) O->tp -= O->Period;

O->MeanSMA = pow(mu/(O->MeanMotion*O->MeanMotion),1.0/3.0);
O->SMA = O->MeanSMA;
O->alpha = 1.0/(O->SMA);
O->SLR = O->SMA*(1.0 - O->ecc*O->ecc);
O->rmin = O->SLR/(1.0 + O->ecc);

O->MeanAnom = O->MeanMotion*(DynTime - O->tp);
O->anom = MeanAnomToTrueAnom(O->MeanAnom,O->ecc);

/* Initialize J2 Drift Parameters (ref Markley and Crassidis, Ch. 10) */
Expand All @@ -492,25 +494,22 @@ void TLE2MeanEph(const char Line1[80], const char Line2[80], double JD,
Coef = 1.5*J2*Re*Re/(O->SLR*O->SLR)*O->MeanMotion;
O->RAANdot = -Coef*cos(O->inc);
O->ArgPdot = Coef*(2.0-2.5*sin(O->inc)*sin(O->inc));
O->RAAN0 = O->RAAN - O->RAANdot*(DynTime-O->Epoch);
O->ArgP0 = O->ArgP - O->ArgPdot*(DynTime-O->Epoch);
/* 10.122 */
O->MeanAnom0 = O->MeanAnom - O->MeanMotion*(DynTime-O->Epoch);
O->RAAN = O->RAAN0 + O->RAANdot*(DynTime-O->Epoch);
while(O->RAAN > PI) O->RAAN -= TWOPI;
while(O->RAAN < -PI) O->RAAN += TWOPI;
O->ArgP = O->ArgP0 + O->ArgPdot*(DynTime-O->Epoch);
while(O->ArgP > PI) O->ArgP -= TWOPI;
while(O->ArgP < -PI) O->ArgP += TWOPI;
/* 10.126 */
O->J2Rw2bya = J2*Re*Re/O->MeanSMA;
}
else {
O->RAANdot = 0.0;
O->ArgPdot = 0.0;
O->RAAN0 = O->RAAN;
O->ArgP0 = O->ArgP;
O->RAAN = O->RAAN0;
O->ArgP = O->ArgP0;
O->J2Rw2bya = 0.0;
/* 10.122 */
O->MeanAnom0 = O->MeanAnom - O->MeanMotion*(DynTime-O->Epoch);
}

#undef TWOPI
#undef D2R
}
/**********************************************************************/
/* Ref: Markley and Crassidis, 10.4.3 */
Expand All @@ -530,6 +529,8 @@ void MeanEph2RV(struct OrbitType *O, double DynTime)

/* 10.122 */
O->MeanAnom = O->MeanAnom0 + O->MeanMotion*(DynTime - O->Epoch);
while(O->MeanAnom > PI) O->MeanAnom -= TWOPI;
while(O->MeanAnom < -PI) O->MeanAnom += TWOPI;

O->anom = MeanAnomToTrueAnom(O->MeanAnom,O->ecc);

Expand Down Expand Up @@ -582,6 +583,9 @@ void MeanEph2RV(struct OrbitType *O, double DynTime)
}

}
#undef TWOPI
#undef PI
#undef D2R
/**********************************************************************/
/* TLEs use UTC. 42 orbits use TT. So LeapSec are needed. */
long LoadTleFromFile(const char *Path, const char *TleFileName,
Expand Down
2 changes: 1 addition & 1 deletion Source/42glut.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1391,7 +1391,7 @@ int HandoffToGui(int argc, char **argv)
}

/* Comment out when OpenGL installation is stable */
/* CheckOpenGLProperties(); */
/* CheckOpenGLProperties(); */

/* .. Dive into Event Loop */
TimerDuration = 0;
Expand Down
8 changes: 8 additions & 0 deletions Source/42report.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,7 @@ void Report(void)
static FILE *IllumFile;
//static FILE *ProjAreaFile;
static FILE *AccFile;
static FILE *GpsFile;
//static FILE *Kepfile;
//static FILE *EHfile;
static char First = TRUE;
Expand Down Expand Up @@ -263,6 +264,9 @@ void Report(void)
AlbedoFile = FileOpen(InOutPath,"Albedo.42","w");
IllumFile = FileOpen(InOutPath,"Illum.42","w");
}
if (SC[0].Ngps > 0) {
GpsFile = FileOpen(InOutPath,"Gps.42","w");
}
}

if (OutFlag) {
Expand Down Expand Up @@ -361,6 +365,10 @@ void Report(void)
fprintf(AccFile,"%le %le ",SC[0].Accel[i].TrueAcc,SC[0].Accel[i].MeasAcc);
fprintf(AccFile,"\n");
}
if (SC[0].Ngps > 0) {
fprintf(GpsFile,"%le %le %le\n",
SC[0].GPS[0].PosN[0],SC[0].GPS[0].PosN[1],SC[0].GPS[0].PosN[2]);
}
if (SC[0].Ncss > 0) {
for(i=0;i<SC[0].Ncss;i++) {
fprintf(IllumFile,"%le ",SC[0].CSS[i].Illum);
Expand Down

0 comments on commit 4be580d

Please sign in to comment.