// ll288a.c
// compute McIlwain L layer versus angle in server sky orbit
// also compute the histogram of probability versus L shell
// include eccentricity
//
// gcc -o ll288a ll288a.c -lm 

#define  NAME      "M288a"

// #define  NPOINTS   37        //
// #define  NBINS     10        //
// #define  NSAMPLE   20        //

#define  NPOINTS   361        //
#define  NBINS     229        //
#define  NSAMPLE   360000     //

#define  DEBUG     1          //
#define  AMIN      -180.0     //  orbit longitude
#define  AMAX       180.0     //

#define  RE        6371.0     //  McIlwain earth radius kilometers
#define  M288      12789.0    //  M288
#define  M360      14441.0    //  M360
#define  M480      16756.0    //  M480
#define  ECC       0.0226     //  Eccentricity
#define  RORB      (M288)     //  Orbit semimajor axis kilometers

#define  ANG_N     10.26      //  field tilt from north
#define  ANG_E    -71.78      //  field tilt to east

#define  OFF_R     552.0      //  field offset in kilometers
#define  OFF_N     22.2       //  field offset tilt from north
#define  OFF_E     141.6      //  field offset tilt towards east

#define  RAD       (3.1415926536/180.0) 

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

//-----------------------------------------------------------------
double  h_o ;
double  x_o ;       // equatorial plane
double  y_o ;       // NORTH (as POVRAY)
double  z_o ;       // equatorial plane (as POVRAY)

double  lm, bm, phim ;
double  lmin, lmax ;
double  ecc1 ;

//-----------------------------------------------------------------
void lshell( double ang ) {

   // position in earth radius=1 coordinates 

   double rorb = ecc1*RORB/RE ;
   double xorb = -rorb * sin( RAD * ang );
   double yorb = 0.0             ;
   double zorb = -rorb * cos( RAD * ang );

   // position relative to dipole offset

   double x1   = xorb - x_o ;
   double y1   = yorb - y_o ;
   double z1   = zorb - z_o ;

   // rotate relative to dipole east-west

   double a2   = RAD * ANG_E ;
   double x2   = x1 * cos(a2) + z1 * sin(a2) ;
   double y2   = y1                          ;
   double z2   = z1 * cos(a2) - x1 * sin(a2) ;

   // rotate relative to dipole north-south 

   double am   = RAD * ANG_N ;
   double xm   = x2 * cos(am) + y2 * sin(am) ;
   double ym   = y2 * cos(am) - x2 * sin(am) ;
   double zm   = z2                          ;
   double xzm  = sqrt( xm*xm + zm*zm )       ;
   double rm   = sqrt( xzm*xzm + ym*ym )     ;

   // find McIlwain parameters               

          phim = atan2( ym , xzm )           ;
   double sinm = xzm / rm                    ;
   double cosm = ym  / rm                    ;
          lm   = rm / ( sinm * sinm )        ;
          bm   = sqrt( 3.0*cosm*cosm+1.0 )
               / ( rm*rm*rm )                ;
  
   /*
   printf( "%10.3f_ang%10.3f_xm%10.3f_ym%10.3f_zm%10.3f_phim\n",
            ang,      xm,      ym,      zm,      phim       );
   */
}

//-----------------------------------------------------------------
int main() {

   FILE    *datfile                     ; //
   char    gnuplot[200]                 ; //
   char    filename[80]                 ; //
   int     iang                         ; //
   int     iecc                         ; //
   
   lmin =  1000.0 ;
   lmax = -1000.0 ;

   h_o = OFF_R * cos( RAD * OFF_N ) ;
   x_o =   h_o * sin( RAD * OFF_E ) / RE ;
   y_o = OFF_R * sin( RAD * OFF_N ) / RE ;
   z_o =  -h_o * cos( RAD * OFF_E ) / RE ;

   // --------------------------------------------------------------------
   // compute angle loop

   sprintf( filename, "%ss.dat", NAME );
   datfile  = fopen( filename, "wb" );
   for( iang=0 ; iang < NPOINTS ; iang++ ) {
      double ang  = AMIN + (AMAX-AMIN)*((double)iang )/((double)(NPOINTS-1));

      fprintf( datfile, "%16.4f", ang );

      for( iecc = -1 ; iecc <= 1 ; iecc++ ) {
         ecc1 = 1.0 + ECC*((double)iecc) ;
         lshell( ang ) ;
         if( lmin > lm ) { lmin = lm ; }
         if( lmax < lm ) { lmax = lm ; }

         double lamm = fabs( phim/RAD ) ;
         fprintf( datfile, "%16.4f%16.4f%16.4f", lm, bm, lamm ) ; 
      }
      fprintf( datfile, "\n" );
   }
   fclose( datfile );

   lmax += 1e-4 ;
   lmin -= 1e-4 ;

   // histogram

   fprintf( stderr, "%12.6f%12.6f\n", lmin, lmax ); fflush( stderr );
  
   double  prob[NBINS], lval     ; //

   // histogram initialization 
   for( iang=0 ; iang < NBINS ; iang++ ) { prob[ iang ] = 0.0 ; }

   double lscale = ((double)NBINS) /(lmax-lmin);

   // histogram compute

   double  ecca[360] ;
   double  delta = ((double) NBINS )/(360.0*((double) NSAMPLE ));
   double  ascale = (AMAX-AMIN)/((double)(NSAMPLE)) ;
   int     nb1   = NBINS-1 ;

   for( iecc = 0 ; iecc < 360 ; iecc++ ) {
     ecca[iecc] = 1.0 + ECC * cos( RAD * ( (double) iecc ) );
   }

   for( iang=0 ; iang < NSAMPLE ; iang++ ) {
      double ang  = AMIN + ascale*((double)iang );
      for( iecc = 0 ; iecc < 360 ; iecc++ ) {
         ecc1    = ecca[iecc] ;
         lshell( ang ) ;
         if( ( lm > lmax ) || ( lm < lmin ) ) {
            fprintf( stdout, "%12.6f%12.6f%12.6f\n", ang, ecc1, lm );
            fflush(  stdout );
         }
         int bin = (int) ( (lm-lmin)*lscale ) ;
         if( bin < 0   ) { bin = 0   ; }
         if( bin > nb1 ) { bin = nb1 ; }
         prob[ bin ] += delta ;
      }
   }

   // histogram printout
   sprintf( filename, "%sh.dat", NAME );
   datfile  = fopen( filename, "wb" );
 
   double  dlval = (lmax-lmin)/((double) NBINS );

   lval = lmin + 0.5*dlval ;
   fprintf( datfile, "%16.4f%16.4f%16.4f\n", 0.0, lval ) ; 

   for( iang=0 ; iang < NBINS ; iang++ ) {
      lval = lmin + dlval * (0.5+(double)iang) ;
      fprintf( datfile, "%16.4f%16.4f%16.4f\n", prob[ iang ], lval ) ; 
   }

   lval = lmax - 0.5*dlval ;
   fprintf( datfile, "%16.4f%16.4f%16.4f\n", 0.0, lval ) ; 

   fclose( datfile );

   // prepare gnuplot controp file
   sprintf( gnuplot, "cat XXXXX3.gp | sed -e \"s/XXXXX/%s/\" > %s.gp",
                                                        NAME,  NAME );
   system(  gnuplot );

   // run gnuplot
   sprintf( gnuplot, "/usr/local/bin/gp %s.gp", NAME );
   system(  gnuplot );
   return(0);
}