Attachment 'nl02.c'
Download 1 // nl02.c night illumination
2 // random, full, partial, and non-illuminating orientations
3 // circular orbits - we will do elliptical later, as well as solar inclination
4 //
5 // cc -o nl02 nl02.c -lm
6 //
7 // This produces two graphs:
8 // 1) Illumination versus hour of the night, in W/m^4^
9 // multiply by m^2^ of thinsat
10 // 2) Power versus orbital angle
11 //
12 // It produces three numbers per orientation:
13 // 1) Peak equatorial illumination at 7pm
14 // 2) Average power around orbit
15 // 3)
16 //
17 // sunwards is 0=noon, orbit 0 degrees
18 //
19 // note: I tried to use a min() function, that demands a -std=c99 flag
20 // (which is very silly)
21 //
22 // NOTE - - INCLUDE FAR SIDE ILLUMINATION
23
24 // #define DEBUG 1
25
26 #define NAME "nl02"
27
28 #define HOURSTART 6.0 // Start of night
29 #define HOUREND 12.0 // Middle of night (use symmetry)
30 #define DHOUR 0.0625 // Time steps through the night
31 #define DORB 0.01 // Orbital angle degree steps
32 // #define DHOUR 1.0 // Time steps through the night
33 // #define DORB 5.0 // Orbital angle degree steps
34
35 #define ALBEDO 0.5 // Worst case thinsat albedo for night light
36 #define EFF 0.123 // efficiency of solar cell
37 #define FILL 0.7 // solar cell fill
38 #define SUNPOWER 1367.0 // W/m^2 Sun illumination
39 #define TW 1E12 // Terawatt
40 #define AXIAL 22.3 // Axial tilt degrees
41 #define FULLMOON 27E-3 // W/m^2 Full moon illumination
42
43 #define C 2.997E8 // Speed of light KM per sec
44 #define KM 1000.0 // 1000m/km
45
46 #define RE (6378.0*KM) // Earth equatorial radius kilometers
47 #define M288 (12789.0*KM) // Server Sky orbit radius kilometers
48 #define RORB (M288)
49
50 #define PI 3.1415926535897932
51 #define PI2 (PI/2.0)
52 #define RAD (PI/180.0)
53
54 //----------------------------------------------------------------
55 // full moon equivalents per terawatt
56
57 #define SATAREA (1E12/(FILL*EFF*SUNPOWER))
58 #define REFLECT (SATAREA*ALBEDO*SUNPOWER)
59 #define FET ((DORB/360.0)*(REFLECT/FULLMOON))
60 //-----------------------------------------------------------------
61 #include <math.h>
62 #include <string.h>
63 #include <stdio.h>
64 #include <stdlib.h>
65 //-----------------------------------------------------------------
66
67 int main() {
68
69 FILE *datfile ; //
70 char gnuplot[200] ; //
71 char filename[80] ; //
72
73 double hour ; //
74 double orbit ; //
75
76 double asoe = asin( RE / RORB ) ; // shadow angle of earth on orbit
77 double acoe = acos( RE / RORB ) ; // shadow angle of earth on orbit
78 double adorb = RAD*DORB ; // orbital increment angle
79
80 // first, compute power level ----------------------------------
81
82 double pmax = 90.0/DORB ;
83 double pmin = 90.0/DORB ;
84 double pzero = 90.0/DORB ;
85
86 sprintf( filename, "%spow.dat", NAME ); // power
87 datfile = fopen( filename, "wb" );
88
89 for( orbit=90.0+DORB/2.0 ; orbit < 180.0 ; orbit += DORB ) {
90 double aorbit = RAD*orbit ;
91 double xorbit = RORB*sin(aorbit);
92 double yorbit = -RORB*cos(aorbit);
93
94 double imax = 0.0 ;
95 double imin = 0.0 ;
96 double izero = 0.0 ;
97
98 if( xorbit > RE ) {
99 imax = 1.0 ;
100 imin = sin( aorbit );
101 izero = cos( atan2( yorbit, xorbit-RE ));
102 }
103
104 pmax += imax ;
105 pmin += imin ;
106 pzero += izero ;
107
108 fprintf( datfile, "%8.3f%9.5f%9.5f%9.5f\n",
109 orbit, imax, imin, izero );
110 }
111
112 fclose( datfile );
113
114 // normalize power
115
116 pmax *= DORB/180.0 ;
117 pmin *= DORB/180.0 ;
118 pzero *= DORB/180.0 ;
119
120 printf( "power max=%6.4f min=%6.4f zero=%6.4f\n",
121 pmax, pmin, pzero );
122
123 // Scaling factor for full moon equivalents.
124 // this provides backend scaling for the night
125 // time scaling integtrated segment of orbit.
126 //
127 // FET = SATAREA // for 1 terawatt of thinsats,
128 // = 1E12/( fillfactor
129 // * efficiency
130 // * sunpower )
131 // * ALBEDO * sunpower
132 // * DORB/360 // scaling for integration
133 // / FULLMOON // full moon power, small
134 //
135 // see the #define statements above
136
137 double femax = FET/pmax ;
138 double femin = FET/pmin ;
139 double fezero = FET/pzero ;
140 double ferand = FET/pzero ; // pzero smallest
141
142 printf( "%12.3e=SATAREA%12.3e=ALBEDO%12.3e=DORB360%12.3e=FULLMOON\n",
143 SATAREA, ALBEDO, DORB/360.0, FULLMOON );
144 printf( "%12.3e=femax%12.3e=femin%12.3e=fezero%12.3e=ferand\n",
145 femax, femin, fezero, ferand );
146
147 // second, compute night time illumination ----------------------
148
149 sprintf( filename, "%snlp.dat", NAME ); // power
150 datfile = fopen( filename, "wb" );
151
152 for( hour=HOURSTART ; hour <= HOUREND ; hour += DHOUR ) {
153 double ahour = RAD*15.0*hour ;
154 double xhour = RE*sin(ahour); // 0 radians towards sun
155 double yhour = -RE*cos(ahour);
156
157 double nlrand = 0.0;
158 double nlmax = 0.0;
159 double nlmin = 0.0;
160 double nlzero = 0.0;
161
162 double aostart = 0.5*adorb + ahour - acoe ;
163 double aostop = ahour + acoe ;
164
165 // integrate over visible thinsats
166
167 double aorbit ;
168
169 for( aorbit=aostart ; aorbit < aostop ; aorbit += adorb ) {
170 double xorbit = RORB*sin(aorbit); // 0 radians towards sun
171 double yorbit = -RORB*cos(aorbit);
172
173 double distx = xorbit-xhour ;
174 double disty = yorbit-yhour ;
175 double dist2 = distx*distx + disty*disty ;
176
177 double dist = sqrt( dist2 );
178 double adist = atan2( distx, -disty ); // 0 radians towards sun
179
180 // assume in shadow, make perpendicular (fake) for 0 light
181 double amax = RAD*90.0 ;
182 double amin = RAD*90.0 ;
183 double azero = RAD*90.0 ;
184 double irand = 0.00 ;
185
186 if( yorbit < 0.0 ) { // in front of earth, daylight sky
187 amax = 0.0 ;
188 amin = 0.0 ;
189 azero = 0.0 ;
190 irand = 0.25 ;
191
192 } else if( xorbit > RE ) { // westward of earth
193 amax = 0.0 ;
194 amin = aorbit-RAD*90.0 ;
195 azero = atan2( yorbit, xorbit-RE );
196 irand = 0.25 ;
197
198 } else if( xorbit < -RE ) { // eastward of earth
199 amax = 0.0 ;
200 amin = aorbit-RAD*270.0 ;
201 azero = atan2( -yorbit, -xorbit-RE );
202 irand = 0.25 ;
203 }
204
205 // thinsat illumination
206
207 double imax = cos( amax );
208 double imin = cos( amin );
209 double izero = cos( azero );
210
211 // distance attenuation, receive angle
212 // cos( ahour-adist -pi/2 ) = sin( ahour-adist )
213
214 // Lambert, transmit angle
215 double cadmax = -cos(adist - amax ) ;
216 double cadmin = -cos(adist - amin ) ;
217 double cadzero = -cos(adist - azero ) ;
218
219 // receive angle
220 double cp = cos(ahour-adist)/dist2 ;
221
222 if( cp > 0.0 ) {
223 // add increment, unnormalized
224 if( cadmax > 0.0 ) { nlmax += cp * imax * cadmax ; }
225 if( cadmin > 0.0 ) { nlmin += cp * imin * cadmin ; }
226 if( cadzero > 0.0 ) { nlzero += cp * izero * cadzero ; }
227 nlrand += cp * irand ;
228 }
229
230 #ifdef DEBUG
231 printf( "X%9.4f=hour %9.4f=hstart %9.4f=hstop %9.4f=horbit\n",
232 hour, aostart/(RAD*15.0), aostop/(RAD*15.0), aorbit/(RAD*15.0) );
233 printf( "X%9.4f=aorbit%9.4f=amax%11.4f=amin",
234 aorbit/RAD, amax/RAD, amin/RAD );
235 printf( "%11.3f=azero%9.3f=ahour\n",
236 azero/RAD, ahour/RAD );
237 printf( "X%10.1f=xorbit%10.1f=yorbit", xorbit/KM, yorbit/KM );
238 printf( "%9.1f=distx%10.1f=disty%9.1f=dist\n",
239 distx/KM, disty/KM, dist/KM );
240 printf( "X%9.4f=adist %9.4f=cadmax%9.4f=cadmin%9.4f=cadzero\n",
241 adist/RAD, cadmax, cadmin, cadzero );
242 printf( "X%12.3e= irand %12.3e= imax%13.3e= imin%13.3e= izero\n",
243 irand, imax, imin, izero );
244 printf( "X%11.5f=cos(recv)%13.3e=dist2%13.3e=cp\n",
245 cp*dist2, dist2, cp );
246 printf( "X%12.3e=nlrand%13.3e=nlmax%13.3e=nlmin%13.3e=nlzero\n\n",
247 nlrand, nlmax, nlmin, nlzero );
248 #endif // DEBUG
249
250 } // end of orbit loop
251
252 nlmax *= femax ;
253 nlmin *= femin ;
254 nlzero *= fezero ;
255 nlrand *= ferand ;
256
257 double nlmoon = -cos( ahour );
258
259 fprintf( datfile, "%8.4f %11.4e %11.4e %11.4e %11.4e %11.4e\n",
260 hour, nlrand, nlmax, nlmin, nlzero, nlmoon );
261 }
262 fclose( datfile );
263 }
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.