Attachment 'exp02.c'
Download 1 // exp02.c
2 // make
3 // produces exp02.dat file
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 double tdec[3] ;
54 tdec[0] = 0.0 ; tdec[1] = 0.0 ; tdec[2] = 0.0 ;
55 int fi ;
56
57 double denm ; // density dalt km down
58 double dens ;
59 double denp ; // density dalt km up
60 // double temp ;
61 double laps ;
62
63 /* input values */
64 for (i=0;i<7;i++)
65 aph.a[i]=100;
66
67 for (i=0;i<24;i++)
68 flags.switches[i]=1;
69
70 input.doy = 79 ; // March 21 vernal equinox
71 input.year = 0 ; // without effect
72 input.g_lat = 0.0 ;
73 input.g_long = 0.0 ;
74 input.ap = 4.0 ; // magnetic index
75
76 // header
77
78 printf( "#f107" );
79 for( i=0 ; i < 3 ; i++) printf( "%8.0f%9.0f%9.0f ",f107[i],f107[i],f107[i] ) ;
80 printf( "\n#alt");
81 for( i=0 ; i < 3 ; i++) printf( " exper. decay vdrop" ) ;
82 printf( "\n" );
83
84 // altitude loop
85
86 for( alt =a0 ; alt <= a1 ; alt += as ) {
87 printf( "%4.0f", alt );
88
89 // decay from geo -------------
90 // starting conditions experimental orbit
91
92 double rp = RE + alt * KM ; // m perigee radius
93 double ra = RGEO ; // m apogee radius
94 double sem = 0.5*(ra+rp) ; // m semimajor axis
95 double ecc = 0.5*(ra-rp)/sem ; // eccentricity
96 double v0 = sqrt( MU*((0.5/rp)+(0.5/ra))); // m/s
97 double vp = (1.0+ecc)*v0 ; // m/s perigee velocity
98 double va = (1.0-ecc)*v0 ; // m/s apogee velocity
99
100 #ifdef TEST
101 printf( "\n");
102 printf( "rp =%14.4e\n", rp );
103 printf( "ra =%14.4e\n", ra );
104 printf( "sem =%14.4e\n", sem );
105 printf( "ecc =%14.4e\n", ecc );
106 printf( "v0 =%14.4e\n", v0 );
107 printf( "vp =%14.4e\n", vp );
108 printf( "va =%14.4e\n", va );
109 #endif
110
111 for( fi = 0 ; fi <= 2 ; fi++ ) {
112 input.f107A=f107[fi];
113 input.f107 =f107[fi];
114
115 // prepare to take average
116 denm = 0.0 ;
117 denp = 0.0 ;
118 dens = 0.0 ;
119
120 for( hour=h0 ; hour < h1 ; hour += hs ) {
121 input.sec = hour * HOUR ;
122 input.lst = hour ;
123
124 input.alt = alt-dalt;
125 gtd7(&input, &flags, &output);
126 denm += output.d[5] ;
127
128 input.alt = alt+dalt;
129 gtd7(&input, &flags, &output);
130 denp += output.d[5] ;
131
132 input.alt = alt;
133 gtd7(&input, &flags, &output);
134 dens += output.d[5] ;
135 }
136
137 denm *= havg ; // normalize average
138 denp *= havg ;
139 dens *= havg ;
140
141 laps = 2.0*dalt /(log(denm/denp)); // compute lapse rate, km
142
143 // ------------
144
145 // Compute decay times
146
147 double texp = 0.0 ; // experiment time from geo (42164 km )
148 double H = laps*KM ; // lapse rate in meters
149 double vpt = sqrt( MU*((2.0/rp)-(1.0/(rp+H)))); // m/s vp thresh
150
151 #ifdef TEST
152 printf( "\n\n");
153 printf( "dens=%14.4e\n", dens);
154 printf( "H =%14.4e\n", H );
155 printf( "vpt =%14.4e\n", vpt );
156 #endif
157 double vpx = vp ; ; // experimental vp
158
159 // experimental decay loop
160
161 while( vpx > vpt ) {
162 double rvv = rp*vpx*vpx ;
163 double sq2 = rp / ( 2.0*MU - rvv ) ;
164 double dte = 2.0*pi*MU*sq2*sqrt(sq2) ; // orbit time
165 texp += dte ; // add orbit time
166 double dvp = dens*AE*vpx*sqrt((pi*rp*H)/(1.0-MU/rvv));
167 vpx -= dvp ; // slow down
168 #ifdef TEST
169 printf( "dt=%14.4e t=%14.4e dv=%14.4e vpx=%14.4e\n",
170 dte, texp, dvp, vpx );
171 #endif
172 }
173
174 // circular decay, tumbling --------------
175 // dr/dt = dens * AM/2 * sqrt( mu rp )
176 // dr = KM * as ;
177 // 2.0 <- half the area: disk divided by two faces of sphere
178
179 double vdec = 0.5 * dens * AM * sqrt( MU * rp ) ;
180 tdec[fi] += KM * as / vdec ;
181 printf( "%9.2e%9.2e%9.2e", texp/DAY, tdec[fi]/DAY, vdec );
182 }
183 printf( "\n" );
184 }
185 exit(0);
186 }
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.