// population coverage
// overlay population versus latitude with coverage graph
// assume round earth and circular orbit
//
// gcc -o pc02 pc02.c -lgd -lpng -lm

#define  INFILE    "population1K.png"
#define  OFILE     "popcoverage1K.png"

#define  ELEV        5.0   // elevation above horizon, degrees
#define  RE       6371.0   // earth average radius, kilometers
#define  R288    12789.0   // orbit radius

#define  PI        3.1415926535897932
#define  RAD       (PI/180.00)
#define  FNT       "DejaVuMonoSans"
#define  FS0        30             // title
#define  FS1        22
#define  FS2        16             // smaller

// MAJOR KLUDGE WARNING!  These plot corners depend on original drawing
#define  X0          4
#define  X1       1010
#define  X2       1023
#define  Y0         62
#define  Y1        565
#define  Y2        622
#define  XS0      (X0)
#define  XSCAL     350

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

//-----------------------------------------------------------------
int main() {

   FILE    *datfile                     ; //
   char    gnuplot[200]                 ; //
   char    filename[80]                 ; //

   gdImagePtr im                        ; // working image

   FILE       *pngin;                   ; // output file
   FILE       *pngout;                  ; // output file
   int        black, white, blue, purp  ; // colors
   int        ix, iy                    ; // daytime point
   int        ixo, iyo                  ; // old daytime point
   int        mx                        ; // midnight point
   int        mxo                       ; // old midnight point
   double     elev0 = acos( RE/R288 )   ; //
   double     elev  = elev0 - RAD*ELEV  ; // 
   double     celev = cos( elev )       ; //
   double     midloss = asin( RE/R288 ) ; //

   double     yscale = PI/(Y1-Y0)       ; //
   double     ymid   = 0.5*(Y1+Y0)      ; // middle of graph

   printf( "%9.4f=midloss", midloss/RAD );
   printf( "%9.4f=elev0 %9.4f=elev %9.6f=celev\n", elev0/RAD,  elev/RAD, celev );

   int        y0     = (int) ( ymid     - elev/yscale );
   int        y1     = (int) ( ymid + 1 + elev/yscale );

   pngin = fopen( INFILE, "rb" );
   im = gdImageCreateFromPng(pngin);
   fclose( pngin );
   black   = gdImageColorAllocate (im,   0,   0,   0);
   white   = gdImageColorAllocate (im, 255, 255, 255);
   blue    = gdImageColorAllocate (im,   0,   0, 255);
   purp    = gdImageColorAllocate (im, 151,   0, 152);

   iyo     = y0  ;
   ixo     = XS0 ;
   mxo     = XS0 ;

   gdImageSetThickness( im, 3 );

   for( iy = y0+1 ; iy <=y1 ; iy++ ) {
      double y        = (double) iy ;
      double latitude = yscale*( ymid - iy );   // latitude in radians
      double clat     = cos( latitude );
      double coverage = 0.0 ;
      if( clat > celev ) { coverage = acos( celev / clat ) ; }
      ix = XS0 + (int) ( XSCAL * coverage );
      gdImageLine(im,  ix, iy, ixo, iyo,  blue );

      mx = XS0 ;
      if( coverage > midloss ) {
         mx = XS0 + (int) ( XSCAL * (coverage-midloss) );
      }
      gdImageLine(im,  mx, iy, mxo, iyo,  black );

      ixo = ix ;  iyo = iy ; mxo = mx ;
   }

   gdImageFilledRectangle( im, 0, 0,    X2, Y0-2, white );
   gdImageFilledRectangle( im, 0, Y1+2, X2, Y2,   white );

   gdImageStringFT( im, NULL,                // imagespace, bounding box
                    black, FNT, FS0, 0.0,    // color, font, fontsize, angle
                    X0+40,  Y0-15,           // x, y,
                    "Population and Coverage by Latitude" ); // text

   gdImageStringFT( im, NULL,                // imagespace, bounding box
                    purp, FNT, FS1, 0.0,     // color, font, fontsize, angle
                    X0+XSCAL, Y1-140,        // x, y,
                    "population by latitude" ); // text

   gdImageStringFT( im, NULL,                // imagespace, bounding box
                    blue, FNT, FS1, 0.0,     // color, font, fontsize, angle
                    X0+XSCAL, Y1-100,        // x, y,
                    "6am-6pm coverage by latitude" ); // text

   gdImageStringFT( im, NULL,                // imagespace, bounding box
                    black, FNT, FS2, 0.0,    // color, font, fontsize, angle
                    X0+XSCAL, Y1-70,         // x, y,
                    "linearly decreasing to" ); // text

   gdImageStringFT( im, NULL,                // imagespace, bounding box
                    black, FNT, FS1, 0.0,    // color, font, fontsize, angle
                    X0+XSCAL, Y1-40,         // x, y,
                    "midnight coverage by latitude" ); // text

   gdImageStringFT( im, NULL,                // imagespace, bounding box
                    black, FNT, FS2, 0.0,    // color, font, fontsize, angle
                    X0+20, Y1+23,            // x, y,
                    "assume: 5% elevation above horizon, circular orbit, M288=12789km radius" ); // text
   gdImageStringFT( im, NULL,                // imagespace, bounding box

                    black, FNT, FS2, 0.0,    // color, font, fontsize, angle
                    X0+20, Y1+51,            // x, y,
                    "y2K population graph from http://www.radicalcartography.net/index.html?histpop" ); // text

   pngout = fopen( OFILE, "wb");
   gdImagePngEx( im, pngout, 1 );
   gdImageDestroy(im);
   fclose(pngout);

   return(0);
}