// orbit04w.c
// really fake orbit
// This approximates a Kepler orbit with high eccentricity
// 
// Given the way I set up the palette, I would expect a white
// background, but it makes a dark one.  Huh?
//
// compile with:  gcc -o orbit04w orbit04w.c -lgd -lpng -lm
//
// Uses the libgd library.  For documentation see the file:
//   /usr/share/doc/gd-*/index.html  
//
// or the website:
//   http://www.libgd.org/Reference

#include "gd.h"
#include "math.h"
#include <stdio.h>

#define  PNGFMT      "orbit04w.png"

#define  XSIZE       1000
#define  SEMIMAJOR   300
#define  YSIZE       700
#define  SATSIZE     100
#define  NFRAME      120
#define  SMEAR       4
#define  LWID        7
#define  FUDGE       5.0
#define  ERROR       1e-4
#define  TAIL        1.5
#define  ECCENTRIC   0.20
// #define  KLUDGE1     (1.0/(1.0-ECCENTRIC*ECCENTRIC))
#define  KLUDGE1     (1.0+2.0/NFRAME)

int main() {
	FILE   *pngout    ;                /* Declare output file */
   char   dirname[80]     ;
   char   framename[80]   ;

   int     sun1, white, red, green, blue;  /* Declare color indexes */
   int     gray, dgray, dwhit, trans ;
   int     dred, cyan, dcyan ;
   double  pi2    = 8.0 * atan( 1.0 );
   double  delta  = pi2 / (NFRAME)  ;
   double  angle  ;
   double  dangle ;
   double  r      ;

   // compute orbit sizes
   
   double  semimajor  = SEMIMAJOR ;
   double  eccentric  = ECCENTRIC ;
   double  xsize      = eccentric*semimajor ;
   double  ysize      = sqrt( semimajor*semimajor - xsize*xsize );
   double  fudge      = (1.0-eccentric) ;
   double  r00        = semimajor * ( 1.0 - eccentric*eccentric ) ; 

   int     centery    = YSIZE/2   ;
   int     centerx    = XSIZE/2   ;
   int     dcell      = SATSIZE/2 ;
   int     derth      = 0.25* XSIZE;

   gdImagePtr im             ;                /* Declare the image */
   gdPoint rpoints[2*NFRAME] ;
   gdPoint cpoints[2*NFRAME] ;
   int     point ;
   int     frame ;
   int     ecount=0   ;
   double  object=0.0 ;
   double  error, error0 ;
   int     ox, oy ;  
   int     ii, jj ;
   int     rpoint, cpoint ;

   //------------------------------------------------------------
 
   // loop until end angle closes
   for( error = 100 ; fabs(error) > ERROR  ; ) {
      point = 0 ;

      // test computation
      for( angle=0.0 ; angle <= pi2 ; ) {
         r  = r00/( 1.0 + eccentric * cos( angle ));
         dangle = delta * semimajor * semimajor / ( r * r ) ;
         angle += dangle ;
         point++ ;
      }
      ecount++ ;
      error  = (double)(point-(NFRAME)) - ((angle-pi2)/dangle) ;
      delta *= ( 1.0+fudge*error/(NFRAME) ) ;
   }

   // compute points for ellipse with line segments
   point = 0 ;
   double cangle = 0.0 ;
   double sm     = semimajor ;
   for( angle=0.0 ; angle <= pi2 ; ) {
      r  = r00/( 1.0 + eccentric * cos( angle ));
      dangle = delta * semimajor * semimajor / ( r * r ) ;
      angle  += dangle ;
      cangle += delta*KLUDGE1  ;
      cpoints[point].x = centerx - (int) ( r   * cos(  angle ) ) ;
      cpoints[point].y = centery + (int) ( r   * sin(  angle ) ) ;
      rpoints[point].x = centerx - (int) ( sm  * cos( cangle ) ) ;
      rpoints[point].y = centery + (int) ( sm  * sin( cangle ) ) ;
      point++ ;
   }

   printf( "%5d points computed, %14.6e residual angle error, now render\n",
            point, angle-pi2 );
   

   im = gdImageCreateTrueColor(XSIZE, YSIZE );

   /* Allocate the colors sun1, white, red, green, blue */
   white = gdImageColorAllocate(im, 255, 255, 255);
   sun1  = gdImageColorAllocate(im,  51,  51, 102);
   red   = gdImageColorAllocate(im, 255,   0,   0);
   dred  = gdImageColorAllocate(im,  96,   0,   0);
   green = gdImageColorAllocate(im,   0, 255,   0);
   cyan  = gdImageColorAllocate(im,   0, 255, 255);
   dcyan = gdImageColorAllocate(im,   0, 160, 160);
   gray  = gdImageColorAllocate(im, 128, 128, 128);
   dwhit = gdImageColorAllocate(im, 192, 192, 192);
   dgray = gdImageColorAllocate(im,  48,  48,  48);
   trans = gdImageColorAllocate(im, 1, 1, 1);

   /* line thickness */
   gdImageSetThickness( im, LWID );

   /* draw orbit */
   gdImagePolygon( im, rpoints, point, dcyan );
   gdImagePolygon( im, cpoints, point, gray  );

   /* earth */
   ox = centerx ;
   oy = centery ;
   gdImageFilledArc(im,  ox    , oy   , derth, derth, 0, 360, green, gdArc );


   for( frame = 0 ; frame < NFRAME ; frame += NFRAME/4 ) {
      // server_sat smear
      for( ii=-SMEAR ; ii < 0 ; ii++ ) {
         jj=frame+ii ;
         if( jj < 0 ) jj += NFRAME ;
         ox = rpoints[jj].x ;
         oy = rpoints[jj].y ;
         gdImageFilledArc(im, ox   , oy   , dcell, dcell, 0, 360, dcyan,  gdArc);

	 ox = cpoints[jj].x ;
         oy = cpoints[jj].y ;
         gdImageFilledArc(im, ox   , oy   , dcell, dcell, 0, 360, gray, gdArc);
      }

      // server_sat main
      ox = rpoints[frame].x ;
      oy = rpoints[frame].y ;
      gdImageFilledArc(im, ox   , oy   , dcell, dcell, 0, 360, cyan , gdArc);

      ox = cpoints[frame].x ;
      oy = cpoints[frame].y ;
      gdImageFilledArc(im, ox   , oy   , dcell, dcell, 0, 360, white , gdArc);
   }
   pngout = fopen( PNGFMT, "wb");
   gdImagePngEx( im, pngout, 1 );
   gdImageDestroy(im);
   fclose(pngout);
}