attachment:a6.c of Brown2013

   1 // radio array, 5x5x5, sine warp spacing, power scale
   2 //
   3 // compile with:      gcc -o a6 a6.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 // libgd origin at upper left
  12 
  13 #include "gd.h"
  14 #include "math.h"
  15 #include <stdio.h>
  16 #include <stdlib.h>
  17 
  18 #define  PNGFMT      "a6.png"
  19 
  20 //#define  VSPACE      10.0         // sin variation
  21 #define  VSPACE      0.0          // sin variation
  22 
  23 #define  WAVELENGTH  10.0         // radio wavelength in pixel units
  24 #define  SPACE       40.0         // array spacing in pixel units
  25 
  26 #define  TANGX       0.46891338   // target X angle in radians
  27 #define  TANGY       0.0          // target Y angle in radians
  28 #define  DANGX       0.0          // output display center X angle in degrees
  29 #define  DANGY       0.0          // output display center Y angle in degrees
  30 #define  DANGH       70.0         // output display angle height in degrees
  31 #define  DANGW       70.0         // output display angle width  in degrees
  32 #define  DANGS       10.0         // output display angle step
  33 #define  DANGE       30.0         // output earth display angle
  34 
  35 #define  ZSCALE      400.0        // depth scaling
  36 #define  TILT        0.125663706  // array tilt around X axis, radians
  37 // #define  TILT        0.00         // array tilt around X axis, radians
  38 #define  SPOT        10           // size of drawn spot
  39 #define  ASIZE       5            // number of elements per side
  40 #define  ASIZE3      125          // number of elements
  41 
  42 #define  XSIZE       1280         // display window in pixels
  43 #define  YSIZE       720          // display window in pixels
  44 #define  OX0         600          // output array left
  45 #define  OX1         1250         // output array right
  46 #define  OY0          11          // output array top
  47 #define  OY1         661          // output array bottom
  48 #define  YACENT      360          // center of array object display
  49 #define  XACENT      300          // center of array object display
  50 #define  FNT         "DejaVuMonoSans"
  51 #define  FSZ         28
  52 
  53 // ==========================================================================
  54 
  55 int main() {
  56 
  57    double   pi2 = 8.0*atan(1.0)  ;
  58    double   d2r = pi2/360.0      ; // degrees to radians
  59    double   k   = pi2/WAVELENGTH ; // wavenumber
  60    FILE     *pngout              ; // file handle for PNG output frame
  61    char     labstring[80]        ; // used for labelling
  62    double   ax[ASIZE3]           ; // x position of element
  63    double   ay[ASIZE3]           ; // y position of element
  64    double   az[ASIZE3]           ; // z position of element
  65    double   ap[ASIZE3]           ; // phase of element
  66    double   rx[ASIZE3]           ; // randomized position of element
  67    double   ry[ASIZE3]           ; // randomized position of element
  68    double   rz[ASIZE3]           ; // randomized position of element
  69    int      icolor               ; // index for computing colors
  70    int      i                    ; // general counter
  71    int      fs                   ; // large font size
  72 
  73 // --------------------------------------------------------------------------
  74 
  75 //  dither value generation ----------------------------------------------
  76 
  77    int jcnt=0, jx, jy, jz ;
  78    double radx, rady, radz, sinx, siny, sinz, cosx, cosy, cosz ;
  79    double rads = 8.0*atan(1.0) / ((double) ASIZE );
  80    for( jz=0; jz<ASIZE ; jz++ ) {
  81       radz = rads*(0.5+(double)jz) ;
  82       cosz = cos( radz );
  83       sinz = sin( radz );
  84       for( jy=0; jy<ASIZE ; jy++ ) {
  85       rady = rads*(0.5+(double)jy) ;
  86          cosy = cos( rady );
  87          siny = sin( rady );
  88          for( jx=0; jx<ASIZE ; jx++ ) {
  89             radx = rads*(0.5+(double)jx) ;
  90             cosx = cos( radx );
  91             sinx = sin( radx );
  92             rx[jcnt] = VSPACE * ( siny - cosz );
  93             ry[jcnt] = VSPACE * ( sinz - cosx );
  94             rz[jcnt] = VSPACE * ( sinx - cosy );
  95             jcnt++ ;
  96    } } } 
  97 
  98    double ffrac = 1.0 ;
  99    double rfactor = 0.5*(1.0-cos(0.5*pi2*ffrac));
 100 
 101 //  define array with even spacing ------------------------------------------
 102 
 103    double   stx    = sin( (double) TANGX ) ;
 104    double   ctx    = cos( (double) TANGX ) ;
 105    double   sty    = sin( (double) TANGY ) ;
 106    double   cty    = cos( (double) TANGY ) ;
 107    double   tilt   = (double) TILT         ;
 108    double   stilt  = sin(          tilt  ) ;
 109    double   ctilt  = cos(          tilt  ) ;
 110    double   asize  = (double) ASIZE        ;  
 111    double   offset = 0.5*(1.0-asize)       ;
 112    double   space  = SPACE                 ;
 113    int      ix, iy, iz                     ;
 114 
 115    int      acnt = 0                       ;
 116    for( iz=0; iz<ASIZE ; iz++ ) {
 117       double zz = space*(((double)iz)+offset) ;
 118       for( iy=0; iy<ASIZE ; iy++ ) {
 119          double yy = space*(((double)iy)+offset) ;
 120          double zr = zz*ctilt - yy*stilt;
 121          double yr = yy*ctilt + zz*stilt;
 122          double xr = -2.0*yr ;         // xr offset location moves with orbit
 123          for( ix=0; ix<ASIZE ; ix++ ) {
 124             // compute array positions 
 125             double xx = xr + space*(((double)ix)+offset) ;
 126             xx += rfactor*rx[acnt] ;
 127             yr += rfactor*ry[acnt] ;
 128             zr += rfactor*rz[acnt] ;
 129 
 130             ax[acnt] = xx ;
 131             ay[acnt] = yr ;
 132             az[acnt] = zr ;
 133 
 134             // compute target angle
 135             // first rotation, Y axis, by TANGX
 136             // double bx = xx*ctx - zr*stx ;
 137             // double by = yr ;
 138             // double bz = zr*ctx + xx*stx ;
 139 
 140             // second rotation, X axis, by TANGY 
 141 
 142             // double cy = by*cty - bz*sty ;
 143             // double cz = bz*cty + by*sty ;
 144 
 145             // ap[i] = -k*cz    ;
 146             ap[acnt] = -k*( (zr*ctx+xx*stx)*cty + yr*sty )  ;
 147             acnt++ ;
 148    }  }  }
 149 
 150 
 151    // set up array and define colors -----------------------------------------
 152 
 153    gdImagePtr im = gdImageCreateTrueColor(XSIZE, YSIZE );
 154 
 155    // Allocate standard colors
 156    int black = gdImageColorAllocate(im,   0,   0,   0);
 157    int white = gdImageColorAllocate(im, 255, 255, 255);
 158    int sun1  = gdImageColorAllocate(im,  51,  51, 102);
 159    int red   = gdImageColorAllocate(im, 255,   0,   0);
 160    int green = gdImageColorAllocate(im,   0, 255,   0);
 161    int dgreen= gdImageColorAllocate(im,   0, 128,   0);
 162    int blue  = gdImageColorAllocate(im,   0,   0, 255);
 163    int gray  = gdImageColorAllocate(im, 128, 128, 128);
 164    int dgray = gdImageColorAllocate(im,  48,  48,  48);
 165    int trans = gdImageColorAllocate(im, 1, 1, 1);
 166 
 167 #define C5MAX  256 
 168 #define C2MAX  128
 169 #define C1MAX   64
 170 #define MMTHR  0.4
 171 
 172    gdImageSetThickness( im, 3 );
 173 
 174    int gcolor[C5MAX]              ; // color map numbers
 175 
 176    // allocate gray scale array 
 177    for( icolor=0 ; icolor<C5MAX ; icolor++ ) {
 178       int col = 255-icolor ;
 179       gcolor[icolor] = gdImageColorAllocate(im, col, col, col);
 180    }
 181 
 182    // White background ??
 183 
 184    gdImageFilledRectangle( im, 0, 0, XSIZE-1, YSIZE-1, white );
 185 
 186    // Slide Label -----------------------------------------------------------
 187   
 188    fs= 40 ;
 189    gdImageStringFT( im, NULL,           // imagespace, bounding box
 190       black , FNT, fs, 0.0,             // color, font, fontsize, angle
 191       20     , fs+25  ,                 // x, y
 192       "Thinsat Array" );                // text
 193 
 194    fs= 30 ;
 195    gdImageStringFT( im, NULL,           // imagespace, bounding box
 196       black , FNT, fs, 0.0,             // color, font, fontsize, angle
 197       50     , 120    ,                 // x, y
 198       "Regular Spacing" );              // text
 199 //    "Sin Dithered Spacing" );         // text
 200 
 201    gdImageStringFT( im, NULL,           // imagespace, bounding box
 202 //    dgreen , FNT, fs, 0.0,            // color, font, fontsize, angle
 203       red , FNT, fs, 0.0,            // color, font, fontsize, angle
 204       60     , 540    ,                 // x, y
 205 //    "Dither washes\nout concentrated\ngrating lobes" );         // text
 206       "Uniform array\nmakes concentrated\ngrating lobes" );       // text
 207 
 208 //  draw array with perspective ---------------------------------------------
 209 
 210    for( i=0; i<ASIZE3 ; i++ ) {
 211       double zscale = 0.8 / ( 1.0- az[i]/((double) ZSCALE ) ) ;
 212       int cir       =          (int) ( zscale*SPOT  ) ;
 213       int xdraw     = XACENT + (int) ( zscale*ax[i] ) ;
 214       int ydraw     = YACENT + (int) ( zscale*ay[i] ) ;
 215       icolor        = 128 + (int)(16.0*az[i]/space )  ;
 216       if( icolor > 255 ) { icolor = 255 ; }
 217       gdImageFilledEllipse(im, xdraw, ydraw, cir, cir, gcolor[icolor] );
 218    }  
 219 
 220 // set up scaling
 221 
 222    double mmax   = pow(asize,6.0)             ; // largest expected magnitd
 223    double mmid   = MMTHR*mmax                 ; // log scale threshold
 224    double mmin   = 1e-4*mmax                  ; // smallest plotted magnitd
 225    double cmax   = (double) C5MAX             ; // largest expected color
 226    double cmin   = 0.0                        ; // smallest expected color
 227    double am1    = (cmax-cmin)/(mmax-mmin)    ; // power slope to color
 228    double am0    = cmin - am1*mmin            ; // power offset to color
 229    double cmid   = am1*mmid + am0             ; // transition to log scale
 230    double al1    = (cmid-cmin)/log(mmid/mmin) ; // log slope to color
 231    double al0    = cmin - al1*log(mmin)       ; // log offset to color
 232    double ox0    = (double) OX0               ; // output display left
 233    double ox1    = (double) OX1               ; // output display right
 234    double oy0    = (double) OY0               ; // output display top
 235    double oy1    = (double) OY1               ; // output display bottom
 236    double dangy  = (double) DANGY             ; // display vert angle offset
 237    double dangx  = (double) DANGX             ; // display horz angle offset
 238    double dangw  = (double) DANGW             ; // display width degrees
 239    double dangh  = (double) DANGH             ; // display height degrees
 240    double dangx0 = d2r*(dangx-0.5*dangw)      ; // left pixel value, radians
 241    double dangx1 = d2r*(dangx+0.5*dangw)      ; // right pix value, radians
 242    double dangxs = d2r*dangw/(ox1-ox0)        ; // pixel X slope radians
 243    double dangxf = dangx0-dangxs*ox0          ; // pixel X offset radians
 244    double dangy0 = d2r*(dangy-0.5*dangh)      ; // top pixel value, radians
 245    double dangy1 = d2r*(dangy+0.5*dangh)      ; // bottom pix value, radians
 246    double dangys = d2r*dangh/(oy1-oy0)        ; // pixel Y slope radians
 247    double dangyf = dangy0-dangys*oy0          ; // pixel Y offset radians
 248    int    ox, oy                              ; // output grid indices
 249    int    oxc    = (OX0+OX1)/2                ; // center of output display
 250    int    oyc    = (OY0+OY1)/2                ;
 251 
 252 //   now compute phases
 253 
 254    for( ox=OX0 ; ox<OX1 ; ox++ ) {
 255       double  angx = dangxf + dangxs*((double)ox );
 256       double  spx  = sin( angx );
 257       double  cpx  = cos( angx );
 258 
 259       for( oy=OY0 ; oy<OY1 ; oy++ ) {
 260          double angy =  dangyf + dangys*((double)oy );
 261          double spy  = sin( angy );
 262          double cpy  = cos( angy );
 263 
 264          // sum the components
 265          double s2 = 0.0 ;
 266          double c2 = 0.0 ;
 267 
 268          // compute z distance and phase of each source to plane
 269          for( i=0; i<ASIZE3 ; i++ ) {
 270             // xr = ax[i] ;
 271             // yr = ax[i] ;
 272             // zr = az[i] ;
 273 
 274             // compute pixel angle  -----
 275             // first rotation, Y axis, by TANGX
 276             //   bx = xr*cpx - zr*spx ;
 277             //   by = yr ;
 278             //   bz = zr*cpx + xr*spx ;
 279             // second rotation, X axis, by TANGY 
 280             //   cx = bx ;
 281             //   cy = by*cpy - bz*spy ;
 282             //   cz = bz*cpy + by*spy ;
 283 
 284             // ap[i] = -k*cz    ;
 285             // ap[i] = -k*( (zr*cpx+xr*spx)*cpy + yr*spy);
 286    
 287             double phase = k*( ay[i]*spy + ( az[i]*cpx + ax[i]*spx)*cpy ) ;
 288             phase += ap[i]        ;
 289             s2    += sin( phase ) ;
 290             c2    += cos( phase ) ;
 291          } 
 292 
 293          // power/intensity computation
 294          double mm  = s2*s2 + c2*c2 ;
 295          if( mm < mmin )  mm = mmin ;    // nothing dimmer than mmin
 296          int mag = (int) ( am1*mm+am0 ) ;
 297          if( mm < mmid ) mag = al1*log(mm)+al0 ;
 298          if( mag > C5MAX ) mag = C5MAX  ;    // nothing brighter 
 299 
 300          gdImageSetPixel( im, ox, oy, gcolor[mag] ) ;
 301    }  }
 302 
 303    // draw earth circle ----------------------------------------------------
 304 
 305    int ex = (2*DANGE*(OX1-OX0))/DANGW ;
 306    int ey = (2*DANGE*(OY1-OY0))/DANGH ;
 307 
 308    gdImageArc( im, oxc, oyc, ex, ey, 0, 360, green );
 309  
 310    // draw box  ------------------------------------------------------------
 311    
 312    gdImageLine( im, OX0, OY0, OX1, OY0, black );
 313    gdImageLine( im, OX1, OY1, OX1, OY0, black );
 314    gdImageLine( im, OX1, OY1, OX0, OY1, black );
 315    gdImageLine( im, OX0, OY0, OX0, OY1, black );
 316   
 317    // draw ticmarks on edges ----------------------------------------------
 318    
 319    double dangs = (double) DANGS                  ; // display angle step
 320    int    mx    = (int) (0.001+0.5*(dangw/dangs)) ; // tickmark +/- count
 321    int    my    = (int) (0.001+0.5*(dangh/dangs)) ; // tickmark +/- count
 322    double sx    = (ox1-ox0)*dangs/dangw           ; // tickmark step size
 323    double sy    = (oy0-oy1)*dangs/dangh           ; // tickmark step size
 324 
 325    // draw x tickmarks 
 326 
 327    for( i = -mx ; i <= mx ; i++ ) {
 328       ox = oxc + (int) ( sx*i ) ;
 329       gdImageLine( im, ox, OY1, ox, OY1-5, black );
 330       gdImageLine( im, ox, OY0, ox, OY0+5, black );
 331 
 332       sprintf( labstring, "%3d", ( (int)( i*dangs ) ) );
 333       gdImageStringFT( im, NULL,        // imagespace, bounding
 334          black, FNT, FSZ, 0.0,          // color, font, fontsize, angle
 335          ox-1.6*FSZ, OY1+1.5*FSZ,       // angle, x, y
 336          labstring );                   // the text
 337    }
 338    gdImageLine( im, oxc, OY1, oxc, OY1-10, black );
 339    gdImageLine( im, oxc, OY0, oxc, OY0+10, black );
 340 
 341    // draw y tickmarks
 342    
 343    for( i = -mx ; i <= mx ; i++ ) {
 344       oy = oyc + (int) ( sy*i ) ;
 345       gdImageLine( im, OX1, oy,  OX1-5, oy, black );
 346       gdImageLine( im, OX0, oy,  OX0+5, oy, black );
 347 
 348       sprintf( labstring, "%3d", ( (int)( i*dangs ) ) );
 349       gdImageStringFT( im, NULL,       // imagespace, bounding
 350          black, FNT, FSZ, 0.0,         // color, font, fontsize, angle
 351          OX0-3*FSZ, oy+FSZ/2,          // x, y,
 352          labstring );                  // the text
 353 
 354    }
 355    gdImageLine( im, OX0, oyc, OX0+10, oyc, black );
 356    gdImageLine( im, OX1, oyc, OX1-10, oyc, black );
 357 
 358    // labels on output array ------------------------------------------------
 359    gdImageStringFT( im, NULL,          // imagespace, bounding box
 360       dgreen , FNT, FSZ, 0.0,           // color, font, fontsize, angle
 361       OX0+10, OY1-10,                  // x, y
 362       "Ground Pattern" );              // text
 363 
 364    // output the frame ------------------------------------------------------
 365 
 366    pngout = fopen( PNGFMT, "wb");
 367    gdImagePngEx( im, pngout, 1 );
 368    gdImageDestroy(im);
 369    fclose(pngout);
 370 }