attachment:ecl02.c of Eclipse

   1 // eclipse period versus time of year
   2 // cc ecl02.c -lm -o ecl02
   3 //
   4 //  This plots the eclipse time in second versus time of year 
   5 //  The time of year is expressed in degrees, not days or months,
   6 //  which is a little odd but much easier to visualize
   7 
   8 #define  NAME    "ecl02"
   9 #define  OFILE   "ecl02.txt"
  10 
  11 #define  YSTEP   1.0         // degrees
  12 #define  YEAR    360.000     //
  13 #define  M288    12788866.0  // m288 orbit radius meters
  14 #define  MGEO    42614000.0  // geo orbit radius meters
  15 #define  T288    14400.0     // m288 orbit time
  16 #define  TGEO    86400.0     // geo orbit time
  17 #define  R_M     (MGEO)      // orbit radius
  18 #define  TSUN    (TGEO)      // sun-relative orbit time
  19 
  20 #define  R_P     6356752.0   // polar earth radius meters
  21 #define  R_E     6378137.0   // equatorial earth radius  meters
  22 #define  TILT    23.439281   // degrees
  23 
  24 #define  FNT    "DejaVuMonoSansBold"
  25 #define  FSZ    16
  26 
  27 // MAJOR KLUDGE WARNING!  These plot corners depend on what Gnuplot draws:
  28 #define  XGP_UL    168
  29 #define  YGP_UL     83
  30 
  31 #include <math.h>
  32 #include <gd.h>
  33 #include <string.h>
  34 #include <stdio.h>
  35 #include <stdlib.h>
  36 
  37 FILE     *file          ; //
  38 char     gnuplot[200]   ; //
  39 
  40 // --------------------------------------------------------------------------
  41 
  42 int main() {
  43    double  pi     = 4.0*atan(1.0);
  44    double  d2r    = pi/180.0     ;
  45    double  yang                  ;
  46    double  phi = d2r * TILT ;
  47    double  kr0 = sin(phi)*sin(phi) ;
  48    double  kr1 = -kr0*( 1.0 - ( (R_M*R_M + R_P*R_P) / (R_E*R_E) ) ) ;
  49    
  50    double  xea   = 0.0 ;
  51    double  fra   = 0.0 ;
  52    double  tea   = 0.0 ;
  53 
  54    printf( "%12.9f=kr1  %12.9f=kr0\n", kr1, kr0 );
  55 
  56    file = fopen( OFILE, "w" );
  57 
  58    for( yang=0.0 ; yang < YEAR-0.001 ; yang++ ) {
  59       double beta = d2r * yang ;
  60       double sin2beta = sin(beta)*sin(beta) ;
  61       double xe    ;
  62       double numer = (1.0-kr1*sin2beta) ;
  63       if( numer > 0.0 ) { xe = R_E *sqrt(numer/(1.0-kr0*sin2beta)) ; }
  64       else              { xe = 0.0 ;                                 }
  65     
  66       double frac = asin( xe/R_M )/pi;
  67       double te   = TSUN * frac ;
  68       fprintf( file,"%9.4f%9.6f%12.2f%12.4f\n",
  69                 yang, frac, xe, te ); // add a plot point
  70       xea    += xe   ;
  71       fra    += frac ;
  72       tea    += te   ;
  73    }
  74    fclose( file );
  75 
  76    // plot angle versus time
  77 
  78    sprintf( gnuplot, "/usr/local/bin/gp %s.gp", NAME );
  79    system( gnuplot );
  80 
  81    printf( "avg %11.6ffr%14.2fx%14.4ft\n",
  82                  fra/YEAR, xea/YEAR, tea/YEAR ); 
  83    return(0);
  84 }