attachment:orbit04w.c of Brown2013

   1 // orbit04w.c
   2 // really fake orbit
   3 // This approximates a Kepler orbit with high eccentricity
   4 // 
   5 // Given the way I set up the palette, I would expect a white
   6 // background, but it makes a dark one.  Huh?
   7 //
   8 // compile with:  gcc -o orbit04w orbit04w.c -lgd -lpng -lm
   9 //
  10 // Uses the libgd library.  For documentation see the file:
  11 //   /usr/share/doc/gd-*/index.html  
  12 //
  13 // or the website:
  14 //   http://www.libgd.org/Reference
  15 
  16 #include "gd.h"
  17 #include "math.h"
  18 #include <stdio.h>
  19 
  20 #define  PNGFMT      "orbit04w.png"
  21 
  22 #define  XSIZE       1000
  23 #define  SEMIMAJOR   300
  24 #define  YSIZE       700
  25 #define  SATSIZE     100
  26 #define  NFRAME      120
  27 #define  SMEAR       4
  28 #define  LWID        7
  29 #define  FUDGE       5.0
  30 #define  ERROR       1e-4
  31 #define  TAIL        1.5
  32 #define  ECCENTRIC   0.20
  33 // #define  KLUDGE1     (1.0/(1.0-ECCENTRIC*ECCENTRIC))
  34 #define  KLUDGE1     (1.0+2.0/NFRAME)
  35 
  36 int main() {
  37 	FILE   *pngout    ;                /* Declare output file */
  38    char   dirname[80]     ;
  39    char   framename[80]   ;
  40 
  41    int     sun1, white, red, green, blue;  /* Declare color indexes */
  42    int     gray, dgray, dwhit, trans ;
  43    int     dred, cyan, dcyan ;
  44    double  pi2    = 8.0 * atan( 1.0 );
  45    double  delta  = pi2 / (NFRAME)  ;
  46    double  angle  ;
  47    double  dangle ;
  48    double  r      ;
  49 
  50    // compute orbit sizes
  51    
  52    double  semimajor  = SEMIMAJOR ;
  53    double  eccentric  = ECCENTRIC ;
  54    double  xsize      = eccentric*semimajor ;
  55    double  ysize      = sqrt( semimajor*semimajor - xsize*xsize );
  56    double  fudge      = (1.0-eccentric) ;
  57    double  r00        = semimajor * ( 1.0 - eccentric*eccentric ) ; 
  58 
  59    int     centery    = YSIZE/2   ;
  60    int     centerx    = XSIZE/2   ;
  61    int     dcell      = SATSIZE/2 ;
  62    int     derth      = 0.25* XSIZE;
  63 
  64    gdImagePtr im             ;                /* Declare the image */
  65    gdPoint rpoints[2*NFRAME] ;
  66    gdPoint cpoints[2*NFRAME] ;
  67    int     point ;
  68    int     frame ;
  69    int     ecount=0   ;
  70    double  object=0.0 ;
  71    double  error, error0 ;
  72    int     ox, oy ;  
  73    int     ii, jj ;
  74    int     rpoint, cpoint ;
  75 
  76    //------------------------------------------------------------
  77  
  78    // loop until end angle closes
  79    for( error = 100 ; fabs(error) > ERROR  ; ) {
  80       point = 0 ;
  81 
  82       // test computation
  83       for( angle=0.0 ; angle <= pi2 ; ) {
  84          r  = r00/( 1.0 + eccentric * cos( angle ));
  85          dangle = delta * semimajor * semimajor / ( r * r ) ;
  86          angle += dangle ;
  87          point++ ;
  88       }
  89       ecount++ ;
  90       error  = (double)(point-(NFRAME)) - ((angle-pi2)/dangle) ;
  91       delta *= ( 1.0+fudge*error/(NFRAME) ) ;
  92    }
  93 
  94    // compute points for ellipse with line segments
  95    point = 0 ;
  96    double cangle = 0.0 ;
  97    double sm     = semimajor ;
  98    for( angle=0.0 ; angle <= pi2 ; ) {
  99       r  = r00/( 1.0 + eccentric * cos( angle ));
 100       dangle = delta * semimajor * semimajor / ( r * r ) ;
 101       angle  += dangle ;
 102       cangle += delta*KLUDGE1  ;
 103       cpoints[point].x = centerx - (int) ( r   * cos(  angle ) ) ;
 104       cpoints[point].y = centery + (int) ( r   * sin(  angle ) ) ;
 105       rpoints[point].x = centerx - (int) ( sm  * cos( cangle ) ) ;
 106       rpoints[point].y = centery + (int) ( sm  * sin( cangle ) ) ;
 107       point++ ;
 108    }
 109 
 110    printf( "%5d points computed, %14.6e residual angle error, now render\n",
 111             point, angle-pi2 );
 112    
 113 
 114    im = gdImageCreateTrueColor(XSIZE, YSIZE );
 115 
 116    /* Allocate the colors sun1, white, red, green, blue */
 117    white = gdImageColorAllocate(im, 255, 255, 255);
 118    sun1  = gdImageColorAllocate(im,  51,  51, 102);
 119    red   = gdImageColorAllocate(im, 255,   0,   0);
 120    dred  = gdImageColorAllocate(im,  96,   0,   0);
 121    green = gdImageColorAllocate(im,   0, 255,   0);
 122    cyan  = gdImageColorAllocate(im,   0, 255, 255);
 123    dcyan = gdImageColorAllocate(im,   0, 160, 160);
 124    gray  = gdImageColorAllocate(im, 128, 128, 128);
 125    dwhit = gdImageColorAllocate(im, 192, 192, 192);
 126    dgray = gdImageColorAllocate(im,  48,  48,  48);
 127    trans = gdImageColorAllocate(im, 1, 1, 1);
 128 
 129    /* line thickness */
 130    gdImageSetThickness( im, LWID );
 131 
 132    /* draw orbit */
 133    gdImagePolygon( im, rpoints, point, dcyan );
 134    gdImagePolygon( im, cpoints, point, gray  );
 135 
 136    /* earth */
 137    ox = centerx ;
 138    oy = centery ;
 139    gdImageFilledArc(im,  ox    , oy   , derth, derth, 0, 360, green, gdArc );
 140 
 141 
 142    for( frame = 0 ; frame < NFRAME ; frame += NFRAME/4 ) {
 143       // server_sat smear
 144       for( ii=-SMEAR ; ii < 0 ; ii++ ) {
 145          jj=frame+ii ;
 146          if( jj < 0 ) jj += NFRAME ;
 147          ox = rpoints[jj].x ;
 148          oy = rpoints[jj].y ;
 149          gdImageFilledArc(im, ox   , oy   , dcell, dcell, 0, 360, dcyan,  gdArc);
 150 
 151 	 ox = cpoints[jj].x ;
 152          oy = cpoints[jj].y ;
 153          gdImageFilledArc(im, ox   , oy   , dcell, dcell, 0, 360, gray, gdArc);
 154       }
 155 
 156       // server_sat main
 157       ox = rpoints[frame].x ;
 158       oy = rpoints[frame].y ;
 159       gdImageFilledArc(im, ox   , oy   , dcell, dcell, 0, 360, cyan , gdArc);
 160 
 161       ox = cpoints[frame].x ;
 162       oy = cpoints[frame].y ;
 163       gdImageFilledArc(im, ox   , oy   , dcell, dcell, 0, 360, white , gdArc);
 164    }
 165    pngout = fopen( PNGFMT, "wb");
 166    gdImagePngEx( im, pngout, 1 );
 167    gdImageDestroy(im);
 168    fclose(pngout);
 169 }