attachment:a3d05.c of RadioV01

   1 // radio array, 5x5x5, random spacing
   2 //
   3 // compile with    cc -o a3d05 a3d05.c -lgd -lpng -lm
   4 //
   5 // Uses the libgd library.  For documentation see the file:
   6 //      /usr/share/doc/gd-*/index.html  
   7 //   or the website:  http://www.libgd.org/Reference
   8 //  
   9 // You will need truetype fonts.  If you don't have them, you can
  10 // copy ../fonts/truetype/.. from openoffice.org to /usr/share/
  11 //
  12 // Uses swftools to build a swf movie from individual png images
  13 //        for more information, see http://www.swftools.org/
  14 //
  15 // This constructs half the slides, then mirrors them
  16 
  17 #include "gd.h"
  18 #include "math.h"
  19 #include <stdio.h>
  20 #include <stdlib.h>
  21 #include <unistd.h>
  22 
  23 #define  RMMKDIR     "rm -rf a3d05dir ; mkdir a3d05dir"
  24 #define  PNG2SWF     "png2swf -o a3d05.swf -r 10 a3d05dir/*.png" 
  25 #define  PNGFMT      "a3d05dir/a%04d.png"
  26 
  27 #define  YCENT       130
  28 #define  WAVELENGTH  10.0         // radio wavelength in pixel units
  29 #define  SPACE       40.0         // array spacing in pixel units
  30 #define  NPICS       50           // ( half the ) number of frames
  31 #define  RSPACE      30.0         // random variation
  32 
  33 #define  TANGX       0.0          // target X angle in degrees
  34 #define  TANGY       0.0          // target Y angle in degrees
  35 #define  DANGX       0.0          // display center X angle in degrees
  36 #define  DANGY       0.0          // display center Y angle in degrees
  37 
  38 #define  DANGW       150.0        // display angle width in degrees
  39 #define  DANGS       10.0         // display angle step
  40 
  41 #define  ZSCALE      400.0        // depth scaling
  42 // #define  TILT        0.04         // array tilt around X axis, radians
  43 #define  TILT        0.00         // array tilt around X axis, radians
  44 #define  SPOT        10           // size of drawn spot
  45 #define  ASIZE       5            // number of elements per side
  46 #define  XLEFT       40           // left side pixel count
  47 #define  XSIZE       1000         // display window in pixels
  48 #define  YSIZE       750          // display window in pixels
  49 #define  FNT         "DejaVuMonoSans"
  50 #define  FSZ         10
  51 
  52 // ==========================================================================
  53 
  54 int main() {
  55 
  56    gdImagePtr im                    ; // pixel image map used by gd
  57    double   space  = (double) SPACE ; // nominal spacing between server-sats
  58    double   xleft  = (double) XLEFT ; // left of pixel array
  59    double   yacent = (double) YCENT ; // center of upper drawing
  60    double   xsize  = (double) XSIZE ; // width of pixels
  61    double   ysize  = (double) YSIZE ; // height of pixels
  62    double   dangw  = (double) DANGW ; // degree width of output pixel array
  63    double   dangx  = (double) DANGX ; // degree center of output pixel array
  64    double   dangy  = (double) DANGY ; // degree center of output pixel array
  65    double   dangw2 = 0.5*dangw      ; // half width of output pixel array
  66    double   xwide  = xsize-xleft    ; // pixel width of output pixel array
  67    int      ybot   = 2*YCENT        ; // top of output pixel array
  68    double   dybot  = (double) ybot  ; // real ""
  69    double   ytall  = ysize-dybot    ; // height of output pixel array
  70    double   pi2    = 8.0*atan(1.0)  ; // two times PI
  71    double   d2r    = pi2/360.0      ; // degrees to radians conversion factor
  72    int      oxc    =(XSIZE+XLEFT)/2 ; // center of output pixel array
  73    int      oyc    =(YSIZE+ybot)/2  ; // center of output pixel array
  74    double   xacent = (double) oxc   ; // real ""
  75    double   aspect =ytall/xwide     ; // aspect ration of output pixel array
  76    double   angH   =dangw*aspect    ; // angular height of output pixel array
  77    double   angH2 = 0.5*angH        ; // half height of output pixel array
  78 
  79    double   angx                    ; // pixel X angle radians
  80    double   angx0=d2r*(dangx-dangw2); // left pixel value, radians
  81    double   angxs=d2r*dangw/xwide   ; // pixel X slope radians per pixel
  82    double   angxf=angx0-angxs*xleft ; // pixel X offset radians
  83 
  84    double   angy                    ; // pixel Y angle
  85    double   angy0=d2r*(dangy+angH2) ; // top pixel value, radians
  86    double   angys=d2r*angH/ytall    ; // pixel Y slope radians
  87    double   angyf=angy0-angys*ysize ; // pixel Y offset radians
  88 
  89    double   k     = pi2/WAVELENGTH  ; // wavenumber
  90    double   rspace= RSPACE          ; // random spacing amount
  91    double   tilt  = TILT            ; // tilt of array (around Y) in radians
  92    FILE     *pngout                 ; // file handle for PNG output frame
  93    char     dirname[80]             ; // frame output directory
  94    char     framename[80]           ; // directory/name of frame
  95    char     labstring[80]           ; // used for labellin
  96    char     link0[80]               ; // directory/name of source 
  97    char     link1[80]               ; // directory/name of target hardlink
  98    int      black, sun1, white, red ; // color map numbers
  99    int      green, blue, gray       ; // color map numbers
 100    int      dgray, trans            ; // color map numbers
 101    int      gcolor[256]             ; // color map numbers
 102    double   ax[ASIZE][ASIZE][ASIZE] ; // x position of element
 103    double   ay[ASIZE][ASIZE][ASIZE] ; // y position of element
 104    double   az[ASIZE][ASIZE][ASIZE] ; // z position of element
 105    double   rx[ASIZE][ASIZE][ASIZE] ; // x variation of element
 106    double   ry[ASIZE][ASIZE][ASIZE] ; // y variation of element
 107    double   rz[ASIZE][ASIZE][ASIZE] ; // z variation of element
 108    double   ap[ASIZE][ASIZE][ASIZE] ; // phase of element
 109    double   xx, yy, zz              ; //
 110    double   xr, yr, zr              ; //
 111    double   bx, by, bz              ; // coordinates rotated on y axis
 112    double   cx, cy, cz              ; // coordinates rotated on x axis
 113    int      ix, iy, iz              ; // integer counters
 114    int      ox, oy                  ; // output pixel position
 115    int      i                       ; // general counter
 116    double   zscale                  ; // scaling for perspective depth
 117    int      icolor                  ; //
 118    int      cir                     ; //
 119    double   stilt = sin(tilt)       ; // sin of rotation (around y) of array
 120    double   ctilt = cos(tilt)       ; // cos of rotation (around y) of array
 121    int      mag                     ; // pixel magnitude, 0-255
 122    double   ss, cc                  ; // summed sin and cos components
 123    double   stx = sin(TANGX )       ; // sin of X target angle
 124    double   ctx = cos(TANGX )       ; // cos of X target angle
 125    double   sty = sin(TANGY )       ; // sin of Y target angle
 126    double   cty = cos(TANGY )       ; // cos of Y target angle
 127    double   spx                     ; // sin of pixel X angle
 128    double   cpx                     ; // cos of pixel X angle
 129    double   spy                     ; // sin of pixel Y angle
 130    double   cpy                     ; // cos of pixel Y angle
 131    double   phase                   ; // phase
 132    int      xm                      ; // tickmark max
 133    int      xs                      ; // tickmark scale
 134    int      frame                   ; // animation frame count
 135    int      pct=0                   ; // percentage of complete frame
 136    int      pct0=0                  ; // previous percentage
 137    double   rfactor=0.0             ; // random amount added to position
 138    int      asize = ASIZE           ; // nominal spacing of server-sats
 139    double   offset=0.5*(1.0-asize)  ; // offsets counter integer to center
 140    double   asize3=asize*asize*asize; //
 141    double   norm = 255.0 /asize3    ; // normalizes summed amplitude
 142 #ifdef DEBUG
 143    printf( "%12.6f%12.6f%15.9f\n",  angx0, angxf, angxs );
 144    printf( "%12.6f%12.6f%15.9f\n",  angy0, angyf, angys );
 145 #endif
 146 
 147 //  set up target directory
 148 
 149    system( RMMKDIR );
 150 
 151 //  random value generation ----------------------------------------------
 152 
 153    for( iz=0; iz<asize ; iz++ ) {
 154       for( iy=0; iy<asize ; iy++ ) {
 155          for( ix=0; ix<asize ; ix++ ) {
 156             rx[ix][iy][iz]= rspace* ( drand48()-0.5 ) ;
 157             ry[ix][iy][iz]= rspace* ( drand48()-0.5 ) ;
 158             rz[ix][iy][iz]= rspace* ( drand48()-0.5 ) ;
 159    }  }  }
 160 
 161 //  outer loop, sinusodial random factor change --------------------------
 162 //  from 0 to 1 to 0
 163   
 164    for( frame=0 ; frame<NPICS ; frame++ ) {
 165       rfactor =0.5*(1.0-cos( pi2*( ((double)frame)/((double)(NPICS-1)) )));
 166 
 167 #ifndef SINGLE
 168       printf( "%04d/%04d\r", frame, NPICS ) ;
 169       fflush(stdout);
 170 #endif
 171 
 172 //  define array with even spacing ---------------------------------------
 173 
 174    for( iz=0; iz<asize ; iz++ ) {
 175       zz = space*(((double)iz)+offset );
 176       for( iy=0; iy<asize ; iy++ ) {
 177          yy = space*(((double)iy)+offset) ;
 178          zr = zz*ctilt - yy*stilt;
 179          yr = yy*ctilt + zz*stilt;
 180          xr = -2.0*yr ;
 181          for( ix=0; ix<asize ; ix++ ) {
 182             // compute array positions 
 183             xx = xr + space*(((double)ix)+offset) ;
 184             ax[ix][iy][iz] = xx + rfactor*rx[ix][iy][iz] ;
 185             ay[ix][iy][iz] = yr + rfactor*ry[ix][iy][iz] ;
 186 	    az[ix][iy][iz] = zr + rfactor*rz[ix][iy][iz] ;
 187 
 188             // compute target angle
 189 
 190             // first rotation, Y axis, by TANGX
 191             // bx = xx*ctx - zr*stx ;
 192             // by = yr ;
 193             // bz = zr*ctx + xx*stx ;
 194 
 195             // second rotation, X axis, by TANGY 
 196             // cx = bx ;
 197             // cy = by*cty - bz*sty ;
 198             // cz = bz*cty + by*sty ;
 199 
 200             // ap[ix][iy][iz] = -k*cz    ;
 201 
 202             // combine the math to find element phase
 203             ap[ix][iy][iz] = -k*( ( az[ix][iy][iz]*ctx
 204                                    +ax[ix][iy][iz]*stx)*cty 
 205                                  +  ay[ix][iy][iz]*sty )  ;
 206          }
 207    }  }
 208 #ifdef SINGLE
 209    printf( "array setup done\n");
 210 #endif
 211 
 212 // set up array and define colors -----------------------------------------
 213 
 214    im = gdImageCreateTrueColor(XSIZE, YSIZE );
 215 
 216    // Allocate standard colors
 217    black = gdImageColorAllocate(im,   0,   0,   0);
 218    white = gdImageColorAllocate(im, 255, 255, 255);
 219    sun1  = gdImageColorAllocate(im,  51,  51, 102);
 220    red   = gdImageColorAllocate(im, 255,   0,   0);
 221    green = gdImageColorAllocate(im,   0, 255,   0);
 222    blue  = gdImageColorAllocate(im,   0,   0, 255);
 223    gray  = gdImageColorAllocate(im, 128, 128, 128);
 224    dgray = gdImageColorAllocate(im,  48,  48,  48);
 225    trans = gdImageColorAllocate(im, 1, 1, 1);
 226 
 227    // allocate gray scale array 
 228    for( icolor=0 ; icolor<256 ; icolor++ ) {
 229      gcolor[icolor] = gdImageColorAllocate(im, icolor, icolor, icolor);
 230    }
 231 
 232 #ifdef SINGLE
 233    printf( "color setup done\n");
 234 #endif
 235 
 236 //  draw array with perspective ---------------------------------------------
 237 
 238    for( iz=0; iz<asize ; iz++ ) {
 239       for( iy=0; iy<asize ; iy++ ) {
 240          for( ix=0; ix<asize ; ix++ ) {
 241             zscale = 1.0 / ( 1.0- az[ix][iy][iz]/ZSCALE ) ;
 242             cir    = (int) ( zscale * SPOT );
 243             ox     = (int) ( xacent + zscale*ax[ix][iy][iz] ) ;
 244             oy     = (int) ( yacent + zscale*ay[ix][iy][iz] ) ;
 245             icolor = 180 + (int)(30.0*az[ix][iy][iz]/space ) ;
 246             if( icolor > 255 ) { icolor = 255 ; }
 247             gdImageFilledEllipse(im, ox, oy, cir, cir, gcolor[icolor] );
 248    }  }  }
 249 
 250 #ifdef SINGLE
 251    printf( "array drawn\n");
 252 #endif
 253 
 254 //   now compute phases
 255 
 256    for( ox=XLEFT ; ox<XSIZE ; ox++ ) {
 257       angx  =  angxf + angxs*((double)ox );
 258       spx   = sin( angx );
 259       cpx   = cos( angx );
 260  
 261       pct0  = pct ;
 262       pct   = (100*(ox-XLEFT))/(XSIZE-XLEFT) ;
 263 
 264 #ifndef SINGLE
 265       if( pct != pct0 ) {
 266          printf( "%04d/%04d%4d\r", frame, NPICS, pct ) ;
 267          fflush(stdout);
 268       }
 269 #endif
 270 
 271       for( oy=ybot ; oy<YSIZE ; oy++ ) {
 272          angy  =  angyf + angys*((double)oy );
 273          spy   = sin( angy );
 274          cpy   = cos( angy );
 275 
 276          // sum the components
 277          ss = 0.0 ;
 278          cc = 0.0 ;
 279 
 280 #ifdef DEBUG
 281          if( ( ox==oxc ) && ( oy==oyc ) ) {
 282               printf( "%4d%5d%12.6f%12.6f%12.6f%12.6f\n\n", ox, oy,
 283                               cpx, spx, cpy, spy );
 284          }
 285 #endif
 286 
 287          // compute z distance and phase of each source to plane
 288          for( iz=0; iz<asize ; iz++ ) {
 289             for( iy=0; iy<asize ; iy++ ) {
 290                for( ix=0; ix<asize ; ix++ ) {
 291 
 292                   // xr = ax[ix][iy][iz] ;
 293                   // yr = ax[ix][iy][iz] ;
 294                   // zr = az[ix][iy][iz] ;
 295 
 296                   // compute pixel angle
 297                   // first rotation, Y axis, by TANGX
 298                   // bx = xr*cpx - zr*spx ;
 299                   // by = yr ;
 300 	          // bz = zr*cpx + xr*spx ;
 301    
 302                   // second rotation, X axis, by TANGY 
 303    
 304                   // cx = bx ;
 305                   // cy = by*cpy - bz*spy ;
 306                   // cz = bz*cpy + by*spy ;
 307    
 308                   // ap[ix][iy][iz] = -k*cz    ;
 309                   // ap[ix][iy][iz] = -k*( (zr*cpx+xr*spx)*cpy + yr*spy);
 310    
 311                   phase = k*( ay[ix][iy][iz]*spy +
 312                           ( az[ix][iy][iz]*cpx + ax[ix][iy][iz]*spx)*cpy ) ;
 313                   ss += sin( phase + ap[ix][iy][iz] );
 314                   cc += cos( phase + ap[ix][iy][iz] );
 315 
 316 #ifdef DEBUG
 317                   if( ( ox==oxc ) && ( oy==oyc ) ) {
 318                      printf( "%3d%3d%3d%12.6f%12.6f%12.6f%12.6f\n", ix,iy,iz,
 319                               phase, ap[ix][iy][iz], ss, cc );
 320                   }
 321 #endif
 322                }
 323             }
 324          } 
 325 
 326          mag = (int) (norm*sqrt( ss*ss + cc*cc ));
 327          if( abs(mag-181) < 10 ) {
 328             gdImageSetPixel( im, ox, oy, red         ) ;
 329          }
 330          else {
 331             gdImageSetPixel( im, ox, oy, gcolor[mag] ) ;
 332          }
 333    }  }
 334 
 335 #ifdef SINGLE
 336    printf( "intensity points drawn\n");
 337 #endif
 338 
 339    // draw tickmarks on edges ----------------------------------------------
 340 
 341    // draw x tickmarks 
 342    xm = (int) (0.001+0.5*(DANGW/DANGS)) ; 
 343    xs = (XSIZE-XLEFT)*DANGS/DANGW;
 344    for( ix = -xm ; ix <= xm ; ix++ ) {
 345       ox = oxc + (int) ( xs*ix ) ;
 346       gdImageLine( im, ox, ybot, ox, ybot-5, white );
 347 
 348       sprintf( labstring, "%3d", ( (int)( ix*DANGS ) ) );
 349       gdImageStringFT( im, NULL,                   // imagespace, bounding
 350                        white, FNT, FSZ,            // color, font, fontsize
 351                        0.0, ox+FSZ/2, ybot-FSZ/2,  // angle, x, y
 352                        labstring );                // the text
 353 
 354    }
 355    ox = (int) ( xacent ) ;
 356    gdImageLine( im, xacent, ybot, ox, ybot-10, white );
 357 
 358    // draw y tickmarks
 359    xm = (int) (0.001+0.5*(angH/DANGS)) ; 
 360    for( ix = -xm ; ix <= xm ; ix++ ) {
 361       oy = oyc + (int) ( xs*ix ) ;
 362       gdImageLine( im, xleft, oy,  xleft-5, oy, white );
 363 
 364       sprintf( labstring, "%3d", ( (int)( ix*DANGS ) ) );
 365       gdImageStringFT( im, NULL,                   // imagespace, bounding
 366                        white, FNT, FSZ,            // color, font, fontsize
 367                        0.0, xleft-4*FSZ, oy+FSZ/2, // angle, x, y
 368                        labstring );                // the text
 369 
 370    }
 371    ox = (int) ( xacent ) ;
 372    gdImageLine( im, xleft, oyc, xleft-10, oyc, white );
 373 
 374    // output the frame ------------------------------------------------------
 375 
 376    sprintf( framename, PNGFMT , frame );
 377    pngout = fopen( framename, "wb");
 378    gdImagePngEx( im, pngout, 1 );
 379    gdImageDestroy(im);
 380    fclose(pngout);
 381 
 382 }  // end of single frame
 383 
 384    // link to the other half of the sequence
 385    for( i=0 ; i < NPICS ; i++ ) {
 386       sprintf( link0 , PNGFMT, i );
 387       sprintf( link1 , PNGFMT, 2*NPICS-i );
 388       link( link0, link1 ) ;
 389    }
 390 
 391    system( PNG2SWF );
 392 }