// oneday05.c  - compute insertion delta V from launch loop, space elevator
// time in sday, distances in RE

#define MU     398600.4418   // km3/s2
#define RE     6378.137      // km  earth equatorial radius
#define SDAY   86164.099     // stellar day seconds

#define PER     2000.0       // km  perigee altitude
#define STEP        96       // steps per orbit (897.5 seconds, 14.96 min)
#define ITER        20       // number of iterations
#define MSTEP        4       // larger dots

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

// ignore apsidal precession;  that raises apogee a few kilometers,
// but timing is almost identical

int main() {
   double  pi      = 4.0*atan(1.0)  ;  // onepi
   double  pi2     = 2.0*pi         ;  // twopi
   double  perigee = PER+RE         ;  // km perigee radius
   double  omegae  = pi2 / SDAY     ;  // earth rotation angular velocity
   double  tstep   = SDAY/STEP      ;  // time steps

   // rgeo is geosynchronous radius,
   //   also semimajor axis of construction orbit
   double  rgeo    = pow( (MU/(omegae*omegae)), (1.0/3.0)) ; // km
   double  apogee  = 2*rgeo-perigee       ; // km
   double  eps     = 1.0 - perigee/rgeo   ; // eccentricity

   double  rgeoR   = rgeo/RE              ; // geo in earth radii

   // radius = r0 / 1+e cos θ 
   double  r0      = perigee*apogee/rgeo  ; // km
   double  rp      = r0/(1+eps)           ; // km
   double  ra      = r0/(1-eps)           ; // km

   double  v0      = sqrt(MU/r0 )         ; // km/s 
   double  vp      = v0*(1+eps)           ; // km/s
   double  va      = v0*(1-eps)           ; // km/s

   double  omegap  = vp/rp                ; // rad/s
   double  omegaa  = va/ra                ; // rad/s

   int     index            ; 
   int     STEP2 = STEP / 2 ;  //  invert E for large steps

   double  told  = 0.0      ;  // old time
   double  time  = 0.0      ;  // for first step
   double  theta = 0.0      ;  // for first step
   double  err_t = 0.0      ;  // iteration time error

   printf( "#   theta  time(sday) "   );
   printf( "     xe(re)      ye(re)"  );
   printf( "     xe4(re)     ye4(re) ");
   printf( "     xf(re)      yf(re)"  );
   printf( "     xf4(re)     yf4(re) ");
   printf( "     xg(re)      yg(re)"  );
   printf( "     xg4(re)     yg4(re)" );
   printf( "\n");

   for( index = 0 ; index < STEP ;  index++ ) {

      double radius = r0 / ( 1.0 + eps * cos( theta ) ) ;
      double vtan   = v0 * ( 1.0 + eps * cos( theta ) ) ; // tangential velocity
      double vrad   = v0 *         eps * sin( theta )   ; // radial velocity

      // calculate time and relative angle
      double E      = acos((eps+cos(theta))/(1.0 + eps*cos(theta)));
      if( theta > pi )  E = pi2-E        ; // invert for second half of orbit

      double M      = E-eps*sin(E)       ;
      double thetaM = theta - M          ; // Earth-Sat relative angle, radian

      double time   = M / omegae         ; // time, seconds
      double sday   = time / SDAY ;
      double rday   = sday * pi2 ;

      double radR   = radius/RE          ;

      double xe     =   radR*sin(thetaM) ; // Earth-Sat x distance, km
      double ye     =   radR*cos(thetaM) ; // Earth-Sat y distance, km

      double xf     =   radR*sin(theta)  ; // Earth-Sat x distance, km
      double yf     =  -radR*cos(theta)  ; // Earth-Sat y distance, km

      double xg     =  rgeoR*sin(rday)   ; // GEO radius
      double yg     = -rgeoR*cos(rday)   ; // GEO radius

      double xe4, ye4, xf4, yf4, xg4, yg4 ;
      if( index % MSTEP == 0 ) {
        xe4=xe; ye4=ye;
        xf4=xf; yf4=yf;
        xg4=xg; yg4=yg;
       } 

      // rotating frame plot
      printf( "%9.6f%12.6f%12.6f%12.6f%12.6f%12.6f",
              theta, sday, xe,   ye,   xe4,  ye4  );

      // fixed frame plot 
      printf( "%12.6f%12.6f%12.6f%12.6f", xf, yf, xf4, yf4 );

      // geo plot, van Allen belt plot
      printf( "%12.6f%12.6f%12.6f%12.6f", xg, yg, xg4, yg4 );

      printf( "\n" );

      // now, iterate towards next time step 
      time += tstep   ; // this is the target time

      theta += ( vtan*tstep/radius ) ; // should get us close

      int itn ;
      for( itn = ITER ; itn > 0; itn-- ) { 
         E = acos((eps+cos(theta))/(1.0 + eps*cos(theta)));
         if( theta > pi )  E = pi2-E    ; // invert for second half of orbit
         err_t  = time-(E-eps*sin(E))/omegae  ;
         theta += vtan*err_t/radius            ;
         // printf( "DEBUG%16.8f%16.8f%16.8f\n", theta, err_t, tnext ) ;
      }
   } 
   return (0);
}