// radio array, 5x5x5, rotating
//
// compile with        cc -o a3d04 a3d04.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>

#define  RMMKDIR     "rm -rf a3d04dir ; mkdir a3d04dir"
#define  PNG2SWF     "png2swf -o a3d04.swf -r 5 a3d04dir/*.png" 
#define  PNGFMT      "a3d04dir/a%04d.png"

#define  YCENT       130
#define  WAVELENGTH  10.0         // radio wavelength in pixel units
#define  SPACE       30.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        10           // size of drawn spot
#define  ASIZE       5            // number of elements per side
#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                ; 
   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] ;
   int      black, sun1, white, red, green, blue, gray, dgray, trans ;
   int      gcolor[256];
   double   ax[ASIZE][ASIZE][ASIZE];            // x position of element
   double   ay[ASIZE][ASIZE][ASIZE];            // y position of element
   double   az[ASIZE][ASIZE][ASIZE];            // z position of element
   double   ap[ASIZE][ASIZE][ASIZE];            // 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=0.5*(1.0-ASIZE) ;
   int      xdraw, ydraw ;
   int      icolor ;
   int      cir    ;
   double   stilt             ;
   double   ctilt             ;
   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 = ASIZE*ASIZE*ASIZE ;
   double   norm = 255.0 /asize3 ;
   double   phase             ;
   int      xm                ;  // tickmark max
   int      xs                ;  // tickmark scale
   int      frame             ;
   int      pct=0             ;  // percentage of complete frame
   int      pct0=0            ;

#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

//  set up target directory

   system( RMMKDIR );

//  outer loop, rotate tilt 1/4 circle --------------------------------
//  swf will continue looping, making apparent full circle
  
   for( frame=0 ; frame<NPICS ; frame++ ) {
      tilt  = 0.25 * pi2 * ((double) frame)/((double) NPICS ) ;
      stilt = sin( tilt );
      ctilt = cos( tilt );

#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 ; }
            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 ) {
         printf( "%04d/%04d%4d\r", frame, NPICS, 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 
   xm = (int) (0.001+0.5*(DANGW/DANGS)) ; 
   xs = (XSIZE-XLEFT)*DANGS/DANGW;
   for( ix = -xm ; ix <= xm ; ix++ ) {
      ox = oxc + (int) ( xs*ix ) ;
      gdImageLine( im, ox, ybot, ox, ybot-5, white );

      sprintf( labstring, "%3d", ( (int)( ix*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-10, white );

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

      sprintf( labstring, "%3d", ( (int)( ix*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-10, 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

   system( PNG2SWF );
}