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