File indexing completed on 2024-04-21 03:48:42

 // SPDX-License-Identifier: LGPL-2.1-or-later
 //
 // SPDX-FileCopyrightText: 2012 Gerhard HOLTKAMP
 //
 
 /***************************************************************************
 * Calculate Spacecraft around other planets                                *
 *                                                                          *
 *                                                                          *
 * Open Source Code. License: GNU LGPL Version 2+                          *
 *                                                                          *
 * Author: Gerhard HOLTKAMP,        26-AUG-2012                              *
 ***************************************************************************/
 
 /*------------ include files and definitions -----------------------------*/
 #include "planetarySats.h"
 #include <fstream>
 #include <iostream>
 #include <cstdlib>
 #include <cstring>
 #include <cmath>
 #include <ctime>
 using namespace std;
 
 #include "astrolib.h"
 
 // ################ Planetary Sats Class ####################
 
 PlanetarySats::PlanetarySats()
 {
  plsatinit();
 }
 
 PlanetarySats::~PlanetarySats()
 {
 
 }
 
 double PlanetarySats::atan23 (double y, double x)
  {
   // redefine atan2 so that it doesn't crash when both x and y are 0 
   double result;
 
   if ((x == 0) && (y == 0)) result = 0;
   else result = atan2 (y, x);
 
   return result;
  }
 
 void PlanetarySats::plsatinit()
 {
  // initialize planetary sat data
   pls_moonflg = false;
   pls_day = 1;
   pls_month = 1;
   pls_year = 2012;
   pls_hour = 0;
   pls_minute = 0;
   pls_second = 0;
   pls_del_auto = 1;
   pls_step = 60.0;
   pls_delta_rt = 0.0;
   getTime();
   getMars();
 
   strcpy (pls_satelmfl, "./planetarysats.txt");
 
 }
 
 void PlanetarySats::getTime ()  // Get System Time and Date
 {
   time_t tt;
   int hh, mm, ss;
   double jd, hr;
 
   tt = time(NULL);
   jd = 40587.0 + tt/86400.0; // seconds since 1-JAN-1970
 
   jd = jd + pls_delta_rt / 24.0;
   caldat(jd, hh, mm, ss, hr);
   pls_year = ss;
   pls_month = mm;
   pls_day = hh;
 
   dms(hr, hh, mm, jd);
   pls_hour = hh;
   pls_minute = mm;
   pls_second = int(jd);
   if (pls_del_auto) pls_del_tdut = DefTdUt(pls_year);
   setMJD(pls_year, pls_month, pls_day, hh, mm, jd);
 };
 
 void PlanetarySats::setStepWidth(double s)
 {
   // set the step width (in seconds) used for calculations
   if (s < 0.01) pls_step = 0.01;
   else pls_step = s;
 }
 
 void PlanetarySats::setDeltaRT(double drt)
 {
   pls_delta_rt = drt;
 }
 
 void PlanetarySats::setDeltaTAI_UTC(double d)
 {
   // c is the difference between TAI and UTC according to the IERS
   // we have to add 32.184 sec to get to the difference TT - UT
   // which is used in the calculations here
 
   pls_del_tdut = d + 32.184;
   pls_del_auto = 0;
 }
 
 void PlanetarySats::setAutoTAI_UTC()
 {
   // set the difference between TAI and UTC according to the IERS
   pls_del_auto = true;
   pls_del_tdut = DefTdUt(pls_year);
 }
 
 void PlanetarySats::setMJD(int year, int month, int day, int hour, int min, double sec)
 {
     // set the (MJD-) time currently used for calculations to year, month, day, hour, min, sec
     double jd;
 
     pls_year = year;
     pls_month = month;
     pls_day = day;
     pls_hour = hour;
     pls_minute = min;
     pls_second = sec;
 
     jd = ddd(hour, min, sec);
     jd = mjd(day, month, year, jd);
 
     pls_time = jd;
 
     if (pls_del_auto) pls_del_tdut = DefTdUt(pls_year);
     
 }
 
 void PlanetarySats::getDatefromMJD(double mjd, int &year, int &month, int &day, int &hour, int &min, double &sec)
 {
     // convert times given in Modified Julian Date (MJD) into conventional date and time
 
     double magn;
 
     caldat((mjd), day, month, year, magn);
     dms (magn,hour,min,sec);
     if (sec>30.0) min++;
     if (min>59)
      {
       hour++;
       min=0;
      };
 }
 
 void PlanetarySats::setSatFile(char* fname)
 {
   strcpy (pls_satelmfl, fname);
   
 }
 
 void PlanetarySats::setStateVector(double mjd, double x, double y, double z, double vx, double vy, double vz)
 {
     pls_rep[0] = x;
     pls_rep[1] = y;
     pls_rep[2] = z;
     pls_vep[0] = vx;
     pls_vep[1] = vy;
     pls_vep[2] = vz;
 
     int year = 0;
     int month = 0;
     int day = 0;
     int hour = 0;
     int min = 0;
     double sec = 0;
     getDatefromMJD(mjd, year, month, day, hour, min, sec);
     setMJD(year, month, day, hour, min, sec);
     pls_tepoch = pls_time;
     //pls_tepoch = pls_time + pls_del_tdut / 86400.0;  // epoch in TT
 }
 
 int PlanetarySats::getStateVector(int nsat)
 {
  // read the state vector from the planetary sat file
  // nsat = number of eligible sat to select (1 if first or only sat)
  // RETURN number of eligible sat selected, 0 if no suitable sats of file problems
  
  int fsc, j, k, nst, utc;
  int yr, month, day, hour, min;
  double sec;
  bool searching;
  ifstream cfle;
  char satname[40];  // name of satellite
  char plntname[40]; // name of planet
  
  strcpy(satname, "");
  strcpy(plntname, "");
  nst = 0;
 
  searching = true;
  fsc = 0;
 
  cfle.open(pls_satelmfl, ios::in);
  if (!cfle)
  {
   searching = false;
   cfle.close();
  };
  if(searching)
  {
   while (searching)
   {
     fsc = 1;
 
     if (!cfle.getline(satname,40)) fsc = 0; 
     else
     {
      k = strlen(satname);
      if ((k>1) && (satname[0] == '#'))
      {
       for (j=1; j<k; ++j)
       {
        pls_satname[j-1] = satname[j]; 
        if(pls_satname[j-1] == '\n') pls_satname[j-1] = '\0';
       };
       pls_satname[k-1] = '\0';
      }
      else fsc = 0;
     };
   
     if (cfle.eof())
     {
      fsc = 0;
      searching = false;
     };
 
     if (fsc)
     {
      if (!cfle.getline(plntname, 40)) fsc = 0;
      else
      {
       k = strlen(plntname);
       if (k>0)
       {
        if(plntname[k-1] == '\n') plntname[k-1] = '\0';
        if((k>1) && (plntname[k-2] == '\r')) plntname[k-2] = '\0';
       };
      };
     };
 
     if (fsc)
     {
      cfle >> yr >> month >> day >> hour >> min >> sec >> utc;
      if (cfle.bad()) fsc = 0;
      if (cfle.eof())
      {
        fsc = 0;
        searching = false;
       };
     };
 
     if (fsc)
     {
      cfle >> pls_rep[0] >> pls_rep[1] >> pls_rep[2];
      if (cfle.bad()) fsc = 0;
      if (cfle.eof())
      {
        fsc = 0;
        searching = false;
       };
     };
 
     if (fsc)
     {
      cfle >> pls_vep[0] >> pls_vep[1] >> pls_vep[2];
      if (cfle.bad()) fsc = 0;
     };
 
     if (fsc)
     {
       if (strncmp(pls_plntname, plntname, 4) == 0) nst++;
       if (nst == nsat)
       {
         searching = false;
 
         setMJD(yr, month, day, hour, min, sec);
         pls_tepoch = pls_time;
         if (utc) pls_tepoch = pls_tepoch + pls_del_tdut/86400.0;  // epoch in TT
       };
     };
   };
   cfle.close(); 
  };
 
  if (fsc == 0) nst = 0;
 
  return nst;                
 }   
 
 void PlanetarySats::setPlanet(char* pname)
 {
   pls_moonflg = false;
   strcpy(pls_plntname, pname);
   if (strncmp("Mars", pname, 4) == 0) getMars();
   if (strncmp("Venus", pname, 4) == 0) getVenus();
   if (strncmp("Mercury", pname, 4) == 0) getMercury();
   if (strncmp("Moon", pname, 4) == 0) getMoon();
 }
     
 void PlanetarySats::stateToKepler()
 {
  // convert state vector (mean equatorial J2000.0) into planetary Kepler elements
  
  double t, dt, ag, gm, re, j2;
  double n, c, w, a, ecc, inc;
  Vec3 r1, v1;
  Mat3 mx;
  
  dt = (pls_tepoch - 51544.5) / 36525.0;
  gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
  re = pls_R0;
  j2 = pls_J2;
     
  // convert into planet equatorial reference frame
  if (pls_moonflg)
  {
    mx = mxidn();
    r1 = mxvct (mx, pls_rep);
    v1 = mxvct (mx, pls_vep);
 
  }
  else
  {
    ag = (pls_axl0 + pls_axl1 * dt) * M_PI / 180.0; 
    mx = zrot(ag + M_PI / 2.0);
    r1 = mxvct (mx, pls_rep);
    v1 = mxvct (mx, pls_vep);
    
    ag = (pls_axb0 + pls_axb1 * dt) * M_PI / 180.0; 
    mx = xrot(M_PI / 2.0 - ag);
    r1 = mxvct (mx, r1);
    v1 = mxvct (mx, v1);
  };
 
  v1 *= 86400.0;  // convert into km / day
  
  oscelm(gm, pls_tepoch, r1, v1, t, pls_m0, pls_a, pls_ecc, pls_ra, pls_per, pls_inc);
 
  // now the mean motion
  a = pls_a;
  ecc = pls_ecc;
  inc = pls_inc;
       
  //  preliminary n
  if (a == 0) a = 1.0; // just in case
  if (a < 0) a = -a;   // just in case
  n = sqrt (gm / (a*a*a));
 
  // correct for J2 term
  w = 1.0 - ecc*ecc;
  if (w > 0)
  {
   w = pow (w, -1.5);
   c = sin (inc*M_PI/180.0);
   n = n*(1.0 + 1.5*j2*re*re/(a*a) * w * (1.0 - 1.5*c*c));
  }
  else n = 1.0;   // do something to avoid a domain error
 
  n = n / (2.0*M_PI);
  if (n > 1000.0) n = 1000.0;  // avoid possible errors
 
  pls_n0 = n;
                                 
 }
     
 void PlanetarySats::getKeplerElements(double &perc, double &apoc, double &inc, double &ecc, double &ra, double &tano, double &m0, double &a, double &n0)
 {
   // get Kepler elements of orbit with regard to the planetary equator and prime meridian.
 
   double t, gm;
   Vec3 r1, v1;
   Mat3 mx;
 
 
   if (pls_moonflg)  // for the Moon we normally work in J2000. Now get it into planetary
   {
     gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
     
     mx = getSelenographic(pls_tepoch);
     r1 = mxvct (mx, pls_rep);
     v1 = mxvct (mx, pls_vep);
     v1 *= 86400.0;  // convert into km / day
  
     oscelm(gm, pls_tepoch, r1, v1, t, m0, a, ecc, ra, tano, inc);
 
     // now the mean motion
     if (a == 0) a = 1.0; // just in case
     if (a < 0) a = -a;   // just in case
     n0 = sqrt (gm / (a*a*a));
     n0 = n0 / (2.0*M_PI);
   }
   else
   {
     a = pls_a;
     n0 = pls_n0;
     m0 = pls_m0;
     tano = pls_per;
     ra = pls_ra;
     ecc = pls_ecc;
     inc = pls_inc;
   };
 
   perc = pls_a * (1.0 - pls_ecc) - pls_R0;
   apoc = pls_a * (1.0 + pls_ecc) - pls_R0;  
 
 }
 
 int PlanetarySats::selectSat(char* sname)
 {
   // select specified satellite
   // RETURN 1 if successful, 0 if no suitable satellite found
 
   int nst, res, sl; 
   bool searching;
 
   searching = true;
   nst = 1;
   sl = strlen(sname);
   
   while (searching)
   {
     res = getStateVector(nst);
     if (res)
     {
       if (strncmp(pls_satname, sname, sl) == 0) searching = false;
     }
     else searching = false;
     nst++; 
   }; 
   
   return res;
 }
     
 void PlanetarySats::getSatName(char* sname) const
 {
   strcpy (sname, pls_satname);
 }
     
 void PlanetarySats::currentPos()
 {
   getSatPos(pls_time);
 }
     
 void PlanetarySats::nextStep()
 {
   pls_time = pls_time + pls_step / 86400.0;
   getSatPos(pls_time);
 }
 
 double PlanetarySats::getLastMJD() const
 {
   return pls_time;
 }
 
 void PlanetarySats::getPlanetographic(double &lng, double &lat, double &height) const
 {
   // planetographic coordinates from current state vector
   
   lng = pls_lng;
   lat = pls_lat;
   height = pls_height;
     
 }
 
 void PlanetarySats::getFixedFrame(double &x, double &y, double &z, double &vx, double &vy, double &vz)
 {
   // last state vector coordinates in planetary fixed frame
 
   x = pls_r[0];
   y = pls_r[1];
   z = pls_r[2];
   vx = pls_v[0];
   vy = pls_v[1];
   vz = pls_v[2];
 }
 
 void PlanetarySats::getSatPos (double tutc)
 {
   // Get Position of Satellite at MJD-time t (UTC) 
 
   const double mp2 = 2.0*M_PI;
 
   double t, dt, m0, ran, aper, inc, a, ecc, n0, re;
   double f, e, c, s, k, sh, j2, gm, fac, b1, b2, b3;
   int j;
 
   Vec3 r1, v1, rg1, s2;
   Mat3 m1, m2;
 
   // prepare orbit calculation
 
   t = tutc + pls_del_tdut/86400.0;
   dt = t - pls_tepoch;
 
   ecc = pls_ecc;
   if (ecc >= 1.0) ecc = 0.999;  // to avoid crashes
   a = pls_a;
   n0 = mp2 * pls_n0;
   if (a < 1.0) a = 1.0;  // avoid possible crashes later on
   re = pls_R0;
   f = pls_flat;
   j2 = pls_J2;
   gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
 
   // get current orbit elements ready
   aper = (re / a) / (1.0 - ecc*ecc);
   aper = 1.5 * j2 * aper*aper * n0;
   m0 = pls_inc*M_PI/180.0;
   ran = -aper * cos(m0) * dt;
   m0 = sin (m0);
   aper = aper * (2.0 - 2.5*m0*m0) * dt;
   ran = pls_ra * M_PI/180.0 + ran;
   aper = pls_per * M_PI/180.0 + aper;
   m0 = pls_m0*M_PI/180.0 + n0*dt;
   inc = pls_inc * M_PI/180.0;
 
   // solve Kepler equation
   if (a < 1.0) a = 1.0;  // avoid possible crashes later on
   e = eccanom(m0, ecc);
   fac = sqrt(1.0 - ecc*ecc);
   c = cos(e);
   s = sin(e);
   r1.assign (a*(c - ecc), a*fac*s, 0.0);
 
   m0 = 1.0 - ecc*c;
   k = sqrt(gm/a);
   v1.assign (-k*s/m0, k*fac*c/m0, 0.0);
 
   // convert into reference plane 
   m1 = zrot (-aper);
   m2 = xrot (-inc);
   m1 *= m2;
   m2 = zrot (-ran);
   m2 = m2 * m1;
   r1 = mxvct (m2, r1);
   v1 = mxvct (m2, v1);
   v1 /= 86400.0;
 
   // save state vector in planet-fixed frame
 
   if (pls_moonflg) m1 = getSelenographic(t);
   else m1 = zrot((pls_W + pls_Wd*(t-51544.5))*(M_PI/180.0));
   pls_r = mxvct(m1,r1);
   pls_v = mxvct(m1,v1);
   pls_r *= 1000.0;
   pls_v *= 1000.0;
 
   // get the groundtrack coordinates
 
   rg1 = mxvct(m1,r1);
 
   s2 = carpol(rg1);
   pls_lat = s2[2];   // just preliminary
   pls_lng = s2[1];  // measured with the motion of rotation!
   if (pls_lng > mp2) pls_lng -= mp2;
   if (pls_lng < -M_PI) pls_lng += mp2;
   if (pls_lng > M_PI) pls_lng -= mp2;
 
   // get height above ellipsoid and geodetic latitude
   if (abs(r1) > 0.1)
   {
    if (f == 0) pls_height = abs(r1) - re;
    else
    {
     c = f*(2.0 - f);
     s = r1[0]*r1[0] + r1[1]*r1[1];
     sh = c*r1[2];
 
     for (j=0; j<4; ++j)
      {
       b1 = r1[2] + sh;
       b3 = sqrt(s+b1*b1);
       if (b3 < 1e-5) b3 = sin(pls_lat);  // just in case
       else b3 = b1 / b3;
       b2 = re / sqrt(1.0 - c*b3*b3);
       sh = b2 * c * b3;
      };
 
     sh = r1[2] + sh;
     pls_lat = atan20(sh,sqrt(s));
     sh = sqrt(s+sh*sh) - b2;
     pls_height = sh;
    }
   } 
   else pls_height = 0;  // this should never happen
   
   pls_lat = pls_lat * 180.0 / M_PI;
   pls_lng = pls_lng * 180.0 / M_PI;
 
  }
     
 
 void PlanetarySats::getMars()  // Mars planetary constants
 {
   pls_J2 = 1.964e-3;
   pls_R0 = 3397.2;
   pls_flat = 0.00647630;
   pls_axl0 = 317.681;
   pls_axl1 = -0.108;
   pls_axb0 = 52.886;
   pls_axb1 = -0.061;    
   pls_W = 176.868;
   pls_Wd = 350.8919830;
   pls_GM = 4.282828596416e+13; // 4.282837405582e+13
 }   
 
 void PlanetarySats::getVenus()  // Venus planetary constants
 { 
   pls_J2 = 0.027e-3;
   pls_R0 = 6051.9;
   pls_flat = 0.0;
   pls_axl0 = 272.72;
   pls_axl1 = 0.0;
   pls_axb0 = 67.15;
   pls_axb1 = 0.0;   
   pls_W = 160.26;
   pls_Wd = -1.4813596;
   pls_GM = 3.24858761e+14;  
 }   
 
 void PlanetarySats::getMercury()  // Mercury planetary constants
 { 
   pls_J2 = 0.0;
   pls_R0 = 2439.7;
   pls_flat = 0.0;
   pls_axl0 = 281.01;
   pls_axl1 = -0.033;
   pls_axb0 = 61.45;
   pls_axb1 = -0.005;    
   pls_W = 329.71;
   pls_Wd = 6.1385025;
   pls_GM = 2.20320802e+13;  
 }   
 
 void PlanetarySats::getMoon()  // Moon planetary constants
 {
   pls_moonflg = true; 
   pls_J2 = 0.2027e-3;
   pls_R0 = 1738.0;
   pls_flat = 0.0;
   pls_axl0 = 0.0;
   pls_axl1 = 0.0;
   pls_axb0 = 90.0;
   pls_axb1 = 0.0;   
   pls_W = 0.0;
   pls_Wd = 13.17635898;
   pls_GM = 4.90279412e+12;  
 }   
 
 Mat3 PlanetarySats::getSelenographic (double jd)
 {
   // Calculate the Matrix to transform from Mean of J2000 into selenographic
   // coordinates at MJD time jd.
   
   double t, gam, gmp, l, omg, mln;
   double a, b, c, ic, gn, gp, omp;
   double const degrad = M_PI / 180.0;
   Mat3 m1, m2;
 
   t = (jd - 15019.5) / 36525.0;
   gam = 281.2208333 + ((0.33333333e-5*t + 0.45277778e-3)*t + 1.7191750)*t;
   gmp = 334.3295556 + ((-0.125e-4*t - 0.10325e-1)*t + 4069.0340333)*t;
   l = 279.6966778 + (0.3025e-3*t + 36000.768925)*t;
   omg = 259.1832750 + ((0.22222222e-5*t + 0.20777778e-2)*t - 1934.1420083)*t;
   b = 23.45229444 + ((0.50277778e-6*t -0.16388889e-5)*t - 0.130125e-1)*t;
   mln = 270.4341639 + ((0.1888889e-5*t -0.11333e-2)*t + 481267.8831417)*t;
   ic = 1.535*degrad;
   gn = (l - gam)*degrad;
   gp = (mln - gmp)*degrad;
   omp = (gmp - omg)*degrad;
   a = -107.0*cos(gp) + 37.0*cos(gp+2.0*omp) -11.0*cos(2.0*(gp+omp));
   a = a*0.000277778*degrad + ic;
   c = (-109.0*sin(gp) + 37.0*sin(gp+2.0*omp) - 11.0*sin(2.0*(gp+omp)))/sin(ic);
   c = (c*0.000277778 + omg)*degrad;
   gn = -12.0*sin(gp) + 59.0*sin(gn) + 18.0*sin(2.0*omp);
   gn = gn*0.000277778*degrad;  // tau
 
   b *= degrad;
   gam = cos(a)*cos(b) + sin(a)*sin(b)*cos(c);
   gmp = gam*gam;
   if(gmp > 1.0) gmp = 0;
   else gmp = sqrt(1.0 - gmp);
   gam = atan23(gmp,gam);  // theta
   t = cos(a)*sin(b) - sin(a)*sin(b)*cos(c);
   l = -sin(a)*sin(c);
 
   gmp = atan23(l,t);  // phi
   t = sin(a)*cos(b) - cos(a)*sin(b)*cos(c);
   l = -sin(b)*sin(c);
   a = atan23(l,t);  // delta
   c = a + mln*degrad + gn - c;   // psi
 
   // libration rotation matrix from Mean equator to true selenographic
   m1 = zrot(gmp);
   m2 = xrot(gam);
   m1 = m2 * m1;
   m2 = zrot(c);
   m1 = m2 * m1;
 
   t = julcent(jd);
   m2 = pmatequ(0,t);  // convert from mean of J2000 to mean of epoch
   m1 = m1 * m2;
 
   return m1;
 }