From 4be580d8defc968da97d937ce09ef107c07860cb Mon Sep 17 00:00:00 2001 From: Eric Stoneking Date: Thu, 19 Sep 2024 09:45:16 -0400 Subject: [PATCH] Fixed bug in orbit initialization from TLE --- Demo/Inp_Region.txt | 2 +- InOut/Inp_Sim.txt | 4 ++-- InOut/Optics_Simple.txt | 17 +++++++++++++++++ Kit/Source/glkit.c | 2 +- Kit/Source/iokit.c | 6 ++++++ Kit/Source/orbkit.c | 42 ++++++++++++++++++++++------------------- Source/42glut.c | 2 +- Source/42report.c | 8 ++++++++ 8 files changed, 59 insertions(+), 24 deletions(-) create mode 100644 InOut/Optics_Simple.txt diff --git a/Demo/Inp_Region.txt b/Demo/Inp_Region.txt index 83cd2353..607607b6 100644 --- a/Demo/Inp_Region.txt +++ b/Demo/Inp_Region.txt @@ -1,5 +1,5 @@ ******************** Regions for 42 ******************* -2 ! Number of Regions +1 ! Number of Regions --------------------------------------------------------- TRUE ! Exists "TAG" ! Name diff --git a/InOut/Inp_Sim.txt b/InOut/Inp_Sim.txt index ca1edbd6..5b5cb33e 100755 --- a/InOut/Inp_Sim.txt +++ b/InOut/Inp_Sim.txt @@ -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? @@ -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 diff --git a/InOut/Optics_Simple.txt b/InOut/Optics_Simple.txt new file mode 100644 index 00000000..edd2099d --- /dev/null +++ b/InOut/Optics_Simple.txt @@ -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) diff --git a/Kit/Source/glkit.c b/Kit/Source/glkit.c index d296d45f..19c4c1c8 100755 --- a/Kit/Source/glkit.c +++ b/Kit/Source/glkit.c @@ -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) { diff --git a/Kit/Source/iokit.c b/Kit/Source/iokit.c index e0e45c3e..ef53625c 100755 --- a/Kit/Source/iokit.c +++ b/Kit/Source/iokit.c @@ -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); diff --git a/Kit/Source/orbkit.c b/Kit/Source/orbkit.c index faa0a01c..93cab686 100755 --- a/Kit/Source/orbkit.c +++ b/Kit/Source/orbkit.c @@ -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]; @@ -454,7 +455,7 @@ 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; @@ -462,11 +463,11 @@ void TLE2MeanEph(const char Line1[80], const char Line2[80], double JD, 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; @@ -474,9 +475,9 @@ void TLE2MeanEph(const char Line1[80], const char Line2[80], double JD, 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; @@ -484,6 +485,7 @@ void TLE2MeanEph(const char Line1[80], const char Line2[80], double JD, 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) */ @@ -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 */ @@ -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); @@ -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, diff --git a/Source/42glut.c b/Source/42glut.c index 1a542a9d..08db4695 100755 --- a/Source/42glut.c +++ b/Source/42glut.c @@ -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; diff --git a/Source/42report.c b/Source/42report.c index ac5afeb8..95b8a1cf 100755 --- a/Source/42report.c +++ b/Source/42report.c @@ -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; @@ -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) { @@ -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