## Attachment 'nl02.c'

```   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)
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
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 ;
195             azero  = atan2( yorbit, xorbit-RE );
196             irand  = 0.25 ;
197
198          } else if( xorbit < -RE ) { // eastward of earth
199             amax   = 0.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
213
214 	 // Lambert, transmit angle
218
220          double cp    = cos(ahour-adist)/dist2 ;
221
222          if( cp > 0.0 ) {
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",
233          printf( "X%9.4f=aorbit%9.4f=amax%11.4f=amin",
235          printf( "%11.3f=azero%9.3f=ahour\n",
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 );
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.
• [get | view] (2013-08-06 04:51:46, 9.0 KB) [[attachment:g400.c]]
• [get | view] (2013-08-06 04:52:00, 4004.5 KB) [[attachment:g400.swf]]
• [get | view] (2012-06-04 22:29:30, 8.9 KB) [[attachment:nl02.c]]
• [get | view] (2012-06-04 22:14:08, 11.2 KB) [[attachment:nl02nlp.png]]
• [get | view] (2012-06-04 22:15:23, 10.4 KB) [[attachment:nl02pow.png]]
All files | Selected Files: delete move to page