attachment:ecl01.c of Eclipse

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