// ns03.c
// v0.1 KHL 2009-11-02
// gcc -o ns03 ns03.c -lm ; ./ns03

#define OUTFILE   "ns03.dat"
#define IPOINTS   10000               // integration points
#define LMIN     -90.0                //
#define LMAX      90.0                //
#define LDEL      1.0                 //

#define DEBUG     1
#undef  DEBUG

#define RE        6378000.0           // earth radius
#define RS        12789000.0          // satellite ring distance
#define B         0.075               // server sat radius
#define IS        50000.0             // sun lux times albedo with attenuation
#define NSAT      80E9                // number of satellites
#define PI        3.1415926535897932

#include        <stdlib.h>
#include        <stdio.h>
#include        <math.h>
#include        <string.h>
#include        <time.h>

FILE            *outdat ;

int main () {
   double     xs, ys  ;               // satellite position
   double     xe, ye  ;               // earth position

   double     res2    ;               // satellite to earth distance, squared
   double     xes     ;               // "" x portion
   double     yes     ;               // "" y portion
   double     gamma   ;               // satellite to earth angle
   double     lng     ;               // rangle earth longitude
   double     hr      ;               // hour of longitude angle
   double     il      ;               // earth longitude degrees
   int        is      ;               // satellite integration count
   int        iside   ;               // left or right integration
   double     incr    ;               // intensity increment
   double     isum    ;               // sum of intensity
// two illumination bands
   // horizon
   double     thll    ;               // rangle, left of left illumination
   double     thrr    ;               // rangle, right of right illumination
   double     thw =  acos(RE/RS);     //
   // earth shadow
   double     thrl = asin(RE/RS);     // rangle, right of right illumination
   double     thlr = -thrl      ;     // rangle, right of left illumination

   double     thtot   ;               // total angle
   double     theta   ;               // satellite angle
   double     thdelta ;               // satellite angle increment
   double     d  = 180.0/PI ;         // radians to degrees
   double     k  =0.001     ;         // meters to kilometers

   // open data output file
   outdat = fopen( OUTFILE, "w" ) ;

   for( il=LMIN ; il <= LMAX+1e-6 ; il += LDEL ) {   // longitude loop
      lng  = il/d ;
      thll = lng-thw ;
      thrr = lng+thw ;
 
      thtot = 0.0    ;
      if( thlr > thll ) { thtot += thlr-thll ; }
      if( thrr > thrl ) { thtot += thrr-thrl ; }
      thdelta = thtot / IPOINTS ;
      iside   = (int) ( 0.49 + ( thlr-thll ) / thdelta );
      if (iside < 0 ) { iside = 0 ; }

      xe = RE * sin( lng );
      ye = RE * cos( lng );

      isum = 0 ;

#ifdef DEBUG
         fprintf( outdat,
                  " thtot=%8.2f thdelta=%8.2f iside=%4d\n",
                    thtot*d,    thdelta*d,    iside );
         fprintf( outdat,
                  " thll=%8.2f  thlr=%8.2f  thrl=%8.2f  thrr=%8.2f\n",
                    thll*d,     thlr*d,     thrl*d,     thrr*d );
                         
         fprintf( outdat,
                  " longitude=%9.2f xe=%9.2f ye=%9.2f\n",
                   lng*d,   xe*k,      ye*k );
#endif

      for( is=0 ; is < IPOINTS ; is ++ ) {  // satellite loop

         // find angle
         if( is < iside ) {  // left illumination
            theta = (0.5+(double)is)*thdelta+thll  ;
         }
         else {              // right illumination
            theta = (0.5+(double)(is-iside))*thdelta+thrl ;
         }
  
         xs = RS * sin( theta );
         ys = RS * cos( theta );

         xes = xe-xs ;
         yes = ye-ys ;
         res2 = xes*xes+yes*yes;
	 gamma = atan2( yes, xes );
	 incr  = fabs( sin( 0.5*gamma+0.25*PI ) * sin( gamma-lng ) / res2 ) ;
         isum += incr ;

#ifdef DEBUG
         fprintf( outdat,
                  "th=%8.2f xs=%9.2f ys=%9.2f r=%9.2f g=%8.2f inc=%9.2e\n",
                   theta*d,  xs*k, ys*k, sqrt(res2)*k,  gamma*d,  incr );
#endif

      } // end of satellite integration

      // scaling 
      //  ( thdelta over 2 pi ) angular element
      //  2 over pi^2 scaling times B^2

      isum *= NSAT*IS*B*B*thdelta/(PI*PI*PI) ;

      hr = il/15.0 ;

      fprintf( outdat,  "%7.2f %12.4e\n", hr, isum ); // one data point

   } // end of longitude loop

   fclose( outdat );
   return 0;
}