// uniformly illuminated parabolic dish
//
// OLD compile with:       cc -o adish01A adish01A.c -lgd -lpng -lm -lgsl -lgslcblas
// NEW compile with:       cc -o adish01B adish01B.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/
//
// Uses apngasm to make APNG animation from frames
//
//OLD Uses the GNU Scientific Library to make the J1() bessel function
//OLD for more information, see:
//OLD  http://www.gnu.org/software/gsl/manual/html_node/Regular-Cylindrical-Bessel-Functions.html
//
// now uses j1 function and math.h
//
// This constructs 30 images then dups them forward and reverse.

#include "gd.h"
#include "math.h"
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
// #include <gsl/gsl_sf_bessel.h>

#define  RMMKDIR     "rm -rf adish01dir ; mkdir adish01dir"
#define  PNG2APNG    "/usr/local/bin/apngasm adish01B.png adish01dir/* 1 10" 
#define  PNGFMT      "adish01dir/a%04d.png"

#define  YCENT       130
#define  WAVELENGTH  50.0         // radio wavelength in pixel units
#define  SCALE       10.0         // disk scaling in pixel units
#define  SCMIN       30.0         // minimum dish scale 
#define  SCMAX       60.0         // minimum dish scale 
#define  NPICS       30           // ( half the ) number of frames
#define  RSPACE      30.0         // random variation

#define  DANGW       150          // display angle width in degrees
#define  DANGS       10           // display angle tic step in degrees

#define  XLEFT       40           // left side pixel count
#define  XSIZE       1000         // display window in pixels
#define  YSIZE       750          // display window in pixels
#define  FNT         "/usr/share/fonts/truetype/dejavu/DejaVuSansMono.ttf"
#define  FSZ         10

// #define  SINGLE	     1

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

// int main( ch) {
int main(int argc, char *argv[], char * envp[] ) {

   gdImagePtr im                    ; // pixel image map used by gd
   gdPoint  beam[4]                 ; // trapezoidal output beam
   double   lambda=(double)WAVELENGTH;//
   double   scale  = (double) SCALE ; // dish scale
   double   xleft  = (double) XLEFT ; // left of pixel array
   double   yacent = (double) YCENT ; // center of upper drawing
   double   xsize  = (double) XSIZE ; // width of pixels
   double   ysize  = (double) YSIZE ; // height of pixels
   double   dangw  = (double) DANGW ; // degree width of output pixel array
   double   dangs  = (double) DANGS ; // display angle tic step in degrees
   double   dangw2 = 0.5*dangw      ; // half width of output pixel array
   double   xwide  = xsize-xleft    ; // pixel width of output pixel array
   int      ybot   = 2*YCENT        ; // top of output pixel array
   double   dybot  = (double) ybot  ; // real ""
   double   ytall  = ysize-dybot    ; // height of output pixel array
   double   pi2    = 8.0*atan(1.0)  ; // two times PI
   double   d2r    = pi2/360.0      ; // degrees to radians conversion factor
   int      oxc    =(XSIZE+XLEFT)/2 ; // center of output pixel array
   int      oyc    =(YSIZE+ybot)/2  ; // center of output pixel array
   double   xacent = (double) oxc   ; // real ""
   double   aspect =ytall/xwide     ; // aspect ration of output pixel array
   double   angH   =dangw*aspect    ; // angular height of output pixel array
   double   angH2 = 0.5*angH        ; // half height of output pixel array

////   printf( "ybot%10d  xwide%10d\n", ybot, dybot);  /// temp

   double   angx                    ; // pixel X angle radians
   double   angx0=-d2r*dangw2       ; // left pixel value, radians
   double   angxs=d2r*dangw/xwide   ; // pixel X slope radians per pixel
   double   angxf=angx0-angxs*xleft ; // pixel X offset radians

   double   angy                    ; // pixel Y angle
   double   angy0=-d2r*angH2        ; // top pixel value, radians
   double   angys=d2r*angH/ytall    ; // pixel Y slope radians
   double   angyf=angy0-angys*dybot ; // pixel Y offset radians

   double   ang                     ; // angle from center

   int      brect[8]                ; // this returns a bounding box, not needed, but ...

   FILE     *pngout                 ; // file handle for PNG output frame

   char     dirname[80]             ; // frame output directory
   char     framename[80]           ; // directory/name of frame
   char     labstring[80]           ; // used for labellin
   char     link0[80]               ; // directory/name of source 
   char     link1[80]               ; // directory/name of target hardlink
   int      black, sun1, white, red ; // color map numbers
   int      green, blue, gray       ; // color map numbers
   int      dgray, trans            ; // color map numbers
   int      bcolor[256]             ; // color map numbers
   int      gcolor[256]             ; // color map numbers
   int      ox, oy                  ; // output pixel position
   int      i                       ; // general counter
   double   factr                   ; // beamwidth factor
   int      fs                      ; // big label font size
   int      cir                     ; //
   int      icolor                  ; //
   int      jcolor                  ; //
   int      beamw                   ; //
   int      mag                     ; // pixel magnitude, 0-255
   int      m                       ; // tickmark max
   int      s                       ; // tickmark scale
   int      frame                   ; // animation frame count
   int      pct=0                   ; // percentage of complete frame
   int      pct0=0                  ; // previous percentage
   double   min = 1e-6              ; // min value for bessel function
   double   bx = min                ; // bessel argument
   double   j=j1(bx)/bx             ; // bessel/x near zero value
   double   norm = 255.1/(j*j)      ; // normalizes amplitude
   double   j0=j                    ; //
   double   ffrac                   ; // frame scaling of variable, 0 to 1

#ifdef DEBUG
   printf( "%12.6f\n", j0 );
   printf( "%12.6f%12.6f%15.9f\n",  angx0, angxf, angxs );
   printf( "%12.6f%12.6f%15.9f\n",  angy0, angyf, angys );

//  print environment

    for (i=0 ; envp[i] != NULL; i++) printf("\n%s", envp[i]);
    getchar();
    printf( "\n");
#endif

//  set up target directory

   system( RMMKDIR );

//  outer loop, sinusodial random factor change --------------------------
//  from 0 to 1 to 0
  
   for( frame=0 ; frame<NPICS ; frame++ ) {

      ffrac = 0.5*(1.0-cos(0.5*pi2*((double)frame)/((double)(NPICS-1))));
      scale = SCMIN+(SCMAX-SCMIN)*ffrac ;

#ifndef SINGLE
      printf( "%04d/%04d\r", frame, NPICS ) ;
      fflush(stdout);
#endif

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

      im = gdImageCreateTrueColor(XSIZE, YSIZE );
   
      // Allocate standard colors
      black = gdImageColorAllocate(im,   0,   0,   0);
      white = gdImageColorAllocate(im, 255, 255, 255);
      sun1  = gdImageColorAllocate(im,  51,  51, 102);
      red   = gdImageColorAllocate(im, 255,   0,   0);
      green = gdImageColorAllocate(im,   0, 255,   0);
      blue  = gdImageColorAllocate(im,   0,   0, 255);
      gray  = gdImageColorAllocate(im, 160, 160, 160);
      dgray = gdImageColorAllocate(im,  48,  48,  48);
      trans = gdImageColorAllocate(im,   1,   1,   1);
   
      // allocate gray scale array 
      for( icolor=0 ; icolor<256 ; icolor++ ) {
        gcolor[icolor] = gdImageColorAllocate(im, icolor,   icolor,   icolor);
        bcolor[icolor] = gdImageColorAllocate(im, icolor/3, icolor, icolor/3);
      }
   
   #ifdef SINGLE
      printf( "color setup done\n");
   #endif
   
//   labels -----------------------------------------------------------------
//ERROR labels not appearing ... why?
//
      fs= ybot/6 ;
   
      gdImageStringTTF( im, brect,              // imagespace, bounding box
                       white , FNT, fs, 0.0,   // color, font, fontsize, angle
                       15     , fs+15  ,       //  x, y
                       "Dish Antenna" );       //  text
   
      ox     = oxc + 2*cir +4     ;
      oy     = (int) (0.4*yacent) ;
  
////  printf( "ox%9d      oy%9d    fs%9d\n", ox, oy, fs ); /// TEMP

      fs     = 10*ybot/95 ;
      icolor = gcolor[ (int)( 255.1*(1.0-ffrac )) ] ;
      jcolor = bcolor[ (int)( 255.1*(1.0-ffrac )) ] ;
   
      gdImageStringTTF( im, brect,              // imagespace, bounding box
                       jcolor, FNT, fs, 0.0,   // color, font, fontsize, angle
                       xsize-10*fs, 2*fs+6,     // x, y
                       "Small dish" );         //  text
   
      gdImageStringTTF( im, brect,              // imagespace, bounding box
                       icolor, FNT, fs, 0.0,   // color, font, fontsize, angle
                       xsize-12*fs, 4*fs,      // x, y
                       "Big beam angle" );     //  text
   
      icolor = gcolor[ (int)( 255.1*ffrac ) ] ;
      jcolor = bcolor[ (int)( 255.1*ffrac ) ] ;
   
      gdImageStringTTF( im, brect,              // imagespace, bounding box
                       jcolor, FNT, fs, 0.0,   // color, font, fontsize, angle
                       xsize-10*fs, 6*fs+6,     // x, y
                       "Big dish" );           //  text
     
      //// gdImageLine( im, xsize-10*fs-4, 6*fs-4, ox, oy, jcolor ); // arrow
   
      gdImageStringTTF( im, brect,              // imagespace, bounding box
                       icolor, FNT, fs, 0.0,   // color, font, fontsize, angle
                       xsize-(27*fs/2), 8*fs,  // x, y
                       "Small beam angle");    //  text

//   draw beam --------------------------------------------------------------
   
      // beam 
      for( i=0; i<256; i++) {
   
         factr     = sqrt( (255.0-(double)i)/255.0);
   
         beamw     = (int) (factr*(0.8*scale+0.5*yacent*lambda/scale)) ;
         oy        = (int) (1.7*yacent) ;
         beam[2].x = oxc+beamw ;
         beam[2].y = oy        ;
         beam[3].x = oxc-beamw ;
         beam[3].y = oy        ;
   
         beamw     = (int) ( 0.8 * factr * scale ) ;
         oy        = (int) (0.4*yacent) ;
         beam[0].x = oxc-beamw ;
         beam[0].y = oy        ;
         beam[1].x = oxc+beamw ;
         beam[1].y = oy        ;
         gdImageFilledPolygon( im, beam, 4, gcolor[i] );
      }
 
      oy        = (int) (0.4*yacent) ;
      cir       = (int) (0.39*scale)  ;
   
      // disk
      gdImageFilledEllipse( im, oxc, oy,     4*cir,  cir,   green );
   
      //  square feed 
      beam[0].x = oxc+2     ;
      beam[0].y = oy+cir-1  ;
      beam[1].x = oxc-2     ;
      beam[1].y = oy+cir-1  ;
      beam[2].x = oxc-2     ;
      beam[2].y = oy+cir+5  ;
      beam[3].x = oxc+2     ;
      beam[3].y = oy+cir+5  ;
   
      gdImageFilledPolygon( im, beam, 4,             green );
   
      // feed struts
      gdImageLine( im, oxc-2*cir, oy, oxc,  oy+cir,  green );
      gdImageLine( im, oxc+2*cir, oy, oxc,  oy+cir,  green );
      gdImageLine( im, oxc,       oy, oxc,  oy+cir,  green );
   
//   now compute pixels -----------------------------------------------------
   
      for( ox=XLEFT ; ox<XSIZE ; ox++ ) {
         angx =  angxf + angxs*((double)ox );
   
#ifndef SINGLE
         if( pct != pct0 ) {
            printf( "%04d/%04d%4d\r", frame, NPICS, pct ) ;
            fflush(stdout);
         }
#endif
   
         for( oy=ybot ; oy<YSIZE ; oy++ ) {
            angy = angyf + angys*((double)oy );
            ang  = sqrt( angx*angx + angy*angy) ;
   
            bx = pi2*scale*sin(ang)/lambda  ; // bessel argument
   
            if( fabs(bx) < min ) { j=j0 ; }
            else {
               // OLD j = gsl_sf_bessel_J1(bx)/bx  ; // bessel/x near zero value
               j = j1(bx)/bx                ; // bessel/x near zero value
            }
            mag = (int) (norm*j*j)          ; 
   
#ifdef DEBUG
            if( (ox==oxc) && (oy==oyc) ) {
               printf( "%4d%4d%10.6f%10.6f%10.6f%10.3e%10.3e%4d\n",
                         ox, oy, angx, angy, ang, bx,  j,    mag );
            }
#endif
            gdImageSetPixel( im, ox, oy, gcolor[mag] ) ;
      }  }
   
#ifdef SINGLE
      printf( "intensity points drawn\n");
#endif
   
// draw tickmarks on edges ----------------------------------------------
//
   
////  printf( "oxc%9d    oyc%9d\n", oxc, oyc );  // TEMP

      // draw x tickmarks 
      m = (int) ( 0.001+ dangw2/dangs ); 
      s = (int) ( xwide*dangs/dangw   );

////  printf( "ZSM   s%9d    m%9d\n" , s, m);  // TEMP


      for( i = -m ; i <= m ; i++ ) {
         ox = oxc + s*i ;
         gdImageLine( im, ox, ybot, ox, ybot-5, white );
   
         sprintf( labstring, "%3d", i*DANGS );
         gdImageStringTTF( im, brect,                  // imagespace, bounding
                          white, FNT, FSZ,            // color, font, fontsize
                          0.0, ox+FSZ/2, ybot-FSZ/2,  // angle, x, y
                          labstring );                // the text
////      printf( "AAA FSZ%4d     ox%9d    ybot%9d  %s\n", FSZ, ox, ybot, labstring);  // TEMP
      }
      ox = (int) ( xacent ) ;
      gdImageLine( im, xacent, ybot, ox, ybot-10, white );
   
      // draw y tickmarks
      m = (int) ( 0.001 + angH2/dangs ) ; 
   
      for( i = -m ; i <= m ; i++ ) {
         oy = oyc + s*i ;
         gdImageLine( im, XLEFT, oy,  XLEFT-5, oy, white );
   
         sprintf( labstring, "%3d", i*DANGS );
         gdImageStringTTF( im, brect,                  // imagespace, bounding
                          white, FNT, FSZ,            // color, font, fontsize
                          0.0, XLEFT-4*FSZ, oy+FSZ/2, // angle, x, y
                          labstring );                // the text
   
////      printf( "BBB FSZ%4d  XLEFT%9d    oy  %9d  %s\n", FSZ, XLEFT, oy, labstring);  // TEMP
      }
      ox = (int) ( xacent ) ;
      gdImageLine( im, XLEFT, oyc, XLEFT-10, oyc, white );
   
// labels on array -------------------------------------------------------
      fs= ybot/8 ;
   
      gdImageStringTTF( im, brect,              // imagespace, bounding box
                       white , FNT, fs, 0.0,   // color, font, fontsize, angle
                       XLEFT+15, YSIZE-15,     //  x, y
                       "Ground Pattern" );     //  text
/////      printf( "fs%4d    XLEFT+15%9d    YSIZE-15%9d\n", fs, XLEFT+15, YSIZE-15);  // TEMP

// output the frame ------------------------------------------------------
   
      sprintf( framename, PNGFMT , frame );
      pngout = fopen( framename, "wb");
      gdImagePngEx( im, pngout, 1 );
   
// output mirror frame ---------------------------------------------------
      
      sprintf( framename, PNGFMT , (2*NPICS-frame-1) );
      pngout = fopen( framename, "wb");
      gdImagePngEx( im, pngout, 1 );
   
      gdImageDestroy(im);
      fclose(pngout);
   }

   printf( "frames complete\n");

   // system( "ls -l adish01dir" );
   
   printf( "\n\napng\n");
   printf( "%s\n", PNG2APNG );
   system( PNG2APNG );
}