// r300b.c 
//   2012 Feb 28  KHL
//   reentry of a thinsat
//
#include "math.h"
#include "stdio.h"

// the ballistic parameter B

#define   B          0.4   // kg/m^2
#define   Y0    500000.0   // initial altitude

#define   MU   3.9846E14   // m^3/s^2  gravitation
#define   RE   6378000.0   // m
#define   SDAY   86164.1   // s

#define   ASTEP  10000.0   // m
#define   ANUM        51
#define   TSTEP    0.001   // computation timestep
#define   PSTEP      999   // 1000 compute steps per print step

//--------------------------------------------------------
// cheezy atmosphere model

double adens[ANUM] = { 
   1.225e+00,   //   0 km
   4.135e-01,   //  10 km
   8.891e-02,   //  20 km
   1.841e-02,   //  30 km
   3.996e-03,   //  40 km
   1.027e-03,   //  50 km
   3.097e-04,   //  60 km
   8.283e-05,   //  70 km
   1.570e-05,   //  80 km
   3.416e-06,   //  90 km
   5.604e-07,   // 100 km
   9.708e-08,   // 110 km
   2.222e-08,   // 120 km
   8.152e-09,   // 130 km
   3.831e-09,   // 140 km
   2.076e-09,   // 150 km
   1.233e-09,   // 160 km
   7.815e-10,   // 170 km
   5.194e-10,   // 180 km
   3.581e-10,   // 190 km
   2.541e-10,   // 200 km
   1.846e-10,   // 210 km
   1.367e-10,   // 220 km
   1.029e-10,   // 230 km
   7.858e-11,   // 240 km
   6.073e-11,   // 250 km
   4.742e-11,   // 260 km
   3.738e-11,   // 270 km
   2.971e-11,   // 280 km
   2.378e-11,   // 290 km
   1.916e-11,   // 300 km
   1.552e-11,   // 310 km
   1.264e-11,   // 320 km
   1.035e-11,   // 330 km
   8.503e-12,   // 340 km
   7.014e-12,   // 350 km
   5.805e-12,   // 360 km
   4.820e-12,   // 370 km
   4.013e-12,   // 380 km
   3.350e-12,   // 390 km
   2.803e-12,   // 400 km
   2.350e-12,   // 410 km
   1.975e-12,   // 420 km
   1.662e-12,   // 430 km
   1.402e-12,   // 440 km
   1.184e-12,   // 450 km
   1.002e-12,   // 460 km
   8.492e-13,   // 470 km
   7.208e-13,   // 480 km
   6.127e-13,   // 490 km
   5.215e-13 }; // 500 km

double dens( double alt ) {

   double frac = alt / ASTEP ;
   int    indx = (int) frac  ;

   if( indx < 0      ) { indx =  0     ; }
   if( indx > ANUM-2 ) { indx = ANUM-2 ; }

   frac -= (double) indx ;

   return  adens[indx] * pow( adens[indx+1]/adens[indx], frac );
}


//--------------------------------------------------------
int main() {

// initialize, relative to ground

   double y   = Y0              ; // altitude
   double vy  = 0.0             ; // vertical velocity
   double ay                    ; // vertical acceleration (grav-centrif)
   double ax                    ; // horizontal acceleration (drag)
   double at                    ; // total acceleration m/s2
   double dn                    ; // density
   double dr                    ; // drag parameter
   double time = 0.0            ;

   double pi2 = 8.0*atan(1.0)   ; 
   double omg = pi2 / SDAY      ; // earth rotation angular velocity

   double yr  = y + RE          ; // radius from center of earth
   double vxr = sqrt( MU / yr ) ; // orbital velocity at altitude
   double vxe = omg * yr        ; // earth rotational velocity
   double vx  = vxr - vxe       ; // initial velocity wrt atmosphere
   double vt                    ; // total velocity

   int    pcnt = PSTEP          ;

   printf("#   time   altitude     density     velocity     acceler    Vy\n");

// compute loop
   while( y > 0.0 ) {

      dn =  dens(y)                        ;
      dr =  dn / B                         ;
      ax = -dr*vx*vx                       ; // -drag
      ay =  dr*vy*vy + (vxr*vxr -MU/yr)/yr ; // +drag, +centrif, -grav

      if( pcnt++ == PSTEP ) {
         at = sqrt( ax*ax + ay*ay );
         vt = sqrt( vx*vx + vy*vy );

         printf( "%8.0f%13.3f%12.3e%12.4f%12.2e%12.2e\n",
                  time,    y,   dn,   vt,   at,  -vy   );
         pcnt = 0 ;
      }
  
      time += TSTEP     ;

      y   += TSTEP * vy ;
      yr   = y + RE     ;
      vy  += TSTEP * ay ;

      vx  += TSTEP * ax ;
      vxe  = omg * yr   ;
      vxr  = vxe + vx   ;
   }
   return  0 ;
}
//--------------------------------------------------------