attachment:pc02.c of Brown2013

   1 // population coverage
   2 // overlay population versus latitude with coverage graph
   3 // assume round earth and circular orbit
   4 //
   5 // gcc -o pc02 pc02.c -lgd -lpng -lm
   6 
   7 #define  INFILE    "population1K.png"
   8 #define  OFILE     "popcoverage1K.png"
   9 
  10 #define  ELEV        5.0   // elevation above horizon, degrees
  11 #define  RE       6371.0   // earth average radius, kilometers
  12 #define  R288    12789.0   // orbit radius
  13 
  14 #define  PI        3.1415926535897932
  15 #define  RAD       (PI/180.00)
  16 #define  FNT       "DejaVuMonoSans"
  17 #define  FS0        30             // title
  18 #define  FS1        22
  19 #define  FS2        16             // smaller
  20 
  21 // MAJOR KLUDGE WARNING!  These plot corners depend on original drawing
  22 #define  X0          4
  23 #define  X1       1010
  24 #define  X2       1023
  25 #define  Y0         62
  26 #define  Y1        565
  27 #define  Y2        622
  28 #define  XS0      (X0)
  29 #define  XSCAL     350
  30 
  31 #include <gd.h>
  32 #include <math.h>
  33 #include <string.h>
  34 #include <stdio.h>
  35 #include <stdlib.h>
  36 
  37 //-----------------------------------------------------------------
  38 int main() {
  39 
  40    FILE    *datfile                     ; //
  41    char    gnuplot[200]                 ; //
  42    char    filename[80]                 ; //
  43 
  44    gdImagePtr im                        ; // working image
  45 
  46    FILE       *pngin;                   ; // output file
  47    FILE       *pngout;                  ; // output file
  48    int        black, white, blue, purp  ; // colors
  49    int        ix, iy                    ; // daytime point
  50    int        ixo, iyo                  ; // old daytime point
  51    int        mx                        ; // midnight point
  52    int        mxo                       ; // old midnight point
  53    double     elev0 = acos( RE/R288 )   ; //
  54    double     elev  = elev0 - RAD*ELEV  ; // 
  55    double     celev = cos( elev )       ; //
  56    double     midloss = asin( RE/R288 ) ; //
  57 
  58    double     yscale = PI/(Y1-Y0)       ; //
  59    double     ymid   = 0.5*(Y1+Y0)      ; // middle of graph
  60 
  61    printf( "%9.4f=midloss", midloss/RAD );
  62    printf( "%9.4f=elev0 %9.4f=elev %9.6f=celev\n", elev0/RAD,  elev/RAD, celev );
  63 
  64    int        y0     = (int) ( ymid     - elev/yscale );
  65    int        y1     = (int) ( ymid + 1 + elev/yscale );
  66 
  67    pngin = fopen( INFILE, "rb" );
  68    im = gdImageCreateFromPng(pngin);
  69    fclose( pngin );
  70    black   = gdImageColorAllocate (im,   0,   0,   0);
  71    white   = gdImageColorAllocate (im, 255, 255, 255);
  72    blue    = gdImageColorAllocate (im,   0,   0, 255);
  73    purp    = gdImageColorAllocate (im, 151,   0, 152);
  74 
  75    iyo     = y0  ;
  76    ixo     = XS0 ;
  77    mxo     = XS0 ;
  78 
  79    gdImageSetThickness( im, 3 );
  80 
  81    for( iy = y0+1 ; iy <=y1 ; iy++ ) {
  82       double y        = (double) iy ;
  83       double latitude = yscale*( ymid - iy );   // latitude in radians
  84       double clat     = cos( latitude );
  85       double coverage = 0.0 ;
  86       if( clat > celev ) { coverage = acos( celev / clat ) ; }
  87       ix = XS0 + (int) ( XSCAL * coverage );
  88       gdImageLine(im,  ix, iy, ixo, iyo,  blue );
  89 
  90       mx = XS0 ;
  91       if( coverage > midloss ) {
  92          mx = XS0 + (int) ( XSCAL * (coverage-midloss) );
  93       }
  94       gdImageLine(im,  mx, iy, mxo, iyo,  black );
  95 
  96       ixo = ix ;  iyo = iy ; mxo = mx ;
  97    }
  98 
  99    gdImageFilledRectangle( im, 0, 0,    X2, Y0-2, white );
 100    gdImageFilledRectangle( im, 0, Y1+2, X2, Y2,   white );
 101 
 102    gdImageStringFT( im, NULL,                // imagespace, bounding box
 103                     black, FNT, FS0, 0.0,    // color, font, fontsize, angle
 104                     X0+40,  Y0-15,           // x, y,
 105                     "Population and Coverage by Latitude" ); // text
 106 
 107    gdImageStringFT( im, NULL,                // imagespace, bounding box
 108                     purp, FNT, FS1, 0.0,     // color, font, fontsize, angle
 109                     X0+XSCAL, Y1-140,        // x, y,
 110                     "population by latitude" ); // text
 111 
 112    gdImageStringFT( im, NULL,                // imagespace, bounding box
 113                     blue, FNT, FS1, 0.0,     // color, font, fontsize, angle
 114                     X0+XSCAL, Y1-100,        // x, y,
 115                     "6am-6pm coverage by latitude" ); // text
 116 
 117    gdImageStringFT( im, NULL,                // imagespace, bounding box
 118                     black, FNT, FS2, 0.0,    // color, font, fontsize, angle
 119                     X0+XSCAL, Y1-70,         // x, y,
 120                     "linearly decreasing to" ); // text
 121 
 122    gdImageStringFT( im, NULL,                // imagespace, bounding box
 123                     black, FNT, FS1, 0.0,    // color, font, fontsize, angle
 124                     X0+XSCAL, Y1-40,         // x, y,
 125                     "midnight coverage by latitude" ); // text
 126 
 127    gdImageStringFT( im, NULL,                // imagespace, bounding box
 128                     black, FNT, FS2, 0.0,    // color, font, fontsize, angle
 129                     X0+20, Y1+23,            // x, y,
 130                     "assume: 5% elevation above horizon, circular orbit, M288=12789km radius" ); // text
 131    gdImageStringFT( im, NULL,                // imagespace, bounding box
 132 
 133                     black, FNT, FS2, 0.0,    // color, font, fontsize, angle
 134                     X0+20, Y1+51,            // x, y,
 135                     "y2K population graph from http://www.radicalcartography.net/index.html?histpop" ); // text
 136 
 137    pngout = fopen( OFILE, "wb");
 138    gdImagePngEx( im, pngout, 1 );
 139    gdImageDestroy(im);
 140    fclose(pngout);
 141 
 142    return(0);
 143 }