Attachment 'a3d03.c'

Download

   1 // radio array, grating lobes versus size of array
   2 //
   3 // compile with:  cc -o a3d03 a3d03.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 #include "gd.h"
  16 #include "math.h"
  17 #include <stdio.h>
  18 #include <unistd.h>
  19 
  20 #define  RMMKDIR     "rm -rf a3d03dir ; mkdir a3d03dir"
  21 #define  PNG2SWF     "png2swf -o a3d03.swf -r 2 a3d03dir/*.png" 
  22 #define  PNGFMT      "a3d03dir/a%04d.png"
  23 
  24 #define  YCENT       130
  25 #define  WAVELENGTH  2.33         // radio wavelength in pixel units
  26 #define  SPACE       6.0          // array spacing in pixel units
  27 #define  NPICS       100          // number of frames
  28 
  29 #define  TANGX       0.0          // target X angle in degrees
  30 #define  TANGY       0.0          // target Y angle in degrees
  31 #define  DANGX       0.0          // display center X angle in degrees
  32 #define  DANGY       0.0          // display center Y angle in degrees
  33 
  34 #define  DANGW       150.0        // display angle width in degrees
  35 #define  DANGS       10.0         // display angle step
  36 
  37 #define  ZSCALE      400.0        // depth scaling
  38 // #define  TILT        0.04         // array tilt around X axis, radians
  39 #define  TILT        0.00         // array tilt around X axis, radians
  40 #define  SPOT        4            // size of drawn spot
  41 #define  ASIZE       8            //
  42 #define  ASMIN       2            //
  43 #define  ASMAX       20           //
  44 #define  XLEFT       40           // left side pixel count
  45 #define  XSIZE       1000         // display window in pixels
  46 #define  YSIZE       750          // display window in pixels
  47 #define  FNT         "DejaVuMonoSans"
  48 #define  FSZ         10
  49 
  50 // ==========================================================================
  51 
  52 int main() {
  53 
  54    gdImagePtr im                ; 
  55    int      npics   = ASMAX     ;
  56    int      asize   = ASIZE     ;
  57    double   space   = SPACE     ;
  58    double   xleft   = XLEFT     ;
  59    double   yacent  = YCENT     ;
  60    double   xacent  = 0.5*(XSIZE+XLEFT) ;
  61    int      ybot    = 2* (int)(yacent) ;
  62    double   angH    = DANGW*(YSIZE-ybot)/(XSIZE-XLEFT);
  63    double   pi2     = 8.0 * atan( 1.0 );
  64    double   deg2rad = pi2/360.0 ;
  65 
  66    int      oxc = (XSIZE+XLEFT)/2 ;
  67    int      oyc = (YSIZE+ybot)/2 ;
  68 
  69    double   dangx                                ; // pixel X angle radians
  70    double   dangx0  = deg2rad*(DANGX-0.5*DANGW)  ; // left pixel value, radians
  71    double   dangx1  = deg2rad*(DANGX+0.5*DANGW)  ; // right pixel value, radians
  72    double   dangxs  = deg2rad*DANGW/(XSIZE-XLEFT); // pixel X slope radians
  73    double   dangxf  = dangx0-dangxs*XLEFT        ; // pixel X offset radians
  74    double   dangy                                ; // pixel Y angle
  75    double   dangy0  = deg2rad*(DANGY+0.5*angH)   ; // top pixel value, radians
  76    double   dangy1  = deg2rad*(DANGY-0.5*angH)   ; // bottom pixel value, radians
  77    double   dangys  = deg2rad*angH/(YSIZE-ybot)  ; // pixel Y slope radians
  78    double   dangyf  = dangy0-dangys*YSIZE        ; // pixel Y offset radians
  79    double   k       = pi2/WAVELENGTH             ; // wavenumber
  80    double   tilt    = TILT ;
  81    FILE     *pngout       ;
  82    char     dirname[80]   ;
  83    char     framename[80] ;
  84    char     labstring[80] ;
  85    char     link0[80]     ;
  86    char     link1[80]     ;
  87    int      black, sun1, white, red, green, blue, gray, dgray, trans ;
  88    int      gcolor[256];
  89    double   ax[ASMAX][ASMAX][ASMAX];            // x position of element
  90    double   ay[ASMAX][ASMAX][ASMAX];            // y position of element
  91    double   az[ASMAX][ASMAX][ASMAX];            // z position of element
  92    double   ap[ASMAX][ASMAX][ASMAX];            // phase of element
  93    double   xx, yy, zz ;
  94    double   xr, yr, zr ;
  95    double   bx, by, bz ;    // coordinates rotated on y axis
  96    double   cx, cy, cz ;    // coordinates rotated on x axis
  97    int      ix, iy, iz ;
  98    int      ox, oy     ;    // output grid
  99    double   zscale ;
 100    double   offset ;
 101    int      xdraw, ydraw ;
 102    int      icolor ;
 103    int      cir    ;
 104    double   stilt = sin( tilt );
 105    double   ctilt = cos( tilt );
 106    int      mag               ;  // pixel magnitude, 0-255 ;
 107    double   ss, cc            ;  //summed sin and cos components
 108    double   stx = sin(TANGX ) ;
 109    double   ctx = cos(TANGX ) ;
 110    double   sty = sin(TANGY ) ;
 111    double   cty = cos(TANGY ) ;
 112    double   spx               ;  // pixel angle sin and cos
 113    double   cpx               ;
 114    double   spy               ;
 115    double   cpy               ;
 116    double   asize3 = 0        ;
 117    double   asize3old=0       ;
 118    double   norm = 255.0 /asize3 ;
 119    double   phase             ;
 120    int      m                 ;  // tickmark max
 121    int      s                 ;  // tickmark scale 
 122    int      i                 ;  // general counter
 123    int      frame             ;
 124    int      pct=0             ;  // percentage of complete frame
 125    int      pct0=0            ;
 126 
 127    int      prog0=0           ;
 128    int      prog1=0           ;
 129    int      progmax=0         ;
 130    double   progfrac          ;
 131    double   pctmult           ;
 132 
 133 #ifdef DEBUG
 134    printf( "%12.6f%12.6f%12.6f%15.9f\n",  dangx0, dangx1, dangxf, dangxs );
 135    printf( "%12.6f%12.6f%12.6f%15.9f\n",  dangy0, dangy1, dangyf, dangys );
 136 #endif
 137 
 138 //  estimate progress
 139    for( i=ASMIN ; i <= ASMAX ; i++ ) { progmax += i*i*i ; }
 140 
 141 //  set up target directory
 142 
 143    system( RMMKDIR );
 144 
 145 //  outer loop, sinusodial spacing change --------------------------------
 146 //  from SPACEMIN to SPACE ( then link and fill in )
 147   
 148    for( frame=ASMIN ; frame<=ASMAX ; frame++ ) {
 149       asize = frame ;
 150       asize3old = asize3     ;
 151       asize3 = asize*asize*asize ;
 152       offset = 0.5*(1.0-asize) ;
 153       norm   = 255.0 /asize3 ;
 154       prog0  = prog1 ;
 155       prog1  = prog0 + (int) asize3 ;
 156 
 157 #ifndef SINGLE
 158       printf( "%04d/%04d\r", frame, npics ) ;
 159       fflush(stdout);
 160 #endif
 161 
 162 
 163 //  define array with even spacing ---------------------------------------
 164 
 165    for( iz=0; iz<asize ; iz++ ) {
 166       zz = space*(((double)iz)+offset );
 167       for( iy=0; iy<asize ; iy++ ) {
 168          yy = space*(((double)iy)+offset) ;
 169          zr = zz*ctilt - yy*stilt;
 170          yr = yy*ctilt + zz*stilt;
 171          xr = -2.0*yr ;
 172          for( ix=0; ix<asize ; ix++ ) {
 173             // compute array positions 
 174             xx = xr + space*(((double)ix)+offset) ;
 175             ax[ix][iy][iz] = xx ;
 176             ay[ix][iy][iz] = yr ;
 177 	    az[ix][iy][iz] = zr ;
 178 
 179             // compute target angle
 180             // first rotation, Y axis, by TANGX
 181             // bx = xx*ctx - zr*stx ;
 182             // by = yr ;
 183             // bz = zr*ctx + xx*stx ;
 184 
 185             // second rotation, X axis, by TANGY 
 186 
 187             // cx = bx ;
 188             // cy = by*cty - bz*sty ;
 189             // cz = bz*cty + by*sty ;
 190 
 191             // ap[ix][iy][iz] = -k*cz    ;
 192             ap[ix][iy][iz] = -k*( (zr*ctx+xx*stx)*cty + yr*sty )  ;
 193 
 194          }
 195 #ifdef DEBUG
 196          printf( "%3d%3d%18.3f%18.3f%18.3f%18.3f\n", iy, iz,
 197 	    ax[2][iy][iz], ay[2][iy][iz], az[2][iy][iz], ap[2][iy][iz] );
 198 #endif
 199    }  }
 200 #ifdef SINGLE
 201    printf( "array setup done\n");
 202 #endif
 203 
 204 // set up array and define colors -----------------------------------------
 205 
 206    im = gdImageCreateTrueColor(XSIZE, YSIZE );
 207 
 208    // Allocate standard colors
 209    black = gdImageColorAllocate(im,   0,   0,   0);
 210    white = gdImageColorAllocate(im, 255, 255, 255);
 211    sun1  = gdImageColorAllocate(im,  51,  51, 102);
 212    red   = gdImageColorAllocate(im, 255,   0,   0);
 213    green = gdImageColorAllocate(im,   0, 255,   0);
 214    blue  = gdImageColorAllocate(im,   0,   0, 255);
 215    gray  = gdImageColorAllocate(im, 128, 128, 128);
 216    dgray = gdImageColorAllocate(im,  48,  48,  48);
 217    trans = gdImageColorAllocate(im,   1,   1,   1);
 218 
 219    // allocate gray scale array 
 220    for( icolor=0 ; icolor<256 ; icolor++ ) {
 221      gcolor[icolor] = gdImageColorAllocate(im, icolor, icolor, icolor);
 222    }
 223 
 224 #ifdef SINGLE
 225    printf( "color setup done\n");
 226 #endif
 227 
 228 //  draw array with perspective ---------------------------------------------
 229 
 230    for( iz=0; iz<asize ; iz++ ) {
 231       for( iy=0; iy<asize ; iy++ ) {
 232          for( ix=0; ix<asize ; ix++ ) {
 233             zscale = 1.0 / ( 1.0- az[ix][iy][iz]/ZSCALE ) ;
 234             cir    = (int) ( zscale * SPOT );
 235             xdraw  = (int) ( xacent + zscale*ax[ix][iy][iz] ) ;
 236             ydraw  = (int) ( yacent + zscale*ay[ix][iy][iz] ) ;
 237             icolor = 180 + (int)(30.0*az[ix][iy][iz]/space ) ;
 238             if( icolor > 255 ) { icolor = 255 ; }
 239             if( icolor <   0 ) { icolor =   0 ; }
 240             if( ( ydraw > 0 ) && ( ydraw < ybot ) ) { 
 241               gdImageFilledEllipse(im, xdraw, ydraw, cir, cir, gcolor[icolor] );
 242             }
 243    }  }  }
 244 
 245 #ifdef SINGLE
 246    printf( "array drawn\n");
 247 #endif
 248 
 249 //   now compute phases
 250 
 251    for( ox=XLEFT ; ox<XSIZE ; ox++ ) {
 252       dangx =  dangxf + dangxs*((double)ox );
 253       spx   = sin( dangx );
 254       cpx   = cos( dangx );
 255  
 256       pct0  = pct ;
 257       pct   = (100*(ox-XLEFT))/(XSIZE-XLEFT) ;
 258 
 259 #ifndef SINGLE
 260       if( pct != pct0 ) {
 261          pctmult   = ((double)(prog1-prog0))*((double)pct);
 262          progfrac  = (100.0*((double)prog0)+pctmult)/((double)progmax) ;
 263                  
 264          printf( "%03d/%03d%6.2f%3d\r", frame, npics, progfrac, pct ) ;
 265          fflush(stdout);
 266       }
 267 #endif
 268 
 269       for( oy=ybot ; oy<YSIZE ; oy++ ) {
 270          dangy =  dangyf + dangys*((double)oy );
 271          spy   = sin( dangy );
 272          cpy   = cos( dangy );
 273 
 274          // sum the components
 275          ss = 0.0 ;
 276          cc = 0.0 ;
 277 
 278 #ifdef DEBUG
 279          if( ( ox==oxc ) && ( oy==oyc ) ) {
 280               printf( "%4d%5d%12.6f%12.6f%12.6f%12.6f\n\n", ox, oy,
 281                               cpx, spx, cpy, spy );
 282          }
 283 #endif
 284 
 285          // compute z distance and phase of each source to plane
 286          for( iz=0; iz<asize ; iz++ ) {
 287             for( iy=0; iy<asize ; iy++ ) {
 288                for( ix=0; ix<asize ; ix++ ) {
 289 
 290                   // xr = ax[ix][iy][iz] ;
 291                   // yr = ax[ix][iy][iz] ;
 292                   // zr = az[ix][iy][iz] ;
 293 
 294                   // compute pixel angle
 295                   // first rotation, Y axis, by TANGX
 296                   // bx = xr*cpx - zr*spx ;
 297                   // by = yr ;
 298 	          // bz = zr*cpx + xr*spx ;
 299    
 300                   // second rotation, X axis, by TANGY 
 301    
 302                   // cx = bx ;
 303                   // cy = by*cpy - bz*spy ;
 304                   // cz = bz*cpy + by*spy ;
 305    
 306                   // ap[ix][iy][iz] = -k*cz    ;
 307                   // ap[ix][iy][iz] = -k*( (zr*cpx+xr*spx)*cpy + yr*spy);
 308    
 309                   phase = k*( ay[ix][iy][iz]*spy +
 310                           ( az[ix][iy][iz]*cpx + ax[ix][iy][iz]*spx)*cpy ) ;
 311                   ss += sin( phase + ap[ix][iy][iz] );
 312                   cc += cos( phase + ap[ix][iy][iz] );
 313 
 314 #ifdef DEBUG
 315                   if( ( ox==oxc ) && ( oy==oyc ) ) {
 316                      printf( "%3d%3d%3d%12.6f%12.6f%12.6f%12.6f\n", ix,iy,iz,
 317                               phase, ap[ix][iy][iz], ss, cc );
 318                   }
 319 #endif
 320                }
 321             }
 322          } 
 323 
 324          mag = (int) (norm*sqrt( ss*ss + cc*cc ));
 325          if( abs(mag-181) < 10 ) {
 326             gdImageSetPixel( im, ox, oy, red         ) ;
 327          }
 328          else {
 329             gdImageSetPixel( im, ox, oy, gcolor[mag] ) ;
 330          }
 331    }  }
 332 
 333 #ifdef SINGLE
 334    printf( "intensity points drawn\n");
 335 #endif
 336 
 337    // draw tickmarks on edges ----------------------------------------------
 338 
 339    // draw x tickmarks 
 340    m = (int) (0.001+0.5*(DANGW/DANGS)) ; 
 341    s = (XSIZE-XLEFT)*DANGS/DANGW;
 342    for( i = -m ; i <= m ; i++ ) {
 343       ox = oxc + (int) ( s*i ) ;
 344       gdImageLine( im, ox, ybot, ox, ybot-5, white );
 345 
 346       sprintf( labstring, "%3d", ( (int)( i*DANGS ) ) );
 347       gdImageStringFT( im, NULL,                   // imagespace, bounding
 348                        white, FNT, FSZ,            // color, font, fontsize
 349                        0.0, ox+FSZ/2, ybot-FSZ/2,  // angle, x, y
 350                        labstring );                // the text
 351 
 352    }
 353    ox = (int) ( xacent ) ;
 354    gdImageLine( im, xacent, ybot, ox, ybot-20, white );
 355 
 356    // draw y tickmarks
 357    m = (int) (0.001+0.5*(angH/DANGS)) ; 
 358    for( i = -m ; i <= m ; i++ ) {
 359       oy = oyc + (int) ( s*i ) ;
 360       gdImageLine( im, xleft, oy,  xleft-5, oy, white );
 361 
 362       sprintf( labstring, "%3d", ( (int)( i*DANGS ) ) );
 363       gdImageStringFT( im, NULL,                   // imagespace, bounding
 364                        white, FNT, FSZ,            // color, font, fontsize
 365                        0.0, xleft-4*FSZ, oy+FSZ/2, // angle, x, y
 366                        labstring );                // the text
 367 
 368    }
 369    ox = (int) ( xacent ) ;
 370    gdImageLine( im, xleft, oyc, xleft-20, oyc, white );
 371 
 372    // output the frame ------------------------------------------------------
 373 
 374    sprintf( framename, PNGFMT , frame );
 375    pngout = fopen( framename, "wb");
 376    gdImagePngEx( im, pngout, 1 );
 377    gdImageDestroy(im);
 378    fclose(pngout);
 379 
 380 }  // end of single frame
 381 
 382    // now, hard link some frames to existing frames to "fill in the picture"
 383    for( i=0 ; i < ASMIN ; i++ ) {
 384       sprintf( link0 , PNGFMT, ASMIN );
 385       sprintf( link1 , PNGFMT, i );
 386       link( link0, link1 ) ;
 387       sprintf( link0 , PNGFMT, ASMAX );
 388       sprintf( link1 , PNGFMT, ASMAX+1+i );
 389       link( link0, link1 ) ;
 390    }
 391 
 392    // fill in the reverse
 393    for( i=0 ; i <= ASMIN+ASMAX ; i++ ) {
 394       sprintf( link0 , PNGFMT, i );
 395       sprintf( link1 , PNGFMT, (1+2*(ASMIN+ASMAX))-i );
 396       link( link0, link1 ) ;
 397    }
 398 
 399    system( PNG2SWF );
 400 }

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
  • [get | view] (2009-06-09 20:26:13, 12.2 KB) [[attachment:a3d02.c]]
  • [get | view] (2009-06-09 20:27:52, 12453.1 KB) [[attachment:a3d02.swf]]
  • [get | view] (2009-06-09 20:26:26, 13.5 KB) [[attachment:a3d03.c]]
  • [get | view] (2009-06-09 20:28:29, 3594.1 KB) [[attachment:a3d03.swf]]
  • [get | view] (2009-06-09 20:26:36, 12.1 KB) [[attachment:a3d04.c]]
  • [get | view] (2009-06-09 20:29:46, 14017.0 KB) [[attachment:a3d04.swf]]
  • [get | view] (2009-06-09 20:25:59, 15.1 KB) [[attachment:a3d05.c]]
  • [get | view] (2009-06-09 20:31:27, 20001.4 KB) [[attachment:a3d05.swf]]
  • [get | view] (2009-06-09 20:25:47, 12.1 KB) [[attachment:adish01.c]]
  • [get | view] (2009-06-09 20:25:34, 2843.1 KB) [[attachment:adish01.swf]]
  • [get | view] (2009-03-23 17:51:06, 228.2 KB) [[attachment:inter_array.png]]
 All files | Selected Files: delete move to page

You are not allowed to attach a file to this page.