attachment:oneday05.c of HighApogeeConstruction

   1 // oneday05.c  - compute insertion delta V from launch loop, space elevator
   2 // time in sday, distances in RE
   3 
   4 #define MU     398600.4418   // km3/s2
   5 #define RE     6378.137      // km  earth equatorial radius
   6 #define SDAY   86164.099     // stellar day seconds
   7 
   8 #define PER     2000.0       // km  perigee altitude
   9 #define STEP        96       // steps per orbit (897.5 seconds, 14.96 min)
  10 #define ITER        20       // number of iterations
  11 #define MSTEP        4       // larger dots
  12 
  13 #include <stdio.h>
  14 #include <math.h>
  15 
  16 // ignore apsidal precession;  that raises apogee a few kilometers,
  17 // but timing is almost identical
  18 
  19 int main() {
  20    double  pi      = 4.0*atan(1.0)  ;  // onepi
  21    double  pi2     = 2.0*pi         ;  // twopi
  22    double  perigee = PER+RE         ;  // km perigee radius
  23    double  omegae  = pi2 / SDAY     ;  // earth rotation angular velocity
  24    double  tstep   = SDAY/STEP      ;  // time steps
  25 
  26    // rgeo is geosynchronous radius,
  27    //   also semimajor axis of construction orbit
  28    double  rgeo    = pow( (MU/(omegae*omegae)), (1.0/3.0)) ; // km
  29    double  apogee  = 2*rgeo-perigee       ; // km
  30    double  eps     = 1.0 - perigee/rgeo   ; // eccentricity
  31 
  32    double  rgeoR   = rgeo/RE              ; // geo in earth radii
  33 
  34    // radius = r0 / 1+e cos θ 
  35    double  r0      = perigee*apogee/rgeo  ; // km
  36    double  rp      = r0/(1+eps)           ; // km
  37    double  ra      = r0/(1-eps)           ; // km
  38 
  39    double  v0      = sqrt(MU/r0 )         ; // km/s 
  40    double  vp      = v0*(1+eps)           ; // km/s
  41    double  va      = v0*(1-eps)           ; // km/s
  42 
  43    double  omegap  = vp/rp                ; // rad/s
  44    double  omegaa  = va/ra                ; // rad/s
  45 
  46    int     index            ; 
  47    int     STEP2 = STEP / 2 ;  //  invert E for large steps
  48 
  49    double  told  = 0.0      ;  // old time
  50    double  time  = 0.0      ;  // for first step
  51    double  theta = 0.0      ;  // for first step
  52    double  err_t = 0.0      ;  // iteration time error
  53 
  54    printf( "#   theta  time(sday) "   );
  55    printf( "     xe(re)      ye(re)"  );
  56    printf( "     xe4(re)     ye4(re) ");
  57    printf( "     xf(re)      yf(re)"  );
  58    printf( "     xf4(re)     yf4(re) ");
  59    printf( "     xg(re)      yg(re)"  );
  60    printf( "     xg4(re)     yg4(re)" );
  61    printf( "\n");
  62 
  63    for( index = 0 ; index < STEP ;  index++ ) {
  64 
  65       double radius = r0 / ( 1.0 + eps * cos( theta ) ) ;
  66       double vtan   = v0 * ( 1.0 + eps * cos( theta ) ) ; // tangential velocity
  67       double vrad   = v0 *         eps * sin( theta )   ; // radial velocity
  68 
  69       // calculate time and relative angle
  70       double E      = acos((eps+cos(theta))/(1.0 + eps*cos(theta)));
  71       if( theta > pi )  E = pi2-E        ; // invert for second half of orbit
  72 
  73       double M      = E-eps*sin(E)       ;
  74       double thetaM = theta - M          ; // Earth-Sat relative angle, radian
  75 
  76       double time   = M / omegae         ; // time, seconds
  77       double sday   = time / SDAY ;
  78       double rday   = sday * pi2 ;
  79 
  80       double radR   = radius/RE          ;
  81 
  82       double xe     =   radR*sin(thetaM) ; // Earth-Sat x distance, km
  83       double ye     =   radR*cos(thetaM) ; // Earth-Sat y distance, km
  84 
  85       double xf     =   radR*sin(theta)  ; // Earth-Sat x distance, km
  86       double yf     =  -radR*cos(theta)  ; // Earth-Sat y distance, km
  87 
  88       double xg     =  rgeoR*sin(rday)   ; // GEO radius
  89       double yg     = -rgeoR*cos(rday)   ; // GEO radius
  90 
  91       double xe4, ye4, xf4, yf4, xg4, yg4 ;
  92       if( index % MSTEP == 0 ) {
  93         xe4=xe; ye4=ye;
  94         xf4=xf; yf4=yf;
  95         xg4=xg; yg4=yg;
  96        } 
  97 
  98       // rotating frame plot
  99       printf( "%9.6f%12.6f%12.6f%12.6f%12.6f%12.6f",
 100               theta, sday, xe,   ye,   xe4,  ye4  );
 101 
 102       // fixed frame plot 
 103       printf( "%12.6f%12.6f%12.6f%12.6f", xf, yf, xf4, yf4 );
 104 
 105       // geo plot, van Allen belt plot
 106       printf( "%12.6f%12.6f%12.6f%12.6f", xg, yg, xg4, yg4 );
 107 
 108       printf( "\n" );
 109 
 110       // now, iterate towards next time step 
 111       time += tstep   ; // this is the target time
 112 
 113       theta += ( vtan*tstep/radius ) ; // should get us close
 114 
 115       int itn ;
 116       for( itn = ITER ; itn > 0; itn-- ) { 
 117          E = acos((eps+cos(theta))/(1.0 + eps*cos(theta)));
 118          if( theta > pi )  E = pi2-E    ; // invert for second half of orbit
 119          err_t  = time-(E-eps*sin(E))/omegae  ;
 120          theta += vtan*err_t/radius            ;
 121          // printf( "DEBUG%16.8f%16.8f%16.8f\n", theta, err_t, tnext ) ;
 122       }
 123    } 
 124    return (0);
 125 }