Attachment 'drop02.c'
Download 1 // drop02.c
2 // > make
3 // produces drop02.dat files
4
5 #define NAME "drop02"
6
7 #include "nrlmsise-00.h"
8 #include <stdio.h>
9 #include <string.h>
10 #include <stdlib.h>
11 #include <math.h>
12
13 #define AMAG 4.0 // magnetic activity
14 #define A0 450.0 // km bottom altitude simulated
15 #define AISS 415.0 // km iss altitude
16 #define AFAST 50.0 // km fast descent altitude
17 #define VZ0 -0.636 // m/s beginning falling velocity
18 #define TSTEP 1.0 // s print time step
19 #define PSTEP 16 // compute steps per print step
20 #define TFAST 60 // s fast descent timestep
21 #define P0 340 // W/m2 nighttime IR illumination
22
23 #define RE 6378137.0 // m equatorial radius
24 #define MU 3.98600448e14 // m3/s2 gravitational parameter
25 #define BB 5.67E-8 // W/m2K4 black body radiation, em=0.5, A=2
26 #define KM 1000.0 // m/km
27
28
29 double f107[3] = { 70.0, 150.0, 250.0 } ; // solar activity
30 double aam[2] = { 2.5, 5.0 } ; // m2/kg
31 double aiss = AISS ; // km ISS altitude
32
33 char fname[32];
34 char os[1024] ; // output string
35 char as[128] ; // added string
36
37 int main() {
38
39 // double pi = 4.0*atan(1.0); // compute pi
40
41 int sf107 = sizeof( f107 ) / sizeof( f107[0] );
42 int sam = sizeof( aam ) / sizeof( aam[0] );
43
44 struct nrlmsise_output output ;
45 struct nrlmsise_input input ;
46 struct nrlmsise_flags flags ;
47 struct ap_array aph ;
48 int i ;
49
50 // input values - copied from nrlmsise-test, use kg/m^3 for density
51 for (i=0;i<7 ;i++) aph.a[i] = AMAG ; // magnetic activity
52 for (i=0;i<24;i++) flags.switches[i] = 1 ;
53
54 input.doy = 79 ; // March 21 vernal equinox
55 input.year = 0 ; // without effect
56 input.g_lat = 0.0 ; // equator
57 input.g_long = 0.0 ; // prime meridian
58 input.ap = 4.0 ; // magnetic index, moderate
59 input.sec = 0.0 ; // assume midnight
60 input.lst = 0.0 ; // assume midnight
61
62 // directions
63 // r radial (was z)
64 // t tangential (was x)
65 // vv, aa (was total)
66
67 double alt ; // km altitude
68 double r ; // m radius
69 double vv, vv2 ; // m/s total velocity, squared
70 double aa ; // m/s2 total drag deceleration
71 double at ; // m/s2 orbital tangential acceleration
72 double ar ; // m/s2 radial acceleration (+ upwards)
73 double ag ; // m/s2 gravitational acceleration
74 double pow ; // W/m2 deceleration power
75 double temp ; // K temperature
76 double t4 ; // K^4 fourth power of temperature
77 double dt = TSTEP/PSTEP ; // sec compute step
78 double dens ; // kg/m3 density
79 double vt ; // m/s orbital tangential velocity
80 double vr ; // m/s vertical velocity (neg down))
81 int if107, iam ; // loop through
82
83 double tmax, altmax, vtmax, vrmax, aamax, powmax, densmax;
84
85 for( if107=0 ; if107 < sf107 ; if107++ ) {
86 for( iam=0 ; iam < sam ; iam++ ) {
87
88 double am = aam[iam] ; // m2/kg area to mass ratio
89 double f107x = f107[if107] ; // solar activity
90 int pcnt = 0 ;
91 int pstep = PSTEP ;
92 double tempmax = 300 ; // K temperature Kelvin
93 double time = 0.0 ; // s time in seconds
94
95 sprintf(fname,"%s_%03d_%02d.sum", NAME, (int)(f107x), (int)(10*am));
96 FILE * sum = fopen( fname, "w" );
97 sprintf(fname,"%s_%03d_%02d.dat", NAME, (int)(f107x), (int)(10*am));
98 FILE * dat = fopen( fname, "w" );
99
100 // print header
101 sprintf( os, "# time altitude vt vr");
102 sprintf( as, " accel. Power temp dens\n");
103 strcat( os, as );
104 sprintf( as, "# sec km m/s m/s");
105 strcat( os, as );
106 sprintf( as, " m/s2 Watts/m2 K kg/m3\n");
107 strcat( os, as );
108
109 fprintf( dat, "%s", os ); // to data file
110 fprintf( sum, "%s", os ); // to summary file
111
112 // initial orbit
113
114 #define KV1 0.6667 // no idea why
115
116 alt = A0 ; // km starting altitude
117 r = KM*alt + RE ; // m radius
118
119 input.f107 = f107x ; // solar activity
120 input.f107A = f107x ; // solar activity
121 input.alt = alt ; // km altitude
122
123 gtd7(&input, &flags, &output) ; // operate nrlmsise model
124 dens = output.d[5] ; // kg/m3 density
125
126 vt = sqrt(MU/r) ; // m/s orbital tangential velocity
127 vr = -KV1*dens*am*sqrt(MU*r) ; // m/s radial vertical velocity
128
129 while( alt > 0.0 ) {
130 input.alt = alt ; // km altitude
131
132 gtd7(&input, &flags, &output) ; // nrlmsise model
133
134 r = RE + alt*KM ; // m perigee radius
135 dens = output.d[5] ; // kg/m3 density
136 ag = vt*vt/r ; // m/s centripedal acceleration
137 ag -= MU/(r*r) ; // m/s2 gravitational acceleration
138 vv2 = vt*vt+vr*vr ; //
139 vv = sqrt(vv2) ; // m/s total velocity
140 pow = dens*vv2*vv ; // W/m2 drag power per area
141 pow += P0 ; // W/m2 add night earth IR temperature
142 t4 = pow/BB ; // K^4 fourth power of temperature
143 temp = sqrt(sqrt(t4)) ; // K black body temperature
144 aa = -dens*am*vv2 ; // m/s2 drag deceleration
145 at = aa*(vt/vv) ; // m/s2 horizontal deceleration
146 at -= 2.0*vr*vt/r ; // m/s2 vr negative, forward acceleration
147 ar = ag+aa*(vr/vv) ; // m/s2 acceleration
148
149 if( alt < aiss ) {
150 // to stdout and summary file
151 fprintf( dat, "#%8.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n",
152 time,alt, vt, vr, aa, pow, temp, dens );
153 fprintf( sum, "#%8.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n",
154 time,alt, vt, vr, aa, pow, temp, dens );
155 aiss = -99.9 ;
156 }
157
158 if( pcnt-- == 0 ) {
159 // to stdout only
160 fprintf( dat, "%9.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n",
161 time,alt, vt, vr, aa, pow, temp, dens );
162
163 if( ( alt < AFAST ) && ( fmod( time, 60.0) < 0.05 ) )
164 pstep = PSTEP * TFAST ;
165 pcnt += pstep ;
166 }
167 if( temp > tempmax ) {
168 tmax = time ;
169 altmax = alt ;
170 vtmax = vt ;
171 vrmax = vr ;
172 aamax = aa ;
173 powmax = pow ;
174 tempmax = temp ;
175 densmax = dens ;
176 }
177
178 alt += dt*vr/KM ; // km vr negative, altitude reduced
179 vr += dt*ar ; // m/s ar negative, velocity down increase
180 vt += dt*at ; // m/s at negative, orbit velocity decrease
181 time+= dt ; // sec time in seconds
182 }
183 // to stdout and summary file
184
185 sprintf( os, "#%8.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n",
186 time,alt, vt, vr, aa, pow, temp, dens );
187 sprintf( as, "#%8.1f%9.3f%9.1f%9.3f%9.3f%11.1f%9.1f%12.1e\n",
188 tmax,altmax,vtmax,vrmax,aamax,powmax,tempmax,densmax );
189 strcat( os, as );
190 sprintf( as, "# f107=%6.1f\n", f107x );
191 strcat( os, as );
192 sprintf( as, "# Area/Mass=%6.3f m2/kg\n", am );
193 strcat( os, as );
194 sprintf( as, "# start altitude=%9.3f km\n", A0 );
195 strcat( os, as );
196 fprintf( dat, "%s", os );
197 fprintf( sum, "%s", os );
198 fclose( dat );
199 fclose( sum );
200 }
201 }
202 exit(0);
203 }
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.