// orb01.c  - compute insertion delta V from launch loop, space elevator


#define MUKM   398600.4418   // km3/s2
#define RE     6378.0        // km  earth equatorial radius
#define LOOP   80.0          // km  loop altitude
#define PI2    6.28318530718 
#define SDAT   86164.0905    // sidereal day seconds

#define RB      6600.0       // km  orbit sweep bottom
#define RT     50000.0       // km  orbit sweep top
#define RSTEP    100.0       // km  orbit sweep step

#define RAB    29700.0       // km  lowest ra search 
#define RAT    53000.0       // km  highest ra search

#define NLINES    50         // number of lines to print

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

int main() {
   double  omegae = PI2 / SDAT  ;  // earth rotation angular velocity
   double  rloop  = RE + LOOP   ;  // km launch loop radius
   double  rgeo   = pow( (MUKM/(omegae*omegae)), (1.0/3.0)) ;
   double  muom   = 2.0 * MUKM / ( omegae * omegae ) ; 
   int     iter ;

   double  rdest ;
   int lines = 0 ;

   for( rdest = RB ; rdest <= RT ; rdest += RSTEP ) {
      if( lines-- == 0 ) {
         printf( "#\n# KMradius  Eradius   DVloop   DVelev   RAelev    Vdest\n" ) ;
         lines = NLINES ;
      }
      double vdest = sqrt( MUKM / rdest ) ;    // destination orbit velocity
      double rre   = rdest/RE             ;    // destination earth units

      // launch loop transfer orbit
      double lva   = sqrt( 2.0*MUKM / ( rdest * ( 1.0 + rdest/rloop ) ) ) ;
      double ldv   = fabs( vdest - lva );

      // space elevator transfer orbit
      // set up binary search
      double rab   = RAB ;
      double rat   = RAT ;
      double sra         ;
      double srp         ;

      for( iter = 20 ; iter-- > 0 ; ) {
         sra = 0.5 * (rab + rat) ;
         srp = sra / ( ( muom / ( sra*sra*sra ) ) - 1.0 ) ;
         if ( srp < rdest ) { rab = sra ; }
         else               { rat = sra ; }
      }
      double svd   = sra*sra*omegae / rdest ;
      double sdv   = fabs( vdest - svd )    ;

      printf( "%10.0f%9.5f%9.5f%9.5f%9.2f%9.5f\n",
               rdest, rre, ldv, sdv, sra, vdest );
   } 
   return (0);
}