attachment:ns03.c of BiologicalEffects

   1 // ns03.c
   2 // v0.1 KHL 2009-11-02
   3 // gcc -o ns03 ns03.c -lm ; ./ns03
   4 
   5 #define OUTFILE   "ns03.dat"
   6 #define IPOINTS   10000               // integration points
   7 #define LMIN     -90.0                //
   8 #define LMAX      90.0                //
   9 #define LDEL      1.0                 //
  10 
  11 #define DEBUG     1
  12 #undef  DEBUG
  13 
  14 #define RE        6378000.0           // earth radius
  15 #define RS        12789000.0          // satellite ring distance
  16 #define B         0.075               // server sat radius
  17 #define IS        50000.0             // sun lux times albedo with attenuation
  18 #define NSAT      80E9                // number of satellites
  19 #define PI        3.1415926535897932
  20 
  21 #include        <stdlib.h>
  22 #include        <stdio.h>
  23 #include        <math.h>
  24 #include        <string.h>
  25 #include        <time.h>
  26 
  27 FILE            *outdat ;
  28 
  29 int main () {
  30    double     xs, ys  ;               // satellite position
  31    double     xe, ye  ;               // earth position
  32 
  33    double     res2    ;               // satellite to earth distance, squared
  34    double     xes     ;               // "" x portion
  35    double     yes     ;               // "" y portion
  36    double     gamma   ;               // satellite to earth angle
  37    double     lng     ;               // rangle earth longitude
  38    double     hr      ;               // hour of longitude angle
  39    double     il      ;               // earth longitude degrees
  40    int        is      ;               // satellite integration count
  41    int        iside   ;               // left or right integration
  42    double     incr    ;               // intensity increment
  43    double     isum    ;               // sum of intensity
  44 // two illumination bands
  45    // horizon
  46    double     thll    ;               // rangle, left of left illumination
  47    double     thrr    ;               // rangle, right of right illumination
  48    double     thw =  acos(RE/RS);     //
  49    // earth shadow
  50    double     thrl = asin(RE/RS);     // rangle, right of right illumination
  51    double     thlr = -thrl      ;     // rangle, right of left illumination
  52 
  53    double     thtot   ;               // total angle
  54    double     theta   ;               // satellite angle
  55    double     thdelta ;               // satellite angle increment
  56    double     d  = 180.0/PI ;         // radians to degrees
  57    double     k  =0.001     ;         // meters to kilometers
  58 
  59    // open data output file
  60    outdat = fopen( OUTFILE, "w" ) ;
  61 
  62    for( il=LMIN ; il <= LMAX+1e-6 ; il += LDEL ) {   // longitude loop
  63       lng  = il/d ;
  64       thll = lng-thw ;
  65       thrr = lng+thw ;
  66  
  67       thtot = 0.0    ;
  68       if( thlr > thll ) { thtot += thlr-thll ; }
  69       if( thrr > thrl ) { thtot += thrr-thrl ; }
  70       thdelta = thtot / IPOINTS ;
  71       iside   = (int) ( 0.49 + ( thlr-thll ) / thdelta );
  72       if (iside < 0 ) { iside = 0 ; }
  73 
  74       xe = RE * sin( lng );
  75       ye = RE * cos( lng );
  76 
  77       isum = 0 ;
  78 
  79 #ifdef DEBUG
  80          fprintf( outdat,
  81                   " thtot=%8.2f thdelta=%8.2f iside=%4d\n",
  82                     thtot*d,    thdelta*d,    iside );
  83          fprintf( outdat,
  84                   " thll=%8.2f  thlr=%8.2f  thrl=%8.2f  thrr=%8.2f\n",
  85                     thll*d,     thlr*d,     thrl*d,     thrr*d );
  86                          
  87          fprintf( outdat,
  88                   " longitude=%9.2f xe=%9.2f ye=%9.2f\n",
  89                    lng*d,   xe*k,      ye*k );
  90 #endif
  91 
  92       for( is=0 ; is < IPOINTS ; is ++ ) {  // satellite loop
  93 
  94          // find angle
  95          if( is < iside ) {  // left illumination
  96             theta = (0.5+(double)is)*thdelta+thll  ;
  97          }
  98          else {              // right illumination
  99             theta = (0.5+(double)(is-iside))*thdelta+thrl ;
 100          }
 101   
 102          xs = RS * sin( theta );
 103          ys = RS * cos( theta );
 104 
 105          xes = xe-xs ;
 106          yes = ye-ys ;
 107          res2 = xes*xes+yes*yes;
 108 	 gamma = atan2( yes, xes );
 109 	 incr  = fabs( sin( 0.5*gamma+0.25*PI ) * sin( gamma-lng ) / res2 ) ;
 110          isum += incr ;
 111 
 112 #ifdef DEBUG
 113          fprintf( outdat,
 114                   "th=%8.2f xs=%9.2f ys=%9.2f r=%9.2f g=%8.2f inc=%9.2e\n",
 115                    theta*d,  xs*k, ys*k, sqrt(res2)*k,  gamma*d,  incr );
 116 #endif
 117 
 118       } // end of satellite integration
 119 
 120       // scaling 
 121       //  ( thdelta over 2 pi ) angular element
 122       //  2 over pi^2 scaling times B^2
 123 
 124       isum *= NSAT*IS*B*B*thdelta/(PI*PI*PI) ;
 125 
 126       hr = il/15.0 ;
 127 
 128       fprintf( outdat,  "%7.2f %12.4e\n", hr, isum ); // one data point
 129 
 130    } // end of longitude loop
 131 
 132    fclose( outdat );
 133    return 0;
 134 }