attachment:r500c.c of ReEntry

   1 // r300b.c 
   2 //   2012 Feb 28  KHL
   3 //   reentry of a thinsat
   4 //
   5 #include "math.h"
   6 #include "stdio.h"
   7 
   8 // the ballistic parameter B
   9 
  10 #define   B          0.4   // kg/m^2
  11 #define   Y0    500000.0   // initial altitude
  12 
  13 #define   MU   3.9846E14   // m^3/s^2  gravitation
  14 #define   RE   6378000.0   // m
  15 #define   SDAY   86164.1   // s
  16 
  17 #define   ASTEP  10000.0   // m
  18 #define   ANUM        51
  19 #define   TSTEP    0.001   // computation timestep
  20 #define   PSTEP      999   // 1000 compute steps per print step
  21 
  22 //--------------------------------------------------------
  23 // cheezy atmosphere model
  24 
  25 double adens[ANUM] = { 
  26    1.225e+00,   //   0 km
  27    4.135e-01,   //  10 km
  28    8.891e-02,   //  20 km
  29    1.841e-02,   //  30 km
  30    3.996e-03,   //  40 km
  31    1.027e-03,   //  50 km
  32    3.097e-04,   //  60 km
  33    8.283e-05,   //  70 km
  34    1.570e-05,   //  80 km
  35    3.416e-06,   //  90 km
  36    5.604e-07,   // 100 km
  37    9.708e-08,   // 110 km
  38    2.222e-08,   // 120 km
  39    8.152e-09,   // 130 km
  40    3.831e-09,   // 140 km
  41    2.076e-09,   // 150 km
  42    1.233e-09,   // 160 km
  43    7.815e-10,   // 170 km
  44    5.194e-10,   // 180 km
  45    3.581e-10,   // 190 km
  46    2.541e-10,   // 200 km
  47    1.846e-10,   // 210 km
  48    1.367e-10,   // 220 km
  49    1.029e-10,   // 230 km
  50    7.858e-11,   // 240 km
  51    6.073e-11,   // 250 km
  52    4.742e-11,   // 260 km
  53    3.738e-11,   // 270 km
  54    2.971e-11,   // 280 km
  55    2.378e-11,   // 290 km
  56    1.916e-11,   // 300 km
  57    1.552e-11,   // 310 km
  58    1.264e-11,   // 320 km
  59    1.035e-11,   // 330 km
  60    8.503e-12,   // 340 km
  61    7.014e-12,   // 350 km
  62    5.805e-12,   // 360 km
  63    4.820e-12,   // 370 km
  64    4.013e-12,   // 380 km
  65    3.350e-12,   // 390 km
  66    2.803e-12,   // 400 km
  67    2.350e-12,   // 410 km
  68    1.975e-12,   // 420 km
  69    1.662e-12,   // 430 km
  70    1.402e-12,   // 440 km
  71    1.184e-12,   // 450 km
  72    1.002e-12,   // 460 km
  73    8.492e-13,   // 470 km
  74    7.208e-13,   // 480 km
  75    6.127e-13,   // 490 km
  76    5.215e-13 }; // 500 km
  77 
  78 double dens( double alt ) {
  79 
  80    double frac = alt / ASTEP ;
  81    int    indx = (int) frac  ;
  82 
  83    if( indx < 0      ) { indx =  0     ; }
  84    if( indx > ANUM-2 ) { indx = ANUM-2 ; }
  85 
  86    frac -= (double) indx ;
  87 
  88    return  adens[indx] * pow( adens[indx+1]/adens[indx], frac );
  89 }
  90 
  91 
  92 //--------------------------------------------------------
  93 int main() {
  94 
  95 // initialize, relative to ground
  96 
  97    double y   = Y0              ; // altitude
  98    double vy  = 0.0             ; // vertical velocity
  99    double ay                    ; // vertical acceleration (grav-centrif)
 100    double ax                    ; // horizontal acceleration (drag)
 101    double at                    ; // total acceleration m/s2
 102    double dn                    ; // density
 103    double dr                    ; // drag parameter
 104    double time = 0.0            ;
 105 
 106    double pi2 = 8.0*atan(1.0)   ; 
 107    double omg = pi2 / SDAY      ; // earth rotation angular velocity
 108 
 109    double yr  = y + RE          ; // radius from center of earth
 110    double vxr = sqrt( MU / yr ) ; // orbital velocity at altitude
 111    double vxe = omg * yr        ; // earth rotational velocity
 112    double vx  = vxr - vxe       ; // initial velocity wrt atmosphere
 113    double vt                    ; // total velocity
 114 
 115    int    pcnt = PSTEP          ;
 116 
 117    printf("#   time   altitude     density     velocity     acceler    Vy\n");
 118 
 119 // compute loop
 120    while( y > 0.0 ) {
 121 
 122       dn =  dens(y)                        ;
 123       dr =  dn / B                         ;
 124       ax = -dr*vx*vx                       ; // -drag
 125       ay =  dr*vy*vy + (vxr*vxr -MU/yr)/yr ; // +drag, +centrif, -grav
 126 
 127       if( pcnt++ == PSTEP ) {
 128          at = sqrt( ax*ax + ay*ay );
 129          vt = sqrt( vx*vx + vy*vy );
 130 
 131          printf( "%8.0f%13.3f%12.3e%12.4f%12.2e%12.2e\n",
 132                   time,    y,   dn,   vt,   at,  -vy   );
 133          pcnt = 0 ;
 134       }
 135   
 136       time += TSTEP     ;
 137 
 138       y   += TSTEP * vy ;
 139       yr   = y + RE     ;
 140       vy  += TSTEP * ay ;
 141 
 142       vx  += TSTEP * ax ;
 143       vxe  = omg * yr   ;
 144       vxr  = vxe + vx   ;
 145    }
 146    return  0 ;
 147 }
 148 //--------------------------------------------------------
 149