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

#define  NAME      "M288"

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

#define  NPOINTS   361        //
#define  NBINS     100        //
#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  RORB      (M288)     //  Orbit radius 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 ;

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

   // position in earth radius=1 coordinates 

   double rorb = 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                         ; //
   
   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));

      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%16.4f\n", ang, lm, bm, lamm ) ; 
   }
   fclose( datfile );

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

   // histogram

   double  prob[NBINS], lval     ; //
   double  delta = ((double) NBINS )/((double) NSAMPLE );

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

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

   // histogram compute
   for( iang=0 ; iang < NSAMPLE ; iang++ ) {
      double ang  = AMIN + (AMAX-AMIN)*((double)iang )/((double)(NSAMPLE));
      lshell( ang ) ;
      int bin = (int) ( (lm-lmin)*lscale ) ;
      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 XXXXX.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);
}