attachment:ll02.c of Radiation

   1 // ll02.c
   2 // compute McIlwain L layer versus angle in server sky orbit
   3 // also compute the histogram of probability versus L shell
   4 //
   5 // gcc -o ll02 ll02.c -lm 
   6 
   7 #define  NAME      "ll02"
   8 
   9 // #define  NPOINTS   37        //
  10 // #define  NBINS     10        //
  11 // #define  NSAMPLE   20        //
  12 
  13 #define  NPOINTS   361        //
  14 #define  NBINS     100        //
  15 #define  NSAMPLE   100000     //
  16 
  17 #define  DEBUG     1          //
  18 #define  AMIN      0.0        //  orbit angle
  19 #define  AMAX      360.0      //
  20 
  21 #define  RE        6371.0     //  McIlwain earth radius kilometers
  22 #define  RORB      12789.0    //  Server Sky orbit radius kilometers 
  23 
  24 #define  ANG_N     10.26      //  field tilt from north
  25 #define  ANG_E    -71.78      //  field tilt to east
  26 
  27 #define  OFF_R     552.0      //  field offset in kilometers
  28 #define  OFF_N     22.2       //  field offset tilt from north
  29 #define  OFF_E     141.6      //  field offset tilt towards east
  30 
  31 #define  RAD       (3.1415926536/180.0) 
  32 
  33 #include <math.h>
  34 #include <string.h>
  35 #include <stdio.h>
  36 #include <stdlib.h>
  37 
  38 //-----------------------------------------------------------------
  39 double  h_o ;
  40 double  x_o ;
  41 double  y_o ;
  42 double  z_o ;
  43 
  44 double  lm, bm ;
  45 
  46 double  lmin, lmax ;
  47 
  48 //-----------------------------------------------------------------
  49 void lshell( double ang ) {
  50 
  51    // position in earth radius=1 coordinates 
  52 
  53    double rorb = RORB/RE ;
  54    double xorb = rorb * sin( RAD * ang );
  55    double yorb = 0.0             ;
  56    double zorb = rorb * cos( RAD * ang );
  57 
  58    // position relative to dipole offset
  59 
  60    double x1   = xorb - x_o ;
  61    double y1   = yorb - y_o ;
  62    double z1   = zorb - z_o ;
  63 
  64    // rotate relative to dipole east-west
  65 
  66    double a2   = RAD * ANG_E ;
  67    double x2   = x1 * cos(a2) + z1 * sin(a2) ;
  68    double y2   = y1                          ;
  69    double z2   = z1 * cos(a2) - x1 * sin(a2) ;
  70 
  71    // rotate relative to dipole north-south 
  72 
  73    double am   = RAD * ANG_N ;
  74    double xm   = x2 * cos(am) + y2 * sin(am) ;
  75    double ym   = y2 * cos(am) - x2 * sin(am) ;
  76    double zm   = z2                          ;
  77    double xzm  = sqrt( xm*xm + zm*zm )       ;
  78    double rm   = sqrt( xzm*xzm + ym*ym )     ;
  79   
  80    // find McIlwain parameters               
  81 
  82    double phim = atan2( ym , xzm )           ;
  83    double sinm = xzm / rm                    ;
  84    double cosm = ym  / rm                    ;
  85           lm   = rm / ( sinm * sinm )        ;
  86           bm   = sqrt( 3.0*cosm*cosm+1.0 )
  87                / ( rm*rm*rm )                ;
  88 }
  89 
  90 //-----------------------------------------------------------------
  91 int main() {
  92 
  93    FILE    *datfile                     ; //
  94    char    gnuplot[200]                 ; //
  95    char    filename[80]                 ; //
  96    int     iang                         ; //
  97    
  98    lmin =  1000.0 ;
  99    lmax = -1000.0 ;
 100 
 101    h_o = OFF_R * cos( RAD * OFF_N ) ;
 102    x_o =   h_o * sin( RAD * OFF_E ) / RE ;
 103    y_o = OFF_R * sin( RAD * OFF_N ) / RE ;
 104    z_o =  -h_o * cos( RAD * OFF_E ) / RE ;
 105 
 106    // --------------------------------------------------------------------
 107 
 108    // compute angle loop
 109 
 110    datfile  = fopen( "tmp.dat", "wb" );
 111    for( iang=0 ; iang < NPOINTS ; iang++ ) {
 112       double ang  = AMIN + (AMAX-AMIN)*((double)iang )/((double)(NPOINTS-1));
 113 
 114       lshell( ang ) ;
 115       if( lmin > lm ) { lmin = lm ; }
 116       if( lmax < lm ) { lmax = lm ; }
 117 
 118       fprintf( datfile, "%16.4f%16.4f%16.4f\n", ang, lm, bm ) ; 
 119    }
 120    fclose( datfile );
 121 
 122    lmax += 1e-4 ;
 123    lmin -= 1e-4 ;
 124 
 125    // histogram
 126 
 127    double  prob[NBINS], lval     ; //
 128    double  delta = ((double) NBINS )/((double) NSAMPLE );
 129 
 130    // histogram initialization 
 131    for( iang=0 ; iang < NBINS ; iang++ ) { 
 132       prob[ iang ] = 0.0 ;
 133    }
 134 
 135    double lscale = ((double)NBINS) /(lmax-lmin);
 136 
 137    // histogram compute
 138    for( iang=0 ; iang < NSAMPLE ; iang++ ) {
 139       double ang  = AMIN + (AMAX-AMIN)*((double)iang )/((double)(NSAMPLE));
 140       lshell( ang ) ;
 141       int bin = (int) ( (lm-lmin)*lscale ) ;
 142       prob[ bin ] += delta ;
 143    }
 144 
 145    // histogram printout
 146    datfile  = fopen( "tmp1.dat", "wb" );
 147    for( iang=0 ; iang < NBINS ; iang++ ) {
 148       lval = lmin + (lmax-lmin)*(0.5+(double)iang)/((double) NBINS );
 149       fprintf( datfile, "%16.4f%16.4f%16.4f\n", prob[ iang ], lval ) ; 
 150    }
 151    fclose( datfile );
 152 
 153    sprintf( gnuplot, "/usr/local/bin/gp %s.gp", NAME );
 154    system(  gnuplot );
 155    return(0);
 156 }