attachment:ecl03.c of Eclipse

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