Attachment 'exp01.c'
Download 1 // exp01.c
2 // > make
3 // printd exp01.dat
4
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <math.h>
8 #include "nrlmsise-00.h"
9
10 // #define TEST 1
11
12 #define A0 120.0 // km bottom altitude simulated
13 #define A1 600.0 // km top altitude simulated
14 #define ASP 100.0 // km space altitude
15
16 #define RGEO 42164170.0 // m geosynchronous radius
17 #define RE 6378137.0 // m equatorial radius
18 #define MU 3.98600448e14 // m3/s2 gravitational parameter
19 #define AM 5.0 // m2/kg thinsat sail ratio
20 #define AE 2.0 // m2/kg experiment sail ratio, edge on?
21
22 #define KM 1000.0 // m/km
23 #define HOUR 3600.0 // seconds per hour
24 #define DAY 86400.0 // seconds per day
25
26 int main() {
27 struct nrlmsise_output output;
28 struct nrlmsise_input input;
29 struct nrlmsise_flags flags;
30 struct ap_array aph;
31 int i;
32
33 double pi = 4.0*atan(1.0); // compute pi
34
35 // KHL loop
36 // time of day
37 double hour ;
38 double h0 = 0.0 ;
39 double h1 = 24.0 ;
40 double hs = 0.25 ;
41 double havg = hs/(h1-h0) ;
42
43 // altitude above equator
44 double alt ;
45 double a0 = A0 ;
46 double a1 = A1 ;
47 double as = 1.0 ;
48 double dalt = 0.01 ;
49
50 // f07
51 double f107[3] ;
52 f107[0] = 70.0 ; f107[1] = 120.0 ; f107[2] = 200.0 ;
53 int fi ;
54
55 double denm ; // density dalt km down
56 double dens ;
57 double denp ; // density dalt km up
58 // double temp ;
59 double laps ;
60 double tdec[3] ;
61
62 /* input values */
63 for (i=0;i<7;i++)
64 aph.a[i]=100;
65
66 for (i=0;i<24;i++)
67 flags.switches[i]=1;
68
69 input.doy = 79 ; // March 21 vernal equinox
70 input.year = 0 ; // without effect
71 input.g_lat = 0.0 ;
72 input.g_long = 0.0 ;
73 input.ap = 4.0 ; // magnetic index
74
75 // header
76
77 printf( "#f107" );
78 for( i=0 ; i < 3 ; i++) printf( "%6.0f%9.0f%9.0f ", f107[i],f107[i],f107[i] ) ;
79 printf( "\n#alt");
80 for( i=0 ; i < 3 ; i++) printf( " dens. exper. decay" ) ;
81 printf( "\n" );
82
83 // altitude loop
84
85 tdec[0] = 0.0 ; // decay time to atmosphere
86 tdec[1] = 0.0 ; // decay time to atmosphere (100 km )
87 tdec[2] = 0.0 ; // decay time to atmosphere (100 km )
88
89 for( alt =a0 ; alt <= a1 ; alt += as ) {
90 printf( "%4.0f", alt );
91
92 // decay from geo -------------
93 // starting conditions experimental orbit
94
95 double rp = RE + alt * KM ; // m perigee radius
96 double ra = RGEO ; // m apogee radius
97 double sem = 0.5*(ra+rp) ; // m semimajor axis
98 double ecc = 0.5*(ra-rp)/sem ; // eccentricity
99 double v0 = sqrt( MU*((0.5/rp)+(0.5/ra))); // m/s
100 double vp = (1.0+ecc)*v0 ; // m/s perigee velocity
101 double va = (1.0-ecc)*v0 ; // m/s apogee velocity
102
103 #ifdef TEST
104 printf( "\n");
105 printf( "rp =%14.4e\n", rp );
106 printf( "ra =%14.4e\n", ra );
107 printf( "sem =%14.4e\n", sem );
108 printf( "ecc =%14.4e\n", ecc );
109 printf( "v0 =%14.4e\n", v0 );
110 printf( "vp =%14.4e\n", vp );
111 printf( "va =%14.4e\n", va );
112 #endif
113
114 for( fi = 0 ; fi <= 2 ; fi++ ) {
115 input.f107A=f107[fi];
116 input.f107 =f107[fi];
117
118 // prepare to take average
119 denm = 0.0 ;
120 denp = 0.0 ;
121 dens = 0.0 ;
122
123 for( hour=h0 ; hour < h1 ; hour += hs ) {
124 input.sec = hour * HOUR ;
125 input.lst = hour ;
126
127 input.alt = alt-dalt;
128 gtd7(&input, &flags, &output);
129 denm += output.d[5] ;
130
131 input.alt = alt+dalt;
132 gtd7(&input, &flags, &output);
133 denp += output.d[5] ;
134
135 input.alt = alt;
136 gtd7(&input, &flags, &output);
137 dens += output.d[5] ;
138 }
139
140 denm *= havg ; // normalize average
141 denp *= havg ;
142 dens *= havg ;
143
144 laps = 2.0*dalt /(log(denm/denp)); // compute lapse rate, km
145
146 // ------------
147
148 // Compute decay times
149
150 double texp = 0.0 ; // experiment time from geo (42164 km )
151 double H = laps*KM ; // lapse rate in meters
152 double vpt = sqrt( MU*((2.0/rp)-(1.0/(rp+H)))); // m/s vp thresh
153
154 #ifdef TEST
155 printf( "\n\n");
156 printf( "dens=%14.4e\n", dens);
157 printf( "H =%14.4e\n", H );
158 printf( "vpt =%14.4e\n", vpt );
159 #endif
160 double vpx = vp ; ; // experimental vp
161
162 // experimental decay loop
163 while( vpx > vpt ) {
164 double rvv = rp*vpx*vpx ;
165 double sq2 = rp / ( 2.0*MU - rvv ) ;
166 double dte = 2.0*pi*MU*sq2*sqrt(sq2) ; // orbit time
167 texp += dte ; // add orbit time
168 double dvp = dens*AE*vpx*sqrt((pi*rp*H)/(1.0-MU/rvv));
169 vpx -= dvp ; // slow down
170 #ifdef TEST
171 printf( "dt=%14.4e t=%14.4e dv=%14.4e vpx=%14.4e\n",
172 dte, texp, dvp, vpx );
173 #endif
174 } // integration loop
175
176 // circular decay, tumbling --------------
177 // dr/dt = dens * AM/2 * sqrt( mu rp )
178 // 2.0 <- half the area: disk divided by two faces of sphere
179
180 tdec[fi] += 2.0 * KM * as / ( dens * AM * sqrt( MU * rp ) ) ;
181 printf( "%7.0e%9.2e%9.2e", dens, texp/DAY, tdec[fi]/DAY );
182
183 } // F107 loop
184 printf( "\n" );
185 }
186 exit(0);
187 }
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.