// radio array, grating lobes versus size of array
//
// compile with:  cc -o a3d03 a3d03.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 swftools to build a swf movie from individual png images
//        for more information, see http://www.swftools.org/

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

#define  RMMKDIR     "rm -rf a3d03dir ; mkdir a3d03dir"
#define  PNG2SWF     "png2swf -o a3d03.swf -r 2 a3d03dir/*.png" 
#define  PNGFMT      "a3d03dir/a%04d.png"

#define  YCENT       130
#define  WAVELENGTH  2.33         // radio wavelength in pixel units
#define  SPACE       6.0          // array spacing in pixel units
#define  NPICS       100          // number of frames

#define  TANGX       0.0          // target X angle in degrees
#define  TANGY       0.0          // target Y angle in degrees
#define  DANGX       0.0          // display center X angle in degrees
#define  DANGY       0.0          // display center Y angle in degrees

#define  DANGW       150.0        // display angle width in degrees
#define  DANGS       10.0         // display angle step

#define  ZSCALE      400.0        // depth scaling
// #define  TILT        0.04         // array tilt around X axis, radians
#define  TILT        0.00         // array tilt around X axis, radians
#define  SPOT        4            // size of drawn spot
#define  ASIZE       8            //
#define  ASMIN       2            //
#define  ASMAX       20           //
#define  XLEFT       40           // left side pixel count
#define  XSIZE       1000         // display window in pixels
#define  YSIZE       750          // display window in pixels
#define  FNT         "DejaVuMonoSans"
#define  FSZ         10

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

int main() {

   gdImagePtr im                ; 
   int      npics   = ASMAX     ;
   int      asize   = ASIZE     ;
   double   space   = SPACE     ;
   double   xleft   = XLEFT     ;
   double   yacent  = YCENT     ;
   double   xacent  = 0.5*(XSIZE+XLEFT) ;
   int      ybot    = 2* (int)(yacent) ;
   double   angH    = DANGW*(YSIZE-ybot)/(XSIZE-XLEFT);
   double   pi2     = 8.0 * atan( 1.0 );
   double   deg2rad = pi2/360.0 ;

   int      oxc = (XSIZE+XLEFT)/2 ;
   int      oyc = (YSIZE+ybot)/2 ;

   double   dangx                                ; // pixel X angle radians
   double   dangx0  = deg2rad*(DANGX-0.5*DANGW)  ; // left pixel value, radians
   double   dangx1  = deg2rad*(DANGX+0.5*DANGW)  ; // right pixel value, radians
   double   dangxs  = deg2rad*DANGW/(XSIZE-XLEFT); // pixel X slope radians
   double   dangxf  = dangx0-dangxs*XLEFT        ; // pixel X offset radians
   double   dangy                                ; // pixel Y angle
   double   dangy0  = deg2rad*(DANGY+0.5*angH)   ; // top pixel value, radians
   double   dangy1  = deg2rad*(DANGY-0.5*angH)   ; // bottom pixel value, radians
   double   dangys  = deg2rad*angH/(YSIZE-ybot)  ; // pixel Y slope radians
   double   dangyf  = dangy0-dangys*YSIZE        ; // pixel Y offset radians
   double   k       = pi2/WAVELENGTH             ; // wavenumber
   double   tilt    = TILT ;
   FILE     *pngout       ;
   char     dirname[80]   ;
   char     framename[80] ;
   char     labstring[80] ;
   char     link0[80]     ;
   char     link1[80]     ;
   int      black, sun1, white, red, green, blue, gray, dgray, trans ;
   int      gcolor[256];
   double   ax[ASMAX][ASMAX][ASMAX];            // x position of element
   double   ay[ASMAX][ASMAX][ASMAX];            // y position of element
   double   az[ASMAX][ASMAX][ASMAX];            // z position of element
   double   ap[ASMAX][ASMAX][ASMAX];            // phase of element
   double   xx, yy, zz ;
   double   xr, yr, zr ;
   double   bx, by, bz ;    // coordinates rotated on y axis
   double   cx, cy, cz ;    // coordinates rotated on x axis
   int      ix, iy, iz ;
   int      ox, oy     ;    // output grid
   double   zscale ;
   double   offset ;
   int      xdraw, ydraw ;
   int      icolor ;
   int      cir    ;
   double   stilt = sin( tilt );
   double   ctilt = cos( tilt );
   int      mag               ;  // pixel magnitude, 0-255 ;
   double   ss, cc            ;  //summed sin and cos components
   double   stx = sin(TANGX ) ;
   double   ctx = cos(TANGX ) ;
   double   sty = sin(TANGY ) ;
   double   cty = cos(TANGY ) ;
   double   spx               ;  // pixel angle sin and cos
   double   cpx               ;
   double   spy               ;
   double   cpy               ;
   double   asize3 = 0        ;
   double   asize3old=0       ;
   double   norm = 255.0 /asize3 ;
   double   phase             ;
   int      m                 ;  // tickmark max
   int      s                 ;  // tickmark scale 
   int      i                 ;  // general counter
   int      frame             ;
   int      pct=0             ;  // percentage of complete frame
   int      pct0=0            ;

   int      prog0=0           ;
   int      prog1=0           ;
   int      progmax=0         ;
   double   progfrac          ;
   double   pctmult           ;

#ifdef DEBUG
   printf( "%12.6f%12.6f%12.6f%15.9f\n",  dangx0, dangx1, dangxf, dangxs );
   printf( "%12.6f%12.6f%12.6f%15.9f\n",  dangy0, dangy1, dangyf, dangys );
#endif

//  estimate progress
   for( i=ASMIN ; i <= ASMAX ; i++ ) { progmax += i*i*i ; }

//  set up target directory

   system( RMMKDIR );

//  outer loop, sinusodial spacing change --------------------------------
//  from SPACEMIN to SPACE ( then link and fill in )
  
   for( frame=ASMIN ; frame<=ASMAX ; frame++ ) {
      asize = frame ;
      asize3old = asize3     ;
      asize3 = asize*asize*asize ;
      offset = 0.5*(1.0-asize) ;
      norm   = 255.0 /asize3 ;
      prog0  = prog1 ;
      prog1  = prog0 + (int) asize3 ;

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


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

   for( iz=0; iz<asize ; iz++ ) {
      zz = space*(((double)iz)+offset );
      for( iy=0; iy<asize ; iy++ ) {
         yy = space*(((double)iy)+offset) ;
         zr = zz*ctilt - yy*stilt;
         yr = yy*ctilt + zz*stilt;
         xr = -2.0*yr ;
         for( ix=0; ix<asize ; ix++ ) {
            // compute array positions 
            xx = xr + space*(((double)ix)+offset) ;
            ax[ix][iy][iz] = xx ;
            ay[ix][iy][iz] = yr ;
	    az[ix][iy][iz] = zr ;

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

            // second rotation, X axis, by TANGY 

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

            // ap[ix][iy][iz] = -k*cz    ;
            ap[ix][iy][iz] = -k*( (zr*ctx+xx*stx)*cty + yr*sty )  ;

         }
#ifdef DEBUG
         printf( "%3d%3d%18.3f%18.3f%18.3f%18.3f\n", iy, iz,
	    ax[2][iy][iz], ay[2][iy][iz], az[2][iy][iz], ap[2][iy][iz] );
#endif
   }  }
#ifdef SINGLE
   printf( "array setup done\n");
#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, 128, 128, 128);
   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);
   }

#ifdef SINGLE
   printf( "color setup done\n");
#endif

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

   for( iz=0; iz<asize ; iz++ ) {
      for( iy=0; iy<asize ; iy++ ) {
         for( ix=0; ix<asize ; ix++ ) {
            zscale = 1.0 / ( 1.0- az[ix][iy][iz]/ZSCALE ) ;
            cir    = (int) ( zscale * SPOT );
            xdraw  = (int) ( xacent + zscale*ax[ix][iy][iz] ) ;
            ydraw  = (int) ( yacent + zscale*ay[ix][iy][iz] ) ;
            icolor = 180 + (int)(30.0*az[ix][iy][iz]/space ) ;
            if( icolor > 255 ) { icolor = 255 ; }
            if( icolor <   0 ) { icolor =   0 ; }
            if( ( ydraw > 0 ) && ( ydraw < ybot ) ) { 
              gdImageFilledEllipse(im, xdraw, ydraw, cir, cir, gcolor[icolor] );
            }
   }  }  }

#ifdef SINGLE
   printf( "array drawn\n");
#endif

//   now compute phases

   for( ox=XLEFT ; ox<XSIZE ; ox++ ) {
      dangx =  dangxf + dangxs*((double)ox );
      spx   = sin( dangx );
      cpx   = cos( dangx );
 
      pct0  = pct ;
      pct   = (100*(ox-XLEFT))/(XSIZE-XLEFT) ;

#ifndef SINGLE
      if( pct != pct0 ) {
         pctmult   = ((double)(prog1-prog0))*((double)pct);
         progfrac  = (100.0*((double)prog0)+pctmult)/((double)progmax) ;
                 
         printf( "%03d/%03d%6.2f%3d\r", frame, npics, progfrac, pct ) ;
         fflush(stdout);
      }
#endif

      for( oy=ybot ; oy<YSIZE ; oy++ ) {
         dangy =  dangyf + dangys*((double)oy );
         spy   = sin( dangy );
         cpy   = cos( dangy );

         // sum the components
         ss = 0.0 ;
         cc = 0.0 ;

#ifdef DEBUG
         if( ( ox==oxc ) && ( oy==oyc ) ) {
              printf( "%4d%5d%12.6f%12.6f%12.6f%12.6f\n\n", ox, oy,
                              cpx, spx, cpy, spy );
         }
#endif

         // compute z distance and phase of each source to plane
         for( iz=0; iz<asize ; iz++ ) {
            for( iy=0; iy<asize ; iy++ ) {
               for( ix=0; ix<asize ; ix++ ) {

                  // xr = ax[ix][iy][iz] ;
                  // yr = ax[ix][iy][iz] ;
                  // zr = az[ix][iy][iz] ;

                  // 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[ix][iy][iz] = -k*cz    ;
                  // ap[ix][iy][iz] = -k*( (zr*cpx+xr*spx)*cpy + yr*spy);
   
                  phase = k*( ay[ix][iy][iz]*spy +
                          ( az[ix][iy][iz]*cpx + ax[ix][iy][iz]*spx)*cpy ) ;
                  ss += sin( phase + ap[ix][iy][iz] );
                  cc += cos( phase + ap[ix][iy][iz] );

#ifdef DEBUG
                  if( ( ox==oxc ) && ( oy==oyc ) ) {
                     printf( "%3d%3d%3d%12.6f%12.6f%12.6f%12.6f\n", ix,iy,iz,
                              phase, ap[ix][iy][iz], ss, cc );
                  }
#endif
               }
            }
         } 

         mag = (int) (norm*sqrt( ss*ss + cc*cc ));
         if( abs(mag-181) < 10 ) {
            gdImageSetPixel( im, ox, oy, red         ) ;
         }
         else {
            gdImageSetPixel( im, ox, oy, gcolor[mag] ) ;
         }
   }  }

#ifdef SINGLE
   printf( "intensity points drawn\n");
#endif

   // draw tickmarks on edges ----------------------------------------------

   // draw x tickmarks 
   m = (int) (0.001+0.5*(DANGW/DANGS)) ; 
   s = (XSIZE-XLEFT)*DANGS/DANGW;
   for( i = -m ; i <= m ; i++ ) {
      ox = oxc + (int) ( s*i ) ;
      gdImageLine( im, ox, ybot, ox, ybot-5, white );

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

   }
   ox = (int) ( xacent ) ;
   gdImageLine( im, xacent, ybot, ox, ybot-20, white );

   // draw y tickmarks
   m = (int) (0.001+0.5*(angH/DANGS)) ; 
   for( i = -m ; i <= m ; i++ ) {
      oy = oyc + (int) ( s*i ) ;
      gdImageLine( im, xleft, oy,  xleft-5, oy, white );

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

   }
   ox = (int) ( xacent ) ;
   gdImageLine( im, xleft, oyc, xleft-20, oyc, white );

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

   sprintf( framename, PNGFMT , frame );
   pngout = fopen( framename, "wb");
   gdImagePngEx( im, pngout, 1 );
   gdImageDestroy(im);
   fclose(pngout);

}  // end of single frame

   // now, hard link some frames to existing frames to "fill in the picture"
   for( i=0 ; i < ASMIN ; i++ ) {
      sprintf( link0 , PNGFMT, ASMIN );
      sprintf( link1 , PNGFMT, i );
      link( link0, link1 ) ;
      sprintf( link0 , PNGFMT, ASMAX );
      sprintf( link1 , PNGFMT, ASMAX+1+i );
      link( link0, link1 ) ;
   }

   // fill in the reverse
   for( i=0 ; i <= ASMIN+ASMAX ; i++ ) {
      sprintf( link0 , PNGFMT, i );
      sprintf( link1 , PNGFMT, (1+2*(ASMIN+ASMAX))-i );
      link( link0, link1 ) ;
   }

   system( PNG2SWF );
}