// radio array, 5x5x5, sine warp spacing, power scale
//
// compile with:      gcc -o a6 a6.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
//  
// You will need truetype fonts.  If you don't have them, you can
// copy ../fonts/truetype/.. from openoffice.org to /usr/share/
// libgd origin at upper left

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

#define  PNGFMT      "a6.png"

//#define  VSPACE      10.0         // sin variation
#define  VSPACE      0.0          // sin variation

#define  WAVELENGTH  10.0         // radio wavelength in pixel units
#define  SPACE       40.0         // array spacing in pixel units

#define  TANGX       0.46891338   // target X angle in radians
#define  TANGY       0.0          // target Y angle in radians
#define  DANGX       0.0          // output display center X angle in degrees
#define  DANGY       0.0          // output display center Y angle in degrees
#define  DANGH       70.0         // output display angle height in degrees
#define  DANGW       70.0         // output display angle width  in degrees
#define  DANGS       10.0         // output display angle step
#define  DANGE       30.0         // output earth display angle

#define  ZSCALE      400.0        // depth scaling
#define  TILT        0.125663706  // array tilt around X axis, radians
// #define  TILT        0.00         // array tilt around X axis, radians
#define  SPOT        10           // size of drawn spot
#define  ASIZE       5            // number of elements per side
#define  ASIZE3      125          // number of elements

#define  XSIZE       1280         // display window in pixels
#define  YSIZE       720          // display window in pixels
#define  OX0         600          // output array left
#define  OX1         1250         // output array right
#define  OY0          11          // output array top
#define  OY1         661          // output array bottom
#define  YACENT      360          // center of array object display
#define  XACENT      300          // center of array object display
#define  FNT         "DejaVuMonoSans"
#define  FSZ         28

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

int main() {

   double   pi2 = 8.0*atan(1.0)  ;
   double   d2r = pi2/360.0      ; // degrees to radians
   double   k   = pi2/WAVELENGTH ; // wavenumber
   FILE     *pngout              ; // file handle for PNG output frame
   char     labstring[80]        ; // used for labelling
   double   ax[ASIZE3]           ; // x position of element
   double   ay[ASIZE3]           ; // y position of element
   double   az[ASIZE3]           ; // z position of element
   double   ap[ASIZE3]           ; // phase of element
   double   rx[ASIZE3]           ; // randomized position of element
   double   ry[ASIZE3]           ; // randomized position of element
   double   rz[ASIZE3]           ; // randomized position of element
   int      icolor               ; // index for computing colors
   int      i                    ; // general counter
   int      fs                   ; // large font size

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

//  dither value generation ----------------------------------------------

   int jcnt=0, jx, jy, jz ;
   double radx, rady, radz, sinx, siny, sinz, cosx, cosy, cosz ;
   double rads = 8.0*atan(1.0) / ((double) ASIZE );
   for( jz=0; jz<ASIZE ; jz++ ) {
      radz = rads*(0.5+(double)jz) ;
      cosz = cos( radz );
      sinz = sin( radz );
      for( jy=0; jy<ASIZE ; jy++ ) {
      rady = rads*(0.5+(double)jy) ;
         cosy = cos( rady );
         siny = sin( rady );
         for( jx=0; jx<ASIZE ; jx++ ) {
            radx = rads*(0.5+(double)jx) ;
            cosx = cos( radx );
            sinx = sin( radx );
            rx[jcnt] = VSPACE * ( siny - cosz );
            ry[jcnt] = VSPACE * ( sinz - cosx );
            rz[jcnt] = VSPACE * ( sinx - cosy );
            jcnt++ ;
   } } } 

   double ffrac = 1.0 ;
   double rfactor = 0.5*(1.0-cos(0.5*pi2*ffrac));

//  define array with even spacing ------------------------------------------

   double   stx    = sin( (double) TANGX ) ;
   double   ctx    = cos( (double) TANGX ) ;
   double   sty    = sin( (double) TANGY ) ;
   double   cty    = cos( (double) TANGY ) ;
   double   tilt   = (double) TILT         ;
   double   stilt  = sin(          tilt  ) ;
   double   ctilt  = cos(          tilt  ) ;
   double   asize  = (double) ASIZE        ;  
   double   offset = 0.5*(1.0-asize)       ;
   double   space  = SPACE                 ;
   int      ix, iy, iz                     ;

   int      acnt = 0                       ;
   for( iz=0; iz<ASIZE ; iz++ ) {
      double zz = space*(((double)iz)+offset) ;
      for( iy=0; iy<ASIZE ; iy++ ) {
         double yy = space*(((double)iy)+offset) ;
         double zr = zz*ctilt - yy*stilt;
         double yr = yy*ctilt + zz*stilt;
         double xr = -2.0*yr ;         // xr offset location moves with orbit
         for( ix=0; ix<ASIZE ; ix++ ) {
            // compute array positions 
            double xx = xr + space*(((double)ix)+offset) ;
            xx += rfactor*rx[acnt] ;
            yr += rfactor*ry[acnt] ;
            zr += rfactor*rz[acnt] ;

            ax[acnt] = xx ;
            ay[acnt] = yr ;
            az[acnt] = zr ;

            // compute target angle
            // first rotation, Y axis, by TANGX
            // double bx = xx*ctx - zr*stx ;
            // double by = yr ;
            // double bz = zr*ctx + xx*stx ;

            // second rotation, X axis, by TANGY 

            // double cy = by*cty - bz*sty ;
            // double cz = bz*cty + by*sty ;

            // ap[i] = -k*cz    ;
            ap[acnt] = -k*( (zr*ctx+xx*stx)*cty + yr*sty )  ;
            acnt++ ;
   }  }  }


   // set up array and define colors -----------------------------------------

   gdImagePtr im = gdImageCreateTrueColor(XSIZE, YSIZE );

   // Allocate standard colors
   int black = gdImageColorAllocate(im,   0,   0,   0);
   int white = gdImageColorAllocate(im, 255, 255, 255);
   int sun1  = gdImageColorAllocate(im,  51,  51, 102);
   int red   = gdImageColorAllocate(im, 255,   0,   0);
   int green = gdImageColorAllocate(im,   0, 255,   0);
   int dgreen= gdImageColorAllocate(im,   0, 128,   0);
   int blue  = gdImageColorAllocate(im,   0,   0, 255);
   int gray  = gdImageColorAllocate(im, 128, 128, 128);
   int dgray = gdImageColorAllocate(im,  48,  48,  48);
   int trans = gdImageColorAllocate(im, 1, 1, 1);

#define C5MAX  256 
#define C2MAX  128
#define C1MAX   64
#define MMTHR  0.4

   gdImageSetThickness( im, 3 );

   int gcolor[C5MAX]              ; // color map numbers

   // allocate gray scale array 
   for( icolor=0 ; icolor<C5MAX ; icolor++ ) {
      int col = 255-icolor ;
      gcolor[icolor] = gdImageColorAllocate(im, col, col, col);
   }

   // White background ??

   gdImageFilledRectangle( im, 0, 0, XSIZE-1, YSIZE-1, white );

   // Slide Label -----------------------------------------------------------
  
   fs= 40 ;
   gdImageStringFT( im, NULL,           // imagespace, bounding box
      black , FNT, fs, 0.0,             // color, font, fontsize, angle
      20     , fs+25  ,                 // x, y
      "Thinsat Array" );                // text

   fs= 30 ;
   gdImageStringFT( im, NULL,           // imagespace, bounding box
      black , FNT, fs, 0.0,             // color, font, fontsize, angle
      50     , 120    ,                 // x, y
      "Regular Spacing" );              // text
//    "Sin Dithered Spacing" );         // text

   gdImageStringFT( im, NULL,           // imagespace, bounding box
//    dgreen , FNT, fs, 0.0,            // color, font, fontsize, angle
      red , FNT, fs, 0.0,            // color, font, fontsize, angle
      60     , 540    ,                 // x, y
//    "Dither washes\nout concentrated\ngrating lobes" );         // text
      "Uniform array\nmakes concentrated\ngrating lobes" );       // text

//  draw array with perspective ---------------------------------------------

   for( i=0; i<ASIZE3 ; i++ ) {
      double zscale = 0.8 / ( 1.0- az[i]/((double) ZSCALE ) ) ;
      int cir       =          (int) ( zscale*SPOT  ) ;
      int xdraw     = XACENT + (int) ( zscale*ax[i] ) ;
      int ydraw     = YACENT + (int) ( zscale*ay[i] ) ;
      icolor        = 128 + (int)(16.0*az[i]/space )  ;
      if( icolor > 255 ) { icolor = 255 ; }
      gdImageFilledEllipse(im, xdraw, ydraw, cir, cir, gcolor[icolor] );
   }  

// set up scaling

   double mmax   = pow(asize,6.0)             ; // largest expected magnitd
   double mmid   = MMTHR*mmax                 ; // log scale threshold
   double mmin   = 1e-4*mmax                  ; // smallest plotted magnitd
   double cmax   = (double) C5MAX             ; // largest expected color
   double cmin   = 0.0                        ; // smallest expected color
   double am1    = (cmax-cmin)/(mmax-mmin)    ; // power slope to color
   double am0    = cmin - am1*mmin            ; // power offset to color
   double cmid   = am1*mmid + am0             ; // transition to log scale
   double al1    = (cmid-cmin)/log(mmid/mmin) ; // log slope to color
   double al0    = cmin - al1*log(mmin)       ; // log offset to color
   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 degrees
   double dangh  = (double) DANGH             ; // display height degrees
   double dangx0 = d2r*(dangx-0.5*dangw)      ; // left pixel value, radians
   double dangx1 = d2r*(dangx+0.5*dangw)      ; // right pix value, radians
   double dangxs = d2r*dangw/(ox1-ox0)        ; // pixel X slope radians
   double dangxf = dangx0-dangxs*ox0          ; // pixel X offset radians
   double dangy0 = d2r*(dangy-0.5*dangh)      ; // top pixel value, radians
   double dangy1 = d2r*(dangy+0.5*dangh)      ; // bottom pix value, radians
   double dangys = d2r*dangh/(oy1-oy0)        ; // pixel Y slope radians
   double dangyf = dangy0-dangys*oy0          ; // pixel Y offset radians
   int    ox, oy                              ; // output grid indices
   int    oxc    = (OX0+OX1)/2                ; // center of output display
   int    oyc    = (OY0+OY1)/2                ;

//   now compute phases

   for( ox=OX0 ; ox<OX1 ; ox++ ) {
      double  angx = dangxf + dangxs*((double)ox );
      double  spx  = sin( angx );
      double  cpx  = cos( angx );

      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 ;

         // compute z distance and phase of each source to plane
         for( i=0; i<ASIZE3 ; i++ ) {
            // xr = ax[i] ;
            // yr = ax[i] ;
            // zr = az[i] ;

            // compute pixel angle  -----
            // first rotation, Y axis, by TANGX
            //   bx = xr*cpx - zr*spx ;
            //   by = yr ;
            //   bz = zr*cpx + xr*spx ;
            // second rotation, X axis, by TANGY 
            //   cx = bx ;
            //   cy = by*cpy - bz*spy ;
            //   cz = bz*cpy + by*spy ;

            // ap[i] = -k*cz    ;
            // ap[i] = -k*( (zr*cpx+xr*spx)*cpy + yr*spy);
   
            double phase = k*( ay[i]*spy + ( az[i]*cpx + ax[i]*spx)*cpy ) ;
            phase += ap[i]        ;
            s2    += sin( phase ) ;
            c2    += cos( phase ) ;
         } 

         // power/intensity computation
         double mm  = s2*s2 + c2*c2 ;
         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 ) mag = C5MAX  ;    // nothing brighter 

         gdImageSetPixel( im, ox, oy, gcolor[mag] ) ;
   }  }

   // draw earth circle ----------------------------------------------------

   int ex = (2*DANGE*(OX1-OX0))/DANGW ;
   int ey = (2*DANGE*(OY1-OY0))/DANGH ;

   gdImageArc( im, oxc, oyc, ex, ey, 0, 360, green );
 
   // draw box  ------------------------------------------------------------
   
   gdImageLine( im, OX0, OY0, OX1, OY0, black );
   gdImageLine( im, OX1, OY1, OX1, OY0, black );
   gdImageLine( im, OX1, OY1, OX0, OY1, black );
   gdImageLine( im, OX0, OY0, OX0, OY1, black );
  
   // 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 

   for( i = -mx ; i <= mx ; i++ ) {
      ox = oxc + (int) ( sx*i ) ;
      gdImageLine( im, ox, OY1, ox, OY1-5, black );
      gdImageLine( im, ox, OY0, ox, OY0+5, black );

      sprintf( labstring, "%3d", ( (int)( i*dangs ) ) );
      gdImageStringFT( im, NULL,        // imagespace, bounding
         black, FNT, FSZ, 0.0,          // color, font, fontsize, angle
         ox-1.6*FSZ, OY1+1.5*FSZ,       // angle, x, y
         labstring );                   // the text
   }
   gdImageLine( im, oxc, OY1, oxc, OY1-10, black );
   gdImageLine( im, oxc, OY0, oxc, OY0+10, black );

   // draw y tickmarks
   
   for( i = -mx ; i <= mx ; i++ ) {
      oy = oyc + (int) ( sy*i ) ;
      gdImageLine( im, OX1, oy,  OX1-5, oy, black );
      gdImageLine( im, OX0, oy,  OX0+5, oy, black );

      sprintf( labstring, "%3d", ( (int)( i*dangs ) ) );
      gdImageStringFT( im, NULL,       // imagespace, bounding
         black, FNT, FSZ, 0.0,         // color, font, fontsize, angle
         OX0-3*FSZ, oy+FSZ/2,          // x, y,
         labstring );                  // the text

   }
   gdImageLine( im, OX0, oyc, OX0+10, oyc, black );
   gdImageLine( im, OX1, oyc, OX1-10, oyc, black );

   // labels on output array ------------------------------------------------
   gdImageStringFT( im, NULL,          // imagespace, bounding box
      dgreen , FNT, FSZ, 0.0,           // color, font, fontsize, angle
      OX0+10, OY1-10,                  // x, y
      "Ground Pattern" );              // text

   // output the frame ------------------------------------------------------

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