// gcc -o gsr02 gsr02.c -lgd -lpng -lm
#define  NAME   "gsr02"
//
// Keith Lofstrom   2013 March 4
//
// You are free to copy and modify this code under the GNU Public
// License version 2,   http://www.gnu.org/licenses/gpl-2.0.html
//
// I wrote this to run under Linux.  It uses the libgd and the libpng
// libraries.  I'm pretty sure you can get this working with BSD, Mac,
// Windows, SunOS, or just about any other OS with a gcc compiler and
// those libraries.
//
// 99% of the time is spent in the tight loop around line 950
//
// V should probably be a command line parameter.  Indeed, the
// many of the defines could be read from a configuration file.
// Also, this really should be modularized - the contruction of
// the geodesic sphere, the movement animations, and the radio
// computations could be split into three small "library" .c files.

// #define  V      57         // tesselation size , 32492 elements
#define  V      16         // tesselation size ,  1692 elements
// #define  V      3         // tesselation size , 92 elements

// NM sets the array sizes - a smarter person would use a malloc()
#define  NM     600        // max tesselation size, 3,600,002 elements

#define  SPACE       180   // scaling of big drawing
#define  XSIZE      1024   // Y screen size (0=left, 1023=right)
#define  YSIZE       768   // Y screen size (0=top,  599=bottom)
#define  YSCENT     60.0   // Y center of small top images
#define  YBCENT    386.0   // Y center of big bottom image
#define  XBCENT    250.0   // Y center of big bottom image
#define  SIZE          6   // size of spot

// note - these are applied to the display, not the stored
// geodesic sphere display, which is always scaled to unity

#define  XEXP        1.2   // expand array in x (orbital) direction
#define  YEXP        1.0   // expand array in y (NS)      direction
#define  ZEXP        0.2   // expand array in z (radial)  direction

// turn the array just a bit, so objects in back are visible past those
// in front
#define  XZROTATE   0.01

#define  ASKEW       2.0   // apogee skew 

#define  NPICS       240   // number of images in animation
#define  RATE         10   // images per second
#define  ESIZE        80   // diameter of earth disk
#define  ORBSIZE      80   // radius of orbit
#define  ZFACTR      0.2   // scaling for items in back
#define  ASCALE       12   // earth orbiting arrow scale
#define  OSCALE    0.125   // earth orbiting scale
#define  CMID      128.0   // midrange color in daylight
#define  CSCALE    125.0   // color scaling to distance
#define  CDIM       0.25   // dimming factor for eclipse

#define  FNT        "DejaVuMonoSans"
#define  FSZ          40

#define  SPEEDUP     (14400.0*RATE/NPICS)

// radio  --------------------------------------------------------

#define  WAVELENGTH  0.005 // radio wavelength 
#define  TANGX       0.0   // target X angle in radians
#define  TANGY       0.0   // target Y angle in radians
#define  DANGX       0.0   // output display center X angle in kilometers
#define  DANGY       0.0   // output display center Y angle in kilometers
#define  DANGH      14.0   // output display angle height in kilometers
#define  DANGW      14.0   // output display angle width  in kilometers
#define  DANGS       2.0   // output display angle step in kilometers
#define  DANGSCALE  7200.0 // kilometers per radian
#define  OX0         540   // output array left
#define  OX1        1020   // output array right ( 480 wide )
#define  OY0         150   // output array top
#define  OY1         630   // output array bottom ( 360 high )

#define  C5MAX       320
#define  C2MAX       128
#define  C1MAX        64
#define  MMTHR       0.4
#define  MMMIN      1E-6

// --------------------------------------------------------------
// How toroidal orbits work.
//
// An elliptical orbit is faster than a similar circular orbit at
// perigee, and slower at apogee.   If the delta between circular
// and elliptical is ZEXP ( == e R ) then an array element halfway
// between perigee and apogee in an elliptical orbit is displaced
// 2 * ZEXP ahead of its circular neighbor.  We need enough ZEXP
// to keep perigee and apogee objects from colliding, and also to
// provide a reasonably large surface facing forwards or backwards
// in orbit at the 6AM and 6PM positions of the orbit.
// 
// We need enough YEXP (North/South) size to make a reasonable 
// large disk for radio and solar capture.   We need enough XEXP
// to make radio disk, and capture the sun at the noon position.
//
// You can fiddle with the relative values of XEXP, YEXP, and ZEXP
// to see how it shapes the array.  We probably want something a
// little flattened in ZEXP to reduce relative skew. 
//
// This program will arbitrarily draw "right hand" toroidal orbits;
// If your thumb is the orbital direction, the orbits will turn
// in the direction of your fingers.  Viewing from just outside the
// orbit, looking at the center of the earth, the orbit will evolve
// like so:

//  angle   Z (radius)             Y(NS)      X orbit
//    0°    -ZEXP inside,perigee   0          0
//   90°    0                      YEXP top   +2*ZEXP skew forward
//  180°    +ZEXP outside,apogee   0          0
//  270°    0                      YEXP bot   -2*ZEXP skew backward

// --------------------------------------------------------------
// How to draw a geodesic sphere - tesselated icosahedron
//
// make a big array of normalized array points
// for a geodesic sphere radius 1 
//   class 1 icosahedral base
//   V up to whole bunches, possibly millions of points
//   number of vertexes is 10*V^2 + 2
//     
//  12 vertices
//  20 faces
//  30 edges
// 
//  Procedure:
//  (1) define the 12 vertices 
//  (3) define the 30 edges in terms of 2 vertices each (by index)
//  (2) define the 20 faces in terms of 3 vertices each (by index)
//
//  see http://en.wikipedia.org/wiki/Icosahedron
//
//            -1       +1
//  axes:   x left  to  right
//          y down  to  up
//          z back  to  front
//
//  φ = (1 + √5) / 2 is the golden ratio
//               x   y   z
//  0, 1, 2, 3  (±1, ±φ, 0)  
//  4, 5, 6, 7  (0, ±1, ±φ)   
//  8, 9,10,11  (±φ, 0, ±1) 
//
//  The radius is sqrt( 1 + φ^2 )
//  we will normalize that to one, later, and rotate so that 
//  there are vertices top and bottom

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

typedef struct { int a, b, c;    } Index3 ;
typedef struct { int a, b;       } Index2 ;
typedef struct { double x,y,z;   } Point  ;

#define  SZ     ( 10*(NM+2)*NM )  // max array size (divis by 4)

Point    vertex[SZ] ;             // element vertices
double     xdis[SZ], ydis[SZ], zdis[SZ] ; // element location
Index3   colrgb[SZ] ;             // element rgb
int         col[SZ] ;             // element libdg color index
int          ap[SZ] ;             // element phase

#define  RMMKDIR    "rm -rf %sdir ; mkdir %sdir"
#define  PNGFMT     "%sdir/a%05d.png"
#define  PNG2SWF    "png2swf -o %s.swf -r %d %sdir/*.png" 

// --------------------------------------------------------------------
// this array is used to generate the icosahedron main vertices
// 2 is φ ,  1 is 1, 0 is 0

// point 0 is at the left, and 2, 4, 6, 8, and 10 are clockwise around it
// point 1 is at the right, and 3, 5, 7, 9, and 11 are clockwise around it
// 0 is opposite 1, 2 is opposite 3, 4 is opposite 5, etc.

static Index3 vertABC[12] = {
               { 2, 1, 0},   //  0 left point after rotation
               {-2,-1, 0},   //  1 right point after rotation
               { 2,-1, 0},   //  2
               {-2, 1, 0},   //  3
               { 1, 0, 2},   //  4
               {-1, 0,-2},   //  5
               { 0, 2, 1},   //  6
               { 0,-2,-1},   //  7
	       { 0, 2,-1},   //  8
               { 0,-2, 1},   //  9
               { 1, 0,-2},   // 10
               {-1, 0, 2},   // 11
};

// --------------------------------------------------------------------
// this array is used to generate the 30 edges, between the numbered 
// vertices in vertABC

static Index2 edge[30] = {
               {  0,  2 },
               {  0,  4 },
               {  0,  6 },
               {  0,  8 },
               {  0, 10 },
               {  1,  3 },
               {  1,  5 },
               {  1,  7 },
               {  1,  9 },
               {  1, 11 },
               {  2,  4 },
               {  2,  7 },
               {  2,  9 },
               {  2, 10 },
               {  3,  5 },
               {  3,  6 },
               {  3,  8 },
               {  3, 11 },
               {  4,  6 },
               {  4,  9 },
               {  4, 11 },
               {  5,  7 },
               {  5,  8 },
               {  5, 10 },
               {  6,  8 },
               {  6, 11 },
               {  7,  9 },
               {  7, 10 },
               {  8, 10 },
	       {  9, 11 },
};

// --------------------------------------------------------------------
// this array is used to generate the 20 faces, between the numbered
// vertices in vertABC

static Index3 face[20] = {
               {  0,  2,  4 },
               {  0,  4,  6 },
               {  0,  6,  8 },
               {  0,  8, 10 },
               {  0, 10,  2 },
               {  1,  3,  5 },
               {  1,  5,  7 },
               {  1,  7,  9 },
               {  1,  9, 11 },
               {  1, 11,  3 },
               {  2,  7, 10 },
               {  5,  7, 10 },
               {  5,  8, 10 },
               {  5,  8,  3 },
               {  6,  8,  3 },
               {  6, 11,  3 },
	       {  6, 11,  4 },
               {  9, 11,  4 },
               {  9,  2,  4 },
               {  9,  2,  7 },
};

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

double  psi   ; 
double  znorm ;

int     sun1, white, red, green, blue     ; //  color indexes 
int     dred, gray, trans, lgray, cyan    ;
int     dgray, yellow, magenta, black     ; 
int     dcyan                             ;

double  xsize   = (double)   XSIZE        ; // screen width
double  xcenter = (double) ( XSIZE / 4 )  ; // center of big side view
int     ex      = (int) (1.5*ORBSIZE)     ; // drawn earth center x
int     ey      = (int) (1.2*ORBSIZE)     ; // drawn earth center y
int     re      = ESIZE / 2               ; // drawn earth radius 
int     nvert                             ; // vertices in geosphere

double     pi2  = 6.2831853070795864      ;
gdImagePtr im                             ; // Declare the image
FILE       *pngout                        ; // Declare output file

// =============================================================

int  normalize( Point *a ) {
   double len = sqrt( a->x*a->x + a->y*a->y + a->z*a->z );
   if( len == 0.0 ) {
      fprintf(stderr, "Normalize error%12.3f%12.3f%12.3f\n", a->x, a->y,a->z );
      return 1 ;
   }
   a->x /= len ;
   a->y /= len ;
   a->z /= len ;
   // printf( "N %12.3f%12.3f%12.3f%12.3f\n",a->x, a->y, a->z, len );
   return  0  ;
} 

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

double   val( int type ) {
   switch( type ) {
      case -2 :  return -psi ;
      case -1 :  return -1.0 ;
      case  0 :  return  0.0 ;
      case  1 :  return  1.0 ;
      case  2 :  return  psi ;
   }
   return 0.0 ;
}

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

int     gcolor[C5MAX]          ; // color map numbers
int     icolor                 ; //
double  mmin, mmid, mmax       ; //
double  cmin, cmid, cmax       ; //
double  am0, am1               ; //
double  al0, al1               ; //


// given magnitude, return color index
int  mcolor( double mm ) {
   if( mm < mmin )  mm = mmin ;    // nothing dimmer than mmin
   int mag = (int) ( am1*mm+am0 ) ;
   if( mm < mmid ) mag = al1*log(mm)+al0 ;
   if( mag > C5MAX-1 ) mag = C5MAX-1  ;    // nothing brighter 
   return gcolor[ mag ] ;
}


// color setup
void colorstart() {
   int    cp, cq,  ix, of ;
   white   =gdImageColorAllocate( im, 255, 255, 255);
   black   =gdImageColorAllocate( im,   0,   0,   0);
   gray    =gdImageColorAllocate( im, 128, 128, 128);
   dgray   =gdImageColorAllocate( im,  96,  96,  96);
   lgray   =gdImageColorAllocate( im, 192, 192, 192);
   cyan    =gdImageColorAllocate( im,   0, 192, 192);
   dcyan   =gdImageColorAllocate( im,   0,  48,  48);
   red     =gdImageColorAllocate( im, 255,   0,   0);
   dred    =gdImageColorAllocate( im, 128,   0,   0);
   yellow  =gdImageColorAllocate( im, 255, 255,   0);
   green   =gdImageColorAllocate( im,   0, 255,   0);
   blue    =gdImageColorAllocate( im,   0,   0, 255);
   magenta =gdImageColorAllocate( im, 192,   0, 192);
   trans   =gdImageColorAllocate( im,   1,   1,   1);

   // array colors ------

   for( ix = 0 ; ix < nvert ; ix++ ) {
     int ir = colrgb[ix].a ;
     int ig = colrgb[ix].b ;
     int ib = colrgb[ix].c ;
     col[ix]       = gdImageColorAllocate( im, ir   , ig   , ib   );
     col[ix+nvert] = gdImageColorAllocate( im, ir/4 , ig/4 , ib/4 );
   }

   // radio colors ------

   // color definitions for radio scale
   // The colors run from
   //  cmin                     cmax
   //  0        64      128     320
   //  black    red     dkgray  white
   //  mmin     
   //  1e-4                     1.0

   int    bgcolor,rcolor             ; // bluegreen and red
   int    i ;
   char   labstring[32];

   mmax   = pow(((double)nvert),2.0)   ; // largest expected magnitd
   mmid   = MMTHR*mmax                 ; // log scale threshold
   mmin   = MMMIN*mmax                 ; // smallest plotted magnitd
   cmax   = (double) C5MAX             ; // largest expected color
   cmin   = 0.0                        ; // smallest expected color
   am1    = (cmax-cmin)/(mmax-mmin)    ; // power slope to color
   am0    = cmin - am1*mmin            ; // power offset to color
   cmid   = am1*mmid + am0             ; // transition to log scale
   al1    = (cmid-cmin)/log(mmid/mmin) ; // log slope to color
   al0    = cmin - al1*log(mmin)       ; // log offset to color

   // allocate gray scale array 
   for( icolor=0 ; icolor<C1MAX ; icolor++ ) {
     rcolor  = icolor ;
     bgcolor = 0 ;
     gcolor[icolor] = gdImageColorAllocate(im, rcolor, bgcolor, bgcolor);
   }
   for(          ; icolor<C2MAX ; icolor++ ) {
     rcolor  = C1MAX  ;
     bgcolor = icolor-C1MAX ;
     gcolor[icolor] = gdImageColorAllocate(im, rcolor, bgcolor, bgcolor);
   }
   for(          ; icolor<C5MAX ; icolor++ ) {
     rcolor  = icolor-C1MAX  ;
     bgcolor = icolor-C1MAX  ;
     gcolor[icolor] = gdImageColorAllocate(im, rcolor, bgcolor, bgcolor);
   }

   // draw color scale 

   for( i = 0 ; i < 6 ; i += 1 ) {
      double mm=mmax*pow(10.0,-i);

      int dxt = 500 + 88*i ;
      int dx0 = dxt + 47 ;
      int dx1 = dx0 + 1  ;
      int dx2 = dx1 + 20 ;
      int dx3 = dx2 + 1  ;

      int dy3 = OY1+65   ;
      int dy2 = dy3 - 1  ;
      int dy1 = dy2 - 20 ;
      int dy0 = dy1 - 1  ;
      int dyt = dy3 - 6  ;

      gdImageRectangle(im,       dx0, dy0, dx3, dy3, gray        );
      gdImageFilledRectangle(im, dx1, dy1, dx2, dy2, mcolor( mm ));
      
      sprintf( labstring, "1e-%1d", i );
      gdImageStringFT( im, NULL,         // imagespace, bounding
          white, FNT, 14, 0.0,           // color, font, fontsize, angle
          dxt, dyt,                      // x, y
          labstring );                   // the text
   }
}
   

// draw arrowhead (either degrees or radians ) -----------------------------

void arrowhead( double xp, double yp, double asize,
                double deg, double rad, int color ) {
   double angl = rad+pi2*deg/360.0 ;
   gdPoint arrow[3] ;

   arrow[0].x = (int)(xp);
   arrow[0].y = (int)(yp);
   arrow[1].x = (int)(xp+asize*(    -cos(angl)-0.5*sin(angl)));
   arrow[1].y = (int)(yp+asize*( 0.5*cos(angl)-    sin(angl)));
   arrow[2].x = (int)(xp+asize*(    -cos(angl)+0.5*sin(angl)));
   arrow[2].y = (int)(yp+asize*(-0.5*cos(angl)-    sin(angl)));
   gdImageFilledPolygon( im, arrow, 3, color );
}

// draw arrow  -------------------------------------------------------------

void arrow( double x0, double y0, double x1, double y1, int color ) {
   double  asize = 0.02*(x1-x0) ;

   gdImageLine( im, (int)x0, (int)y0, (int)x1, (int)y1, color ) ;
   arrowhead( x1, y1, 0.02*(x1-x0), 0.0, 0.0, color ) ;
}

// draw rotating earth and orbiting body and side arrow  -------------

void earth( double ang ) {
   double  x,  y  ;
   double  xs, ys ;
   int  x0, y0 ;
	   int  x1, y1 ;
   double  an0 ;
   int ie    ;
   int orbit = 2*ORBSIZE ;

   gdImageFilledArc(im, ex, ey, ESIZE, ESIZE,   0, 360, dcyan, gdArc);
   gdImageFilledArc(im, ex, ey, ESIZE, ESIZE,  90, 270,  cyan, gdArc);
   gdImageArc(      im, ex, ey, orbit, orbit,  30, 330, white );
   gdImageArc(      im, ex, ey, orbit, orbit, 330,  30, dgray );

   // draw longitude lines  ( 6 = 180/30 degrees )
   for( ie = 0 ; ie < 6 ; ie ++ ) {
      an0 = (ang/6.0) + pi2*( (double)ie) / 12.0 ;
      x0  = ex - (int)( ((double)re) * cos(an0) );
      x1  = ex + (int)( ((double)re) * cos(an0) );
      y0  = ey + (int)( ((double)re) * sin(an0) );
      y1  = ey - (int)( ((double)re) * sin(an0) );
      gdImageLine( im, x0, y0, x1, y1, black );
   }

   // draw latitude lines
   for( ie = 1 ; ie < 3 ; ie++ ) {
      an0 = pi2*( (double) ie ) / 12.0 ;
      x0  = (int)( 2.0 * ((double)re) * cos(an0) );
      gdImageArc( im, ex, ey, x0, x0, 0, 360, black );
   }
   
   // draw arrowhead to show side point of view
   x= -( ORBSIZE + ASCALE + 2.0 ) ;
   y=  0.0 ;
   xs = x*cos(ang)+y*sin(ang) + (double) ex ;
   ys = y*cos(ang)-x*sin(ang) + (double) ey ;
 
   arrowhead( xs, ys, ASCALE,  0.0, -ang, white );
}

// ==========================================================================

int main() {
   double  angle ;
   int     frame ;
   char    pngfilename[64];
   char    command[128];

   double  test   ;   // depth testing
   int     ox, oy ;
   int     s0     ;

   double  xu, yu ;
   int     xi, yi ;
   int     icol   ;
   double   d2r = pi2/360.0      ; // degrees to radians
   int     i      ;
   double  dv = (double) V ;
   char    labstring[40] ;
   
   // double   k   = pi2/WAVELENGTH ; // wavenumber
   // KLUDGE  approximate scaling for 1 meter thinsat spacing

   double  k = 8.0 * V / ( SPACE * WAVELENGTH );

   psi   = 0.5 * ( 1.0 + sqrt(5.0) ) ;
   znorm = sqrt( psi*psi + 1.0 ) ;

   sprintf( command, RMMKDIR, NAME, NAME );
   system( command );

   // COMPUTE THE ICOSAHEDRON -------------------------------------------
   // we will generate the icosahedron, then interpolate the edge points
   // and face points between them.   These will then be normalized to
   // radius 1

   // compute 12 vertices ------------------------- 

   for( nvert = 0 ; nvert < 12 ; nvert++ ) {
      double xtemp = val( vertABC[nvert].a );
      double ytemp = val( vertABC[nvert].b );
      double ztemp = val( vertABC[nvert].c );
      
      // rotate so x value of 0th vertex is max
      double xr = psi*xtemp + ytemp ;
      double yr = psi*ytemp - xtemp ;
      double zr = znorm*ztemp       ;
      double sxz = XZROTATE         ;
      double cxz = sqrt( 1.0-sxz*sxz );

      // rotate x to z (radial) directions, y (NS) remains same
      vertex[nvert].x = cxz * xr - sxz * zr ;
      vertex[nvert].y = yr                  ;
      vertex[nvert].z = cxz * zr + sxz * xr ;
   }
   // nvert is now 12

   // compute 30 edges ---------------------------
   //  using two points on the icosahedron, in the edge table

   for( i = 0 ; i < 30 ; i++ ) {
      Point va = vertex[ edge[i].a ] ;   // find icosahedron points
      Point vb = vertex[ edge[i].b ] ;   // find icosahedron points
      double da ;

      for( da = 1.0 ; da < dv ; da++ ) {
         double db = dv-da ;
         vertex[ nvert ].x = da*va.x + db*vb.x ;
         vertex[ nvert ].y = da*va.y + db*vb.y ;
         vertex[ nvert ].z = da*va.z + db*vb.z ;
         nvert++ ;   // add more vertices
   }  }

   // number of new edge points = 30*(V-1)
   // nvert is now 30*V - 18

   // compute 20 faces ---------------------------
   //  using three points on the icosahedron, in the face table

   for( i = 0 ; i < 20 ; i++ ) {
      Point va = vertex[ face[i].a ] ;   // find icosahedron points
      Point vb = vertex[ face[i].b ] ;   // find icosahedron points
      Point vc = vertex[ face[i].c ] ;   // find icosahedron points
      double da, db ;

      for( da = 1.0 ; da < dv ; da++ ) {
         for( db = 1.0 ; db < ( dv-da) ; db++ ) {
            double dc = dv-(da+db);
            vertex[ nvert ].x = da*va.x + db*vb.x + dc*vc.x ;
            vertex[ nvert ].y = da*va.y + db*vb.y + dc*vc.y ;
            vertex[ nvert ].z = da*va.z + db*vb.z + dc*vc.z ;
            nvert++ ;   // add more vertices
   }  }  }

   // number of new face points = 10*V^2 - 30*V + 20
   // nvert is now 10*V^2 + 2

   // color -----------------------------------------------
   // first we normalize each vertex point to radius 1.  Then
   // we will choose vertex color based on position in the base array

   for( i = 0 ; i < nvert ; i++ ) {
      normalize( &vertex[i] );
      colrgb[i].a = (int) ( CMID + CSCALE*vertex[i].x ) ;
      colrgb[i].b = (int) ( CMID + CSCALE*vertex[i].y ) ;
      colrgb[i].c = (int) ( CMID + CSCALE*vertex[i].z ) ;
   }

   // ICOSAHEDRON COMPUTED --------------------------------------------

   // graphics loop ---------------------------------------------------

   for( frame = 0 ; frame < NPICS ; frame++ ) {
      im = gdImageCreateTrueColor(XSIZE, YSIZE );
      sprintf( command,  PNGFMT, NAME, frame );
      pngout = fopen( command, "wb");

      angle = frame * pi2 / NPICS ;
      colorstart() ;
      earth(angle) ;

      //  labels -----------------------------------------------------------

      gdImageStringFT( im, NULL,                       // imagespace, bounding 
                       white, FNT, FSZ,                // color, font, fontsize
                       0.0, 0.00*XSIZE, 0.95*YSIZE,    // angle, x, y
                       " Side View" );                 // the text

      sprintf( command, "%7d thinsats", nvert );

      gdImageStringFT( im, NULL,                       // imagespace, bounding 
                       white, FNT, FSZ/2,              // color, font, fontsize
                       0.0, 0.00*XSIZE, 0.85*YSIZE,    // angle, x, y
                       command      );                 // the text


      gdImageStringFT( im, NULL,                       // imagespace, bounding
                       white, FNT, FSZ/2,              // color, font, fontsize
                       0.0, 0.26*XSIZE, 2.2*YSCENT,    // angle, x, y
                       "    Top View" );               // the text

      gdImageStringFT( im, NULL,                       // imagespace, bounding
                       white, FNT, FSZ/2,              // color, font, fontsize
                       0.0, 0.50*XSIZE, 2.2*YSCENT,    // angle, x, y
                       " Front View" );                // the text

      gdImageStringFT( im, NULL,                       // imagespace, bounding
                       white, FNT, FSZ/2,              // color, font, fontsize
                       0.0, 0.71*XSIZE, 2.2*YSCENT,    // angle, x, y
                       "   Side View" );               // the text

      //  time code -------------------------------------------------------
 
      double time0 = ( 4.0 * (double) frame ) / ( (double) NPICS ) ;
      char   timestring[64] ;

      sprintf( timestring, "%3.1f/4.0 hours,%6.0fx speedup", time0, SPEEDUP );

      gdImageStringFT( im, NULL,                       // imagespace, bounding 
                       white, FNT, 0.6*FSZ,            // color, font, fontsize
                       0.0, XSIZE-16.0*FSZ,            // angle, x
                       0.95*YSIZE, timestring );       // y, the text 

      //  direction arrows ------------------------------------------------
          // large side view, small top view, small side view

      arrow( 0.00*XSIZE, YBCENT, 2.0*XBCENT, YBCENT, gray ); // large side view
      arrow( 0.25*XSIZE, YSCENT, 0.52*XSIZE, YSCENT, gray );
      arrow( 0.71*XSIZE, YSCENT, 0.95*XSIZE, YSCENT, gray );

      //  compute element locations in X, Y, Z --------------------------
      //  this is where we skew and distort the array.    The Z scaling
      //  becomes the x distance skew.  The orbit of each element can be
      //  elliptical in Y (north-south)  versus Z (radial), while the X
      //  distance (along the orbit) gets the skew 

      double zm = SPACE * ZEXP ;   // radial
      double ym = SPACE * YEXP ;
      double xm = SPACE * XEXP ;
      double xt = SPACE * ( 1.0 + ASKEW * ZEXP ) ; // range of x values
      double as = ASKEW * zm/ym ;
      double dtest = (257.0/1024.0) ;
      double test0 = (-1025.0/1024.0);
      double test1 = ( 1025.0/1024.0);
      double zfactr = ZFACTR/zm ;

      // offsets from small top view to screen
      int  txo = (int) ( (0.26+0.125) * xsize );
      int  tyo = (int) ( YSCENT );

      // offsets from small front view to screen 
      int  fxo = (int) ( (0.48+0.125) * xsize );
      int  fyo = (int) ( YSCENT );

      // offsets from small side view to screen
      int ssxo = (int) ( (0.71+0.125) * xsize );
      int ssyo = (int) ( YSCENT );

      // offsets from large side view to screen
      int lsxo =         XBCENT ;
      int lsyo = (int) ( YBCENT );

      // compute the points of the rotated, scaled, and skewed array
      for( i=0 ; i < nvert ; i++ ) {
         zdis[i] = zm*(cos(angle)*vertex[i].z+sin(angle)*vertex[i].y ) ;
         ydis[i] = ym*(cos(angle)*vertex[i].y-sin(angle)*vertex[i].z ) ;
         xdis[i] = xm*vertex[i].x            + as*ydis[i] ;
      }


      //  layer from back, so that bigger points in front overlap smaller
      //  points in back.   This will make a large loop
      for( test = test0; test < test1 ; test += dtest ) {

         // compute display slices
         double z0 =  test*zm ;
         double z1 = dtest*zm + z0 ;
         double y0 =  test*ym ;
         double y1 = dtest*ym + y0 ;
         double x0 =  test*xt ;
         double x1 = dtest*xt + x0 ;
 
         // step through all the vertices miltiple times
         for( i=0 ; i < nvert; i++ ) {

            //  draw view from top (Z+) --------------------------------------

            if( ( ydis[i] >= y0 ) && ( ydis[i] < y1 ) ) {

               // small top view

	       ox = txo + (int) ( 0.25*xdis[i] );
               oy = tyo + (int) ( 0.25*zdis[i] );
               gdImageSetPixel( im, ox, oy, col[i] );
               // gdImageFilledArc(im, ox, oy, s0, s0, 0, 360, col[i], gdArc);

               // orbiting tiny top view

               xu = -ORBSIZE - OSCALE * zdis[i] ;
               yu =            OSCALE * xdis[i] ; 

               xi = ex + (int)( xu*cos(angle)+yu*sin(angle)) ;
               yi = ey + (int)( yu*cos(angle)-xu*sin(angle)) ;

               if(  ( xi <   ex      ) 
                  ||( yi < ( ey-re ) )
                  ||( yi > ( ey+re ) ) ) icol = col[      i] ;
               else                      icol = col[nvert+i] ;

               gdImageSetPixel( im, xi, yi, icol );
               // gdImageFilledArc( im, xi, yi, 3, 3, 0, 360, icol, gdArc );
            }

            // side views ---------------------------------------------

            if( ( zdis[i] >= z0 ) && ( zdis[i] < z1 ) ) {

               // BIG side view 
               ox = lsxo + (int) (  xdis[i] );
               oy = lsyo + (int) ( -ydis[i] );
               // gdImageSetPixel( im, ox, oy, col[i] );

               s0 = (int) (     SIZE * (1.00 + zfactr * zdis[i] ) ) ;
               gdImageFilledArc(im, ox, oy, s0, s0, 0, 360, col[i], gdArc);

               // small side view  
               ox = ssxo + (int) (  0.25*xdis[i] );
               oy = ssyo + (int) ( -0.25*ydis[i] );
               gdImageSetPixel( im, ox, oy, col[i] );

               // s0 = (int) ( 0.35*SIZE * (1.00 + zfactr * zdis[i] ) ) ;
               // gdImageFilledArc(im, ox, oy, s0, s0, 0, 360, col[i], gdArc);
            }

            // small front view - apogee to left, closer ---------------

            if( ( xdis[i] >= x0 ) && ( xdis[i] < x1 ) ) {
               ox = fxo + (int) (  0.25*zdis[i] );
               oy = fyo + (int) ( -0.25*ydis[i] );
               gdImageSetPixel( im, ox, oy, col[i] );
               // s0 = (int) ( 0.35*SIZE * (1.30 + 2.0* zfactr * zdis[i] ) ) ;
               // gdImageFilledArc(im, ox, oy, s0, s0, 0, 360, col[i], gdArc);
            }
         }
      }

      // ---------------------------------------------------------------------
      // RADIO COMPUTATION LOOP ----------------------------------------------
      // setup
      // 1) construct the target vector from the center of the array,
      //    length 2pi/wavelength
      // 2) For each element, construct the vector from the center of
      // the array, and compute the dot product of that vector with 
      // the target vector to compute the conjugate phase of the element
      //
      // For each plotted angular point:
      //
      // 3) construct the target vector from the center of the array
      // to the plotted angular point, length 2pi/wavelength
      // 4) Repeat 2 for each element, adding to the phase difference
      // 5) compute the sin and cos of the phase difference, add to
      //    the running total for the angular point
      // 6) compute the magnitude, choose color from palette,
      //    plot point
      // 7) Draw box and key

      double ox0    = (double) OX0               ; // output display left
      double ox1    = (double) OX1               ; // output display right
      double oy0    = (double) OY0               ; // output display top
      double oy1    = (double) OY1               ; // output display bottom
      double dangy  = (double) DANGY             ; // display vert angle offset
      double dangx  = (double) DANGX             ; // display horz angle offset
      double dangw  = (double) DANGW             ; // display width kilometers
      double dangh  = (double) DANGH             ; // display height kilometers
      double dangsc = (double) DANGSCALE         ; // kilometer distance
      double k2r    = 1.0/dangsc                 ; // radians per kilometer
      double dangx0 = k2r*(dangx-0.5*dangw)      ; // left pixel value, radians
      double dangx1 = k2r*(dangx+0.5*dangw)      ; // right pix value, radians
      double dangxs = k2r*dangw/(ox1-ox0)        ; // pixel X slope radians
      double dangxf = dangx0-dangxs*ox0          ; // pixel X offset radians
      double dangy0 = k2r*(dangy-0.5*dangh)      ; // top pixel value, radians
      double dangy1 = k2r*(dangy+0.5*dangh)      ; // bottom pix value, radians
      double dangys = k2r*dangh/(oy1-oy0)        ; // pixel Y slope radians
      double dangyf = dangy0-dangys*oy0          ; // pixel Y offset radians
      int    oxc    = (OX0+OX1)/2                ; // center of output display
      int    oyc    = (OY0+OY1)/2                ;

      // make wavevector 
      double kx = -k * sin( d2r*TANGX )                  ;
      double ky = -k * cos( d2r*TANGX ) * sin( d2r*TANGY );
      double kz = -k * cos( d2r*TANGX ) * cos( d2r*TANGY );

      // conjugate phase for each element

      for( i=0 ; i < nvert; i++ ) {
         ap[i] = kx*xdis[i] + ky*ydis[i] + kz*zdis[i] ;
      }

      // --------------------------------------------------------------------
      // Big Computation Loop -----------------------------------------------
      // a bunch of points, a bunch of emitters
      // The outer ox and oy loops are relatively small, though I can imagine
      // image sizes of 20K by 20K pixels.  
      // The inner nvert loop is thousands to millions,
      // and parallelizing that would save much time.

      int pct  = -1 ;
      int pct0 =  0 ;

      for( ox=OX0 ; ox<OX1 ; ox++ ) {
         double  angx = dangxf + dangxs*((double)ox );
         double  spx  = sin( angx );
         double  cpx  = cos( angx );
   
         pct0  = pct ;
         pct   = (100*(ox-OX0))/(OX1-OX0) ;
   
         if( pct != pct0 ) {
            printf( "%04d/%04d%4d\r", frame, NPICS, pct ) ;
            fflush(stdout);
         }
   
         for( oy=OY0 ; oy<OY1 ; oy++ ) {
            double angy =  dangyf + dangys*((double)oy );
            double spy  = sin( angy );
            double cpy  = cos( angy );
   
            // sum the components
            double s2 = 0.0 ;
            double c2 = 0.0 ;

            kx = k * spx ;
            ky = k * cpx * spy ;
            kz = k * cpx * cpy ;
            
            // CENTRAL COMPUTING LOOP
            // element loop - spends huge time in this loop, nvert is 
            // typically thousands to millions
            for( i=0 ; i < nvert; i++ ) {
               double phase = ap[i] + kx*xdis[i] + ky*ydis[i] + kz*zdis[i] ;
               s2    += sin( phase ) ;
               c2    += cos( phase ) ;
            }
            // END CENTRAL COMPUTING LOOP
           
            // power/intensity computation
            double mm  = s2*s2 + c2*c2 ;
            gdImageSetPixel( im, ox, oy, mcolor(mm) ) ;
      }  } // END big computation triple loop

      // draw box  ----------------------------------------------------------
      
      gdImageLine( im, OX0, OY0, OX1, OY0, white );
      gdImageLine( im, OX1, OY1, OX1, OY0, white );
      gdImageLine( im, OX1, OY1, OX0, OY1, white );
      gdImageLine( im, OX0, OY0, OX0, OY1, white );
     
      // draw ticmarks on edges ----------------------------------------------
      
      double dangs = (double) DANGS                  ; // display angle step
      int    mx    = (int) (0.001+0.5*(dangw/dangs)) ; // tickmark +/- count
      int    my    = (int) (0.001+0.5*(dangh/dangs)) ; // tickmark +/- count
      double sx    = (ox1-ox0)*dangs/dangw           ; // tickmark step size
      double sy    = (oy0-oy1)*dangs/dangh           ; // tickmark step size
   
      // draw x tickmarks 

      int  ft     = 16 ;
   
      for( i = -mx ; i <= mx ; i++ ) {
         ox = oxc + (int) ( sx*i ) ;
         gdImageLine( im, ox, OY1, ox, OY1-5, white );
         gdImageLine( im, ox, OY0, ox, OY0+5, white );
   
         sprintf( labstring, "%3d", ( (int)( i*dangs ) ) );
         gdImageStringFT( im, NULL,        // imagespace, bounding
            white, FNT, ft, 0.0,           // color, font, fontsize, angle
            ox-1.6*ft, OY1+1.5*ft,         // angle, x, y
            labstring );                   // the text
      }
      gdImageLine( im, oxc, OY1, oxc, OY1-10, white );  // big tics
      gdImageLine( im, oxc, OY0, oxc, OY0+10, white );

            // draw y tickmarks
      
      for( i = -mx ; i <= mx ; i++ ) {
         oy = oyc + (int) ( sy*i ) ;
         gdImageLine( im, OX1, oy,  OX1-5, oy, white );
         gdImageLine( im, OX0, oy,  OX0+5, oy, white );
   
         sprintf( labstring, "%3d", ( (int)( i*dangs ) ) );
         gdImageStringFT( im, NULL,       // imagespace, bounding
            white, FNT, ft, 0.0,          // color, font, fontsize, angle
            OX0-3*ft, oy+ft/2,            // x, y,
            labstring );                  // the text
   
      }
      gdImageLine( im, OX0, oyc, OX0+10, oyc, white );
      gdImageLine( im, OX1, oyc, OX1-10, oyc, white );

      gdImageStringFT( im, NULL,         // imagespace, bounding
          white, FNT, ft+4, 0.0,         // color, font, fontsize, angle
          OX0+1*ft, OY1-1*ft,            // x, y
          "Ground spot, km");            // the text
   
      // Output the frame in PNG format.----------------------------

      gdImagePngEx( im, pngout, 1 );
      gdImageDestroy(im);
      fclose(pngout);
   }

   sprintf( command, PNG2SWF, NAME, RATE, NAME );
   system( command );
   return 0 ;
}