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.
  • [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

You are not allowed to attach a file to this page.