// a2.c  space elevator radar slab
//  cc -o a2 a2.c -lm
//  ./a2 > a2.txt

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

#define  RSEMAX  2400.0  // space elevator maximum
#define  PLAT      50.0  // platform altitude
#define  A0       200.0  // minimum space debris altitude
#define  A1       700.0  // Bradley minimum space debris altitude
#define  A2      1700.0  // Bradley maximum space debris altitude
#define  A3      2200.0  // maximum space debris altitude
#define  REY0    -300.0  // minimum earth circle
#define  RE      6378.0  // Earth equatorial radius
#define  DP       100    // points per one side of arc
#define  MU    398600.44 // gravitational parameter

void orbnum( double a ) {
   double  y = RE+PLAT ;
   double  x = sqrt( (RE+a)*(RE+a) - y*y );
   double  v = sqrt( MU/(RE+a) );
   double  s = x / v ;
   double  m = s / 60.0 ; 
   double  b = atan2( x, y );
   double  c = b*(RE+a)*(RE+a)-0.5*x*y ;
   printf( "#||%7.0f ||%6.0f ||%8.2f ||%8.2f ||%7.2f ||%12.4e ||\n",
             a, x, v, s, m, c );
}

int  main() { 
   double   le  = acos((RE+REY0)/ RE    ) ; // longitude of earth circle
   double   la0 = acos((RE+REY0)/(RE+A0)) ; // longitude A0 horizon
   double   la1 = acos((RE+REY0)/(RE+A1)) ; // longitude A1 horizon
   double   la2 = acos((RE+REY0)/(RE+A2)) ; // longitude A2 horizon
   double   la3 = acos((RE+REY0)/(RE+A3)) ; // longitude A3 horizon
   double   r2d = 45.0/atan(1.0) ;

   double   x, y ;

   double   xa3 = (RE+A3)*sin( la3 )  ;
   int        i ;

   for( i=-DP ; i <=DP ; i++ ) {
      double f  = ((double) i )/((double) DP) ;

      // radar horizon line
      x = f * xa3 ;
      y = PLAT    ;
      printf( "%7.1f%7.1f", x, y );

      // earth radius
      x = RE * sin( le*f ) ;
      y = RE * cos( le*f ) - RE   ;
      printf( "%8.1f%7.1f", x, y );

      // A0
      x = (RE+A0)*sin( la0*f ) ;
      y = (RE+A0)*cos( la0*f ) - RE ;
      printf( "%8.1f%7.1f", x, y );

      // A1
      x = (RE+A1)*sin( la1*f ) ;
      y = (RE+A1)*cos( la1*f ) - RE ;
      printf( "%8.1f%7.1f", x, y );

      // A2
      x = (RE+A2)*sin( la2*f ) ;
      y = (RE+A2)*cos( la2*f ) - RE ;
      printf( "%8.1f%7.1f", x, y );

      // A3
      x = (RE+A3)*sin( la3*f ) ;
      y = (RE+A3)*cos( la3*f ) - RE ;
      printf( "%8.1f%7.1f", x, y );

      // space elevator
      y = PLAT + (RSEMAX-PLAT)*fabs(f);
      printf( " 0.0 %7.1f", y );

      printf( "\n");
   }

   printf( "#%8.2f%8.2f%8.2f%8.2f\n", r2d*la0, r2d*la1, r2d*la2, r2d*la3 );

   for( i=100 ; i<2300 ; i+=200 ) orbnum( (double) i ) ;
}