// eclipse period versus time of year
// cc ecl01.c -lm -o ecl01
//
//  This plots the eclipse time in second versus time of year 
//  The time of year is expressed in degrees, not days or months,
//  which is a little odd but much easier to visualize

#define  NAME    "ecl01"
#define  OFILE   "ecl01.txt"

#define  YSTEP   1.0         // degrees
#define  YEAR    360.000     //
#define  M288    12788866.0  // m288 orbit radius meters
#define  R_M     (M288)      // orbit radius
#define  TSUN    14400.0     // orbit time

#define  R_P     6356752.0   // polar earth radius meters
#define  R_E     6378137.0   // equatorial earth radius  meters
#define  TILT    23.439281   // degrees

#define  FNT    "DejaVuMonoSansBold"
#define  FSZ    16

// MAJOR KLUDGE WARNING!  These plot corners depend on what Gnuplot draws:
#define  XGP_UL    168
#define  YGP_UL     83

#include <math.h>
#include <gd.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>

FILE     *file          ; //
char     gnuplot[200]   ; //

// --------------------------------------------------------------------------

int main() {
   double  pi     = 4.0*atan(1.0);
   double  d2r    = pi/180.0     ;
   double  yang                  ;
   double  phi = d2r * TILT ;
   double  kr0 = sin(phi)*sin(phi) ;
   double  kr1 = -kr0*( 1.0 - ( (R_M*R_M + R_P*R_P) / (R_E*R_E) ) ) ;
   
   double  xea   = 0.0 ;
   double  fra   = 0.0 ;
   double  tea   = 0.0 ;

   printf( "%12.9f=kr1  %12.9f=kr0\n", kr1, kr0 );

   file = fopen( OFILE, "w" );

   for( yang=0.0 ; yang < YEAR-0.001 ; yang++ ) {
      double beta = d2r * yang ;
      double sin2beta = sin(beta)*sin(beta) ;
      double xe   = R_E *sqrt((1.0-kr1*sin2beta)/(1.0-kr0*sin2beta)) ;
      double frac = asin( xe/R_M )/pi;
      double te   = TSUN * frac ;
      fprintf( file,"%9.4f%9.6f%12.2f%12.4f\n",
                yang, frac, xe, te ); // add a plot point
      xea    += xe   ;
      fra    += frac ;
      tea    += te   ;
   }
   fclose( file );

   // plot angle versus time

   sprintf( gnuplot, "/usr/local/bin/gp %s.gp", NAME );
   system( gnuplot );

   printf( "avg %11.6ffr%14.2fx%14.4ft\n",
                 fra/YEAR, xea/YEAR, tea/YEAR ); 
   return(0);
}