// drop02.c // > make // produces drop02.dat files #define NAME "drop02" #include "nrlmsise-00.h" #include #include #include #include #define AMAG 4.0 // magnetic activity #define A0 450.0 // km bottom altitude simulated #define AISS 415.0 // km iss altitude #define AFAST 50.0 // km fast descent altitude #define VZ0 -0.636 // m/s beginning falling velocity #define TSTEP 1.0 // s print time step #define PSTEP 16 // compute steps per print step #define TFAST 60 // s fast descent timestep #define P0 340 // W/m2 nighttime IR illumination #define RE 6378137.0 // m equatorial radius #define MU 3.98600448e14 // m3/s2 gravitational parameter #define BB 5.67E-8 // W/m2K4 black body radiation, em=0.5, A=2 #define KM 1000.0 // m/km double f107[3] = { 70.0, 150.0, 250.0 } ; // solar activity double aam[2] = { 2.5, 5.0 } ; // m2/kg double aiss = AISS ; // km ISS altitude char fname[32]; char os[1024] ; // output string char as[128] ; // added string int main() { // double pi = 4.0*atan(1.0); // compute pi int sf107 = sizeof( f107 ) / sizeof( f107[0] ); int sam = sizeof( aam ) / sizeof( aam[0] ); struct nrlmsise_output output ; struct nrlmsise_input input ; struct nrlmsise_flags flags ; struct ap_array aph ; int i ; // input values - copied from nrlmsise-test, use kg/m^3 for density for (i=0;i<7 ;i++) aph.a[i] = AMAG ; // magnetic activity for (i=0;i<24;i++) flags.switches[i] = 1 ; input.doy = 79 ; // March 21 vernal equinox input.year = 0 ; // without effect input.g_lat = 0.0 ; // equator input.g_long = 0.0 ; // prime meridian input.ap = 4.0 ; // magnetic index, moderate input.sec = 0.0 ; // assume midnight input.lst = 0.0 ; // assume midnight // directions // r radial (was z) // t tangential (was x) // vv, aa (was total) double alt ; // km altitude double r ; // m radius double vv, vv2 ; // m/s total velocity, squared double aa ; // m/s2 total drag deceleration double at ; // m/s2 orbital tangential acceleration double ar ; // m/s2 radial acceleration (+ upwards) double ag ; // m/s2 gravitational acceleration double pow ; // W/m2 deceleration power double temp ; // K temperature double t4 ; // K^4 fourth power of temperature double dt = TSTEP/PSTEP ; // sec compute step double dens ; // kg/m3 density double vt ; // m/s orbital tangential velocity double vr ; // m/s vertical velocity (neg down)) int if107, iam ; // loop through double tmax, altmax, vtmax, vrmax, aamax, powmax, densmax; for( if107=0 ; if107 < sf107 ; if107++ ) { for( iam=0 ; iam < sam ; iam++ ) { double am = aam[iam] ; // m2/kg area to mass ratio double f107x = f107[if107] ; // solar activity int pcnt = 0 ; int pstep = PSTEP ; double tempmax = 300 ; // K temperature Kelvin double time = 0.0 ; // s time in seconds sprintf(fname,"%s_%03d_%02d.sum", NAME, (int)(f107x), (int)(10*am)); FILE * sum = fopen( fname, "w" ); sprintf(fname,"%s_%03d_%02d.dat", NAME, (int)(f107x), (int)(10*am)); FILE * dat = fopen( fname, "w" ); // print header sprintf( os, "# time altitude vt vr"); sprintf( as, " accel. Power temp dens\n"); strcat( os, as ); sprintf( as, "# sec km m/s m/s"); strcat( os, as ); sprintf( as, " m/s2 Watts/m2 K kg/m3\n"); strcat( os, as ); fprintf( dat, "%s", os ); // to data file fprintf( sum, "%s", os ); // to summary file // initial orbit #define KV1 0.6667 // no idea why alt = A0 ; // km starting altitude r = KM*alt + RE ; // m radius input.f107 = f107x ; // solar activity input.f107A = f107x ; // solar activity input.alt = alt ; // km altitude gtd7(&input, &flags, &output) ; // operate nrlmsise model dens = output.d[5] ; // kg/m3 density vt = sqrt(MU/r) ; // m/s orbital tangential velocity vr = -KV1*dens*am*sqrt(MU*r) ; // m/s radial vertical velocity while( alt > 0.0 ) { input.alt = alt ; // km altitude gtd7(&input, &flags, &output) ; // nrlmsise model r = RE + alt*KM ; // m perigee radius dens = output.d[5] ; // kg/m3 density ag = vt*vt/r ; // m/s centripedal acceleration ag -= MU/(r*r) ; // m/s2 gravitational acceleration vv2 = vt*vt+vr*vr ; // vv = sqrt(vv2) ; // m/s total velocity pow = dens*vv2*vv ; // W/m2 drag power per area pow += P0 ; // W/m2 add night earth IR temperature t4 = pow/BB ; // K^4 fourth power of temperature temp = sqrt(sqrt(t4)) ; // K black body temperature aa = -dens*am*vv2 ; // m/s2 drag deceleration at = aa*(vt/vv) ; // m/s2 horizontal deceleration at -= 2.0*vr*vt/r ; // m/s2 vr negative, forward acceleration ar = ag+aa*(vr/vv) ; // m/s2 acceleration if( alt < aiss ) { // to stdout and summary file fprintf( dat, "#%8.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n", time,alt, vt, vr, aa, pow, temp, dens ); fprintf( sum, "#%8.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n", time,alt, vt, vr, aa, pow, temp, dens ); aiss = -99.9 ; } if( pcnt-- == 0 ) { // to stdout only fprintf( dat, "%9.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n", time,alt, vt, vr, aa, pow, temp, dens ); if( ( alt < AFAST ) && ( fmod( time, 60.0) < 0.05 ) ) pstep = PSTEP * TFAST ; pcnt += pstep ; } if( temp > tempmax ) { tmax = time ; altmax = alt ; vtmax = vt ; vrmax = vr ; aamax = aa ; powmax = pow ; tempmax = temp ; densmax = dens ; } alt += dt*vr/KM ; // km vr negative, altitude reduced vr += dt*ar ; // m/s ar negative, velocity down increase vt += dt*at ; // m/s at negative, orbit velocity decrease time+= dt ; // sec time in seconds } // to stdout and summary file sprintf( os, "#%8.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n", time,alt, vt, vr, aa, pow, temp, dens ); sprintf( as, "#%8.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n", tmax,altmax,vtmax,vrmax,aamax,powmax,tempmax,densmax ); strcat( os, as ); sprintf( as, "# f107=%6.1f\n", f107x ); strcat( os, as ); sprintf( as, "# Area/Mass=%6.3f m2/kg\n", am ); strcat( os, as ); sprintf( as, "# start altitude=%9.3f km\n", A0 ); strcat( os, as ); fprintf( dat, "%s", os ); fprintf( sum, "%s", os ); fclose( dat ); fclose( sum ); } } exit(0); }