Attachment 'gsr03.c'

Download

   1 // gcc -o gsr03 gsr03.c -lgd -lpng -lm
   2 #define  NAME   "gsr03"
   3 //
   4 // Keith Lofstrom   2013 March 4
   5 //
   6 // You are free to copy and modify this code under the GNU Public
   7 // License version 2,   http://www.gnu.org/licenses/gpl-2.0.html
   8 //
   9 // I wrote this to run under Linux.  It uses the libgd and the libpng
  10 // libraries.  I'm pretty sure you can get this working with BSD, Mac,
  11 // Windows, SunOS, or just about any other OS with a gcc compiler and
  12 // those libraries.
  13 //
  14 // 99% of the time is spent in the tight loop around line 950
  15 //
  16 // V should probably be a command line parameter.  Indeed, the
  17 // many of the defines could be read from a configuration file.
  18 // Also, this really should be modularized - the contruction of
  19 // the geodesic sphere, the movement animations, and the radio
  20 // computations could be split into three small "library" .c files.
  21 
  22 // #define  V      57         // tesselation size , 32492 elements
  23 // #define  V      16         // tesselation size ,  1692 elements
  24 #define  V      3         // tesselation size , 92 elements
  25 
  26 // NM sets the array sizes - a smarter person would use a malloc()
  27 #define  NM     600        // max tesselation size, 3,600,002 elements
  28 
  29 #define  SPACE       180   // scaling of big drawing
  30 #define  XSIZE      1024   // Y screen size (0=left, 1023=right)
  31 #define  YSIZE       768   // Y screen size (0=top,  599=bottom)
  32 #define  YSCENT     60.0   // Y center of small top images
  33 #define  YBCENT    386.0   // Y center of big bottom image
  34 #define  XBCENT    250.0   // Y center of big bottom image
  35 #define  SIZE          6   // size of spot
  36 
  37 // note - these are applied to the display, not the stored
  38 // geodesic sphere display, which is always scaled to unity
  39 
  40 #define  XEXP        1.2   // expand array in x (orbital) direction
  41 #define  YEXP        1.0   // expand array in y (NS)      direction
  42 #define  ZEXP        0.2   // expand array in z (radial)  direction
  43 
  44 // turn the array just a bit, so objects in back are visible past those
  45 // in front
  46 #define  XZROTATE   0.01
  47 
  48 #define  ASKEW       2.0   // apogee skew 
  49 
  50 #define  NPICS       240   // number of images in animation
  51 #define  RATE         10   // images per second
  52 #define  ESIZE        80   // diameter of earth disk
  53 #define  ORBSIZE      80   // radius of orbit
  54 #define  ZFACTR      0.2   // scaling for items in back
  55 #define  ASCALE       12   // earth orbiting arrow scale
  56 #define  OSCALE    0.125   // earth orbiting scale
  57 #define  CMID      128.0   // midrange color in daylight
  58 #define  CSCALE    125.0   // color scaling to distance
  59 #define  CDIM       0.25   // dimming factor for eclipse
  60 
  61 #define  FNT        "DejaVuMonoSans"
  62 #define  FSZ          40
  63 
  64 #define  SPEEDUP     (14400.0*RATE/NPICS)
  65 
  66 // radio  --------------------------------------------------------
  67 
  68 #define  WAVELENGTH  0.005 // radio wavelength 
  69 #define  TANGX       0.0   // target X angle in radians
  70 #define  TANGY       0.0   // target Y angle in radians
  71 #define  DANGX       0.0   // output display center X angle in kilometers
  72 #define  DANGY       0.0   // output display center Y angle in kilometers
  73 #define  DANGH      14.0   // output display angle height in kilometers
  74 #define  DANGW      14.0   // output display angle width  in kilometers
  75 #define  DANGS       2.0   // output display angle step in kilometers
  76 #define  DANGSCALE  7200.0 // kilometers per radian
  77 #define  OX0         540   // output array left
  78 #define  OX1        1020   // output array right ( 480 wide )
  79 #define  OY0         150   // output array top
  80 #define  OY1         630   // output array bottom ( 360 high )
  81 
  82 #define  C5MAX       320
  83 #define  C2MAX       128
  84 #define  C1MAX        64
  85 #define  MMTHR       0.4
  86 #define  MMMIN      1E-6
  87 
  88 // --------------------------------------------------------------
  89 // How toroidal orbits work.
  90 //
  91 // An elliptical orbit is faster than a similar circular orbit at
  92 // perigee, and slower at apogee.   If the delta between circular
  93 // and elliptical is ZEXP ( == e R ) then an array element halfway
  94 // between perigee and apogee in an elliptical orbit is displaced
  95 // 2 * ZEXP ahead of its circular neighbor.  We need enough ZEXP
  96 // to keep perigee and apogee objects from colliding, and also to
  97 // provide a reasonably large surface facing forwards or backwards
  98 // in orbit at the 6AM and 6PM positions of the orbit.
  99 // 
 100 // We need enough YEXP (North/South) size to make a reasonable 
 101 // large disk for radio and solar capture.   We need enough XEXP
 102 // to make radio disk, and capture the sun at the noon position.
 103 //
 104 // You can fiddle with the relative values of XEXP, YEXP, and ZEXP
 105 // to see how it shapes the array.  We probably want something a
 106 // little flattened in ZEXP to reduce relative skew. 
 107 //
 108 // This program will arbitrarily draw "right hand" toroidal orbits;
 109 // If your thumb is the orbital direction, the orbits will turn
 110 // in the direction of your fingers.  Viewing from just outside the
 111 // orbit, looking at the center of the earth, the orbit will evolve
 112 // like so:
 113 
 114 //  angle   Z (radius)             Y(NS)      X orbit
 115 //    0°    -ZEXP inside,perigee   0          0
 116 //   90°    0                      YEXP top   +2*ZEXP skew forward
 117 //  180°    +ZEXP outside,apogee   0          0
 118 //  270°    0                      YEXP bot   -2*ZEXP skew backward
 119 
 120 // --------------------------------------------------------------
 121 // How to draw a geodesic sphere - tesselated icosahedron
 122 //
 123 // make a big array of normalized array points
 124 // for a geodesic sphere radius 1 
 125 //   class 1 icosahedral base
 126 //   V up to whole bunches, possibly millions of points
 127 //   number of vertexes is 10*V^2 + 2
 128 //     
 129 //  12 vertices
 130 //  20 faces
 131 //  30 edges
 132 // 
 133 //  Procedure:
 134 //  (1) define the 12 vertices 
 135 //  (3) define the 30 edges in terms of 2 vertices each (by index)
 136 //  (2) define the 20 faces in terms of 3 vertices each (by index)
 137 //
 138 //  see http://en.wikipedia.org/wiki/Icosahedron
 139 //
 140 //            -1       +1
 141 //  axes:   x left  to  right
 142 //          y down  to  up
 143 //          z back  to  front
 144 //
 145 //  φ = (1 + √5) / 2 is the golden ratio
 146 //               x   y   z
 147 //  0, 1, 2, 3  (±1, ±φ, 0)  
 148 //  4, 5, 6, 7  (0, ±1, ±φ)   
 149 //  8, 9,10,11  (±φ, 0, ±1) 
 150 //
 151 //  The radius is sqrt( 1 + φ^2 )
 152 //  we will normalize that to one, later, and rotate so that 
 153 //  there are vertices top and bottom
 154 
 155 #include  <gd.h>
 156 #include  <math.h>
 157 #include  <stdio.h>
 158 #include  <stdlib.h>
 159 
 160 typedef struct { int a, b, c;    } Index3 ;
 161 typedef struct { int a, b;       } Index2 ;
 162 typedef struct { double x,y,z;   } Point  ;
 163 
 164 #define  SZ     ( 10*(NM+2)*NM )  // max array size (divis by 4)
 165 
 166 Point    vertex[SZ] ;             // element vertices
 167 double     xdis[SZ], ydis[SZ], zdis[SZ] ; // element location
 168 Index3   colrgb[SZ] ;             // element rgb
 169 int         col[SZ] ;             // element libdg color index
 170 int          ap[SZ] ;             // element phase
 171 
 172 #define  RMMKDIR    "rm -rf %sdir ; mkdir %sdir"
 173 #define  PNGFMT     "%sdir/a%05d.png"
 174 #define  PNG2SWF    "png2swf -o %s.swf -r %d %sdir/*.png" 
 175 
 176 // --------------------------------------------------------------------
 177 // this array is used to generate the icosahedron main vertices
 178 // 2 is φ ,  1 is 1, 0 is 0
 179 
 180 // point 0 is at the left, and 2, 4, 6, 8, and 10 are clockwise around it
 181 // point 1 is at the right, and 3, 5, 7, 9, and 11 are clockwise around it
 182 // 0 is opposite 1, 2 is opposite 3, 4 is opposite 5, etc.
 183 
 184 static Index3 vertABC[12] = {
 185                { 2, 1, 0},   //  0 left point after rotation
 186                {-2,-1, 0},   //  1 right point after rotation
 187                { 2,-1, 0},   //  2
 188                {-2, 1, 0},   //  3
 189                { 1, 0, 2},   //  4
 190                {-1, 0,-2},   //  5
 191                { 0, 2, 1},   //  6
 192                { 0,-2,-1},   //  7
 193 	       { 0, 2,-1},   //  8
 194                { 0,-2, 1},   //  9
 195                { 1, 0,-2},   // 10
 196                {-1, 0, 2},   // 11
 197 };
 198 
 199 // --------------------------------------------------------------------
 200 // this array is used to generate the 30 edges, between the numbered 
 201 // vertices in vertABC
 202 
 203 static Index2 edge[30] = {
 204                {  0,  2 },
 205                {  0,  4 },
 206                {  0,  6 },
 207                {  0,  8 },
 208                {  0, 10 },
 209                {  1,  3 },
 210                {  1,  5 },
 211                {  1,  7 },
 212                {  1,  9 },
 213                {  1, 11 },
 214                {  2,  4 },
 215                {  2,  7 },
 216                {  2,  9 },
 217                {  2, 10 },
 218                {  3,  5 },
 219                {  3,  6 },
 220                {  3,  8 },
 221                {  3, 11 },
 222                {  4,  6 },
 223                {  4,  9 },
 224                {  4, 11 },
 225                {  5,  7 },
 226                {  5,  8 },
 227                {  5, 10 },
 228                {  6,  8 },
 229                {  6, 11 },
 230                {  7,  9 },
 231                {  7, 10 },
 232                {  8, 10 },
 233 	       {  9, 11 },
 234 };
 235 
 236 // --------------------------------------------------------------------
 237 // this array is used to generate the 20 faces, between the numbered
 238 // vertices in vertABC
 239 
 240 static Index3 face[20] = {
 241                {  0,  2,  4 },
 242                {  0,  4,  6 },
 243                {  0,  6,  8 },
 244                {  0,  8, 10 },
 245                {  0, 10,  2 },
 246                {  1,  3,  5 },
 247                {  1,  5,  7 },
 248                {  1,  7,  9 },
 249                {  1,  9, 11 },
 250                {  1, 11,  3 },
 251                {  2,  7, 10 },
 252                {  5,  7, 10 },
 253                {  5,  8, 10 },
 254                {  5,  8,  3 },
 255                {  6,  8,  3 },
 256                {  6, 11,  3 },
 257 	       {  6, 11,  4 },
 258                {  9, 11,  4 },
 259                {  9,  2,  4 },
 260                {  9,  2,  7 },
 261 };
 262 
 263 // --------------------------------------------------------------------
 264 
 265 double  psi   ; 
 266 double  znorm ;
 267 
 268 int     sun1, white, red, green, blue     ; //  color indexes 
 269 int     dred, gray, trans, lgray, cyan    ;
 270 int     dgray, yellow, magenta, black     ; 
 271 int     dcyan                             ;
 272 
 273 double  xsize   = (double)   XSIZE        ; // screen width
 274 double  xcenter = (double) ( XSIZE / 4 )  ; // center of big side view
 275 int     ex      = (int) (1.5*ORBSIZE)     ; // drawn earth center x
 276 int     ey      = (int) (1.2*ORBSIZE)     ; // drawn earth center y
 277 int     re      = ESIZE / 2               ; // drawn earth radius 
 278 int     nvert                             ; // vertices in geosphere
 279 
 280 double     pi2  = 6.2831853070795864      ;
 281 gdImagePtr im                             ; // Declare the image
 282 FILE       *pngout                        ; // Declare output file
 283 
 284 // =============================================================
 285 
 286 int  normalize( Point *a ) {
 287    double len = sqrt( a->x*a->x + a->y*a->y + a->z*a->z );
 288    if( len == 0.0 ) {
 289       fprintf(stderr, "Normalize error%12.3f%12.3f%12.3f\n", a->x, a->y,a->z );
 290       return 1 ;
 291    }
 292    a->x /= len ;
 293    a->y /= len ;
 294    a->z /= len ;
 295    // printf( "N %12.3f%12.3f%12.3f%12.3f\n",a->x, a->y, a->z, len );
 296    return  0  ;
 297 } 
 298 
 299 // -------------------------------------------------------------
 300 
 301 double   val( int type ) {
 302    switch( type ) {
 303       case -2 :  return -psi ;
 304       case -1 :  return -1.0 ;
 305       case  0 :  return  0.0 ;
 306       case  1 :  return  1.0 ;
 307       case  2 :  return  psi ;
 308    }
 309    return 0.0 ;
 310 }
 311 
 312 // -------------------------------------------------------------
 313 
 314 int     gcolor[C5MAX]          ; // color map numbers
 315 int     icolor                 ; //
 316 double  mmin, mmid, mmax       ; //
 317 double  cmin, cmid, cmax       ; //
 318 double  am0, am1               ; //
 319 double  al0, al1               ; //
 320 
 321 
 322 // given magnitude, return color index
 323 int  mcolor( double mm ) {
 324    if( mm < mmin )  mm = mmin ;    // nothing dimmer than mmin
 325    int mag = (int) ( am1*mm+am0 ) ;
 326    if( mm < mmid ) mag = al1*log(mm)+al0 ;
 327    if( mag > C5MAX-1 ) mag = C5MAX-1  ;    // nothing brighter 
 328    return gcolor[ mag ] ;
 329 }
 330 
 331 
 332 // color setup
 333 void colorstart() {
 334    int    cp, cq,  ix, of ;
 335    white   =gdImageColorAllocate( im, 255, 255, 255);
 336    black   =gdImageColorAllocate( im,   0,   0,   0);
 337    gray    =gdImageColorAllocate( im, 128, 128, 128);
 338    dgray   =gdImageColorAllocate( im,  96,  96,  96);
 339    lgray   =gdImageColorAllocate( im, 192, 192, 192);
 340    cyan    =gdImageColorAllocate( im,   0, 192, 192);
 341    dcyan   =gdImageColorAllocate( im,   0,  48,  48);
 342    red     =gdImageColorAllocate( im, 255,   0,   0);
 343    dred    =gdImageColorAllocate( im, 128,   0,   0);
 344    yellow  =gdImageColorAllocate( im, 255, 255,   0);
 345    green   =gdImageColorAllocate( im,   0, 255,   0);
 346    blue    =gdImageColorAllocate( im,   0,   0, 255);
 347    magenta =gdImageColorAllocate( im, 192,   0, 192);
 348    trans   =gdImageColorAllocate( im,   1,   1,   1);
 349 
 350    // array colors ------
 351 
 352    for( ix = 0 ; ix < nvert ; ix++ ) {
 353      int ir = colrgb[ix].a ;
 354      int ig = colrgb[ix].b ;
 355      int ib = colrgb[ix].c ;
 356      col[ix]       = gdImageColorAllocate( im, ir   , ig   , ib   );
 357      col[ix+nvert] = gdImageColorAllocate( im, ir/4 , ig/4 , ib/4 );
 358    }
 359 
 360    // radio colors ------
 361 
 362    // color definitions for radio scale
 363    // The colors run from
 364    //  cmin                     cmax
 365    //  0        64      128     320
 366    //  black    red     dkgray  white
 367    //  mmin     
 368    //  1e-4                     1.0
 369 
 370    int    bgcolor,rcolor             ; // bluegreen and red
 371    int    i ;
 372    char   labstring[32];
 373 
 374    mmax   = pow(((double)nvert),2.0)   ; // largest expected magnitd
 375    mmid   = MMTHR*mmax                 ; // log scale threshold
 376    mmin   = MMMIN*mmax                 ; // smallest plotted magnitd
 377    cmax   = (double) C5MAX             ; // largest expected color
 378    cmin   = 0.0                        ; // smallest expected color
 379    am1    = (cmax-cmin)/(mmax-mmin)    ; // power slope to color
 380    am0    = cmin - am1*mmin            ; // power offset to color
 381    cmid   = am1*mmid + am0             ; // transition to log scale
 382    al1    = (cmid-cmin)/log(mmid/mmin) ; // log slope to color
 383    al0    = cmin - al1*log(mmin)       ; // log offset to color
 384 
 385    // allocate gray scale array 
 386    for( icolor=0 ; icolor<C1MAX ; icolor++ ) {
 387      rcolor  = icolor ;
 388      bgcolor = 0 ;
 389      gcolor[icolor] = gdImageColorAllocate(im, rcolor, bgcolor, bgcolor);
 390    }
 391    for(          ; icolor<C2MAX ; icolor++ ) {
 392      rcolor  = C1MAX  ;
 393      bgcolor = icolor-C1MAX ;
 394      gcolor[icolor] = gdImageColorAllocate(im, rcolor, bgcolor, bgcolor);
 395    }
 396    for(          ; icolor<C5MAX ; icolor++ ) {
 397      rcolor  = icolor-C1MAX  ;
 398      bgcolor = icolor-C1MAX  ;
 399      gcolor[icolor] = gdImageColorAllocate(im, rcolor, bgcolor, bgcolor);
 400    }
 401 
 402    // draw color scale 
 403 
 404    for( i = 0 ; i < 6 ; i += 1 ) {
 405       double mm=mmax*pow(10.0,-i);
 406 
 407       int dxt = 500 + 88*i ;
 408       int dx0 = dxt + 47 ;
 409       int dx1 = dx0 + 1  ;
 410       int dx2 = dx1 + 20 ;
 411       int dx3 = dx2 + 1  ;
 412 
 413       int dy3 = OY1+65   ;
 414       int dy2 = dy3 - 1  ;
 415       int dy1 = dy2 - 20 ;
 416       int dy0 = dy1 - 1  ;
 417       int dyt = dy3 - 6  ;
 418 
 419       gdImageRectangle(im,       dx0, dy0, dx3, dy3, gray        );
 420       gdImageFilledRectangle(im, dx1, dy1, dx2, dy2, mcolor( mm ));
 421       
 422       sprintf( labstring, "1e-%1d", i );
 423       gdImageStringFT( im, NULL,         // imagespace, bounding
 424           white, FNT, 14, 0.0,           // color, font, fontsize, angle
 425           dxt, dyt,                      // x, y
 426           labstring );                   // the text
 427    }
 428 }
 429    
 430 
 431 // draw arrowhead (either degrees or radians ) -----------------------------
 432 
 433 void arrowhead( double xp, double yp, double asize,
 434                 double deg, double rad, int color ) {
 435    double angl = rad+pi2*deg/360.0 ;
 436    gdPoint arrow[3] ;
 437 
 438    arrow[0].x = (int)(xp);
 439    arrow[0].y = (int)(yp);
 440    arrow[1].x = (int)(xp+asize*(    -cos(angl)-0.5*sin(angl)));
 441    arrow[1].y = (int)(yp+asize*( 0.5*cos(angl)-    sin(angl)));
 442    arrow[2].x = (int)(xp+asize*(    -cos(angl)+0.5*sin(angl)));
 443    arrow[2].y = (int)(yp+asize*(-0.5*cos(angl)-    sin(angl)));
 444    gdImageFilledPolygon( im, arrow, 3, color );
 445 }
 446 
 447 // draw arrow  -------------------------------------------------------------
 448 
 449 void arrow( double x0, double y0, double x1, double y1, int color ) {
 450    double  asize = 0.02*(x1-x0) ;
 451 
 452    gdImageLine( im, (int)x0, (int)y0, (int)x1, (int)y1, color ) ;
 453    arrowhead( x1, y1, 0.02*(x1-x0), 0.0, 0.0, color ) ;
 454 }
 455 
 456 // draw rotating earth and orbiting body and side arrow  -------------
 457 
 458 void earth( double ang ) {
 459    double  x,  y  ;
 460    double  xs, ys ;
 461    int  x0, y0 ;
 462 	   int  x1, y1 ;
 463    double  an0 ;
 464    int ie    ;
 465    int orbit = 2*ORBSIZE ;
 466 
 467    gdImageFilledArc(im, ex, ey, ESIZE, ESIZE,   0, 360, dcyan, gdArc);
 468    gdImageFilledArc(im, ex, ey, ESIZE, ESIZE,  90, 270,  cyan, gdArc);
 469    gdImageArc(      im, ex, ey, orbit, orbit,  30, 330, white );
 470    gdImageArc(      im, ex, ey, orbit, orbit, 330,  30, dgray );
 471 
 472    // draw longitude lines  ( 6 = 180/30 degrees )
 473    for( ie = 0 ; ie < 6 ; ie ++ ) {
 474       an0 = (ang/6.0) + pi2*( (double)ie) / 12.0 ;
 475       x0  = ex - (int)( ((double)re) * cos(an0) );
 476       x1  = ex + (int)( ((double)re) * cos(an0) );
 477       y0  = ey + (int)( ((double)re) * sin(an0) );
 478       y1  = ey - (int)( ((double)re) * sin(an0) );
 479       gdImageLine( im, x0, y0, x1, y1, black );
 480    }
 481 
 482    // draw latitude lines
 483    for( ie = 1 ; ie < 3 ; ie++ ) {
 484       an0 = pi2*( (double) ie ) / 12.0 ;
 485       x0  = (int)( 2.0 * ((double)re) * cos(an0) );
 486       gdImageArc( im, ex, ey, x0, x0, 0, 360, black );
 487    }
 488    
 489    // draw arrowhead to show side point of view
 490    x= -( ORBSIZE + ASCALE + 2.0 ) ;
 491    y=  0.0 ;
 492    xs = x*cos(ang)+y*sin(ang) + (double) ex ;
 493    ys = y*cos(ang)-x*sin(ang) + (double) ey ;
 494  
 495    arrowhead( xs, ys, ASCALE,  0.0, -ang, white );
 496 }
 497 
 498 // ==========================================================================
 499 
 500 int main() {
 501    double  angle ;
 502    int     frame ;
 503    char    pngfilename[64];
 504    char    command[128];
 505 
 506    double  test   ;   // depth testing
 507    int     ox, oy ;
 508    int     s0     ;
 509 
 510    double  xu, yu ;
 511    int     xi, yi ;
 512    int     icol   ;
 513    double   d2r = pi2/360.0      ; // degrees to radians
 514    int     i      ;
 515    double  dv = (double) V ;
 516    char    labstring[40] ;
 517    
 518    // double   k   = pi2/WAVELENGTH ; // wavenumber
 519    // KLUDGE  approximate scaling for 1 meter thinsat spacing
 520 
 521    double  k = 8.0 * V / ( SPACE * WAVELENGTH );
 522 
 523    psi   = 0.5 * ( 1.0 + sqrt(5.0) ) ;
 524    znorm = sqrt( psi*psi + 1.0 ) ;
 525 
 526    sprintf( command, RMMKDIR, NAME, NAME );
 527    system( command );
 528 
 529    // COMPUTE THE ICOSAHEDRON -------------------------------------------
 530    // we will generate the icosahedron, then interpolate the edge points
 531    // and face points between them.   These will then be normalized to
 532    // radius 1
 533 
 534    // compute 12 vertices ------------------------- 
 535 
 536    for( nvert = 0 ; nvert < 12 ; nvert++ ) {
 537       double xtemp = val( vertABC[nvert].a );
 538       double ytemp = val( vertABC[nvert].b );
 539       double ztemp = val( vertABC[nvert].c );
 540       
 541       // rotate so x value of 0th vertex is max
 542       double xr = psi*xtemp + ytemp ;
 543       double yr = psi*ytemp - xtemp ;
 544       double zr = znorm*ztemp       ;
 545       double sxz = XZROTATE         ;
 546       double cxz = sqrt( 1.0-sxz*sxz );
 547 
 548       // rotate x to z (radial) directions, y (NS) remains same
 549       vertex[nvert].x = cxz * xr - sxz * zr ;
 550       vertex[nvert].y = yr                  ;
 551       vertex[nvert].z = cxz * zr + sxz * xr ;
 552    }
 553    // nvert is now 12
 554 
 555    // compute 30 edges ---------------------------
 556    //  using two points on the icosahedron, in the edge table
 557 
 558    for( i = 0 ; i < 30 ; i++ ) {
 559       Point va = vertex[ edge[i].a ] ;   // find icosahedron points
 560       Point vb = vertex[ edge[i].b ] ;   // find icosahedron points
 561       double da ;
 562 
 563       for( da = 1.0 ; da < dv ; da++ ) {
 564          double db = dv-da ;
 565          vertex[ nvert ].x = da*va.x + db*vb.x ;
 566          vertex[ nvert ].y = da*va.y + db*vb.y ;
 567          vertex[ nvert ].z = da*va.z + db*vb.z ;
 568          nvert++ ;   // add more vertices
 569    }  }
 570 
 571    // number of new edge points = 30*(V-1)
 572    // nvert is now 30*V - 18
 573 
 574    // compute 20 faces ---------------------------
 575    //  using three points on the icosahedron, in the face table
 576 
 577    for( i = 0 ; i < 20 ; i++ ) {
 578       Point va = vertex[ face[i].a ] ;   // find icosahedron points
 579       Point vb = vertex[ face[i].b ] ;   // find icosahedron points
 580       Point vc = vertex[ face[i].c ] ;   // find icosahedron points
 581       double da, db ;
 582 
 583       for( da = 1.0 ; da < dv ; da++ ) {
 584          for( db = 1.0 ; db < ( dv-da) ; db++ ) {
 585             double dc = dv-(da+db);
 586             vertex[ nvert ].x = da*va.x + db*vb.x + dc*vc.x ;
 587             vertex[ nvert ].y = da*va.y + db*vb.y + dc*vc.y ;
 588             vertex[ nvert ].z = da*va.z + db*vb.z + dc*vc.z ;
 589             nvert++ ;   // add more vertices
 590    }  }  }
 591 
 592    // number of new face points = 10*V^2 - 30*V + 20
 593    // nvert is now 10*V^2 + 2
 594 
 595    // color -----------------------------------------------
 596    // first we normalize each vertex point to radius 1.  Then
 597    // we will choose vertex color based on position in the base array
 598 
 599    for( i = 0 ; i < nvert ; i++ ) {
 600       normalize( &vertex[i] );
 601       colrgb[i].a = (int) ( CMID + CSCALE*vertex[i].x ) ;
 602       colrgb[i].b = (int) ( CMID + CSCALE*vertex[i].y ) ;
 603       colrgb[i].c = (int) ( CMID + CSCALE*vertex[i].z ) ;
 604    }
 605 
 606    // ICOSAHEDRON COMPUTED --------------------------------------------
 607 
 608    // graphics loop ---------------------------------------------------
 609 
 610    for( frame = 0 ; frame < NPICS ; frame++ ) {
 611       im = gdImageCreateTrueColor(XSIZE, YSIZE );
 612       sprintf( command,  PNGFMT, NAME, frame );
 613       pngout = fopen( command, "wb");
 614 
 615       angle = frame * pi2 / NPICS ;
 616       colorstart() ;
 617       earth(angle) ;
 618 
 619       //  labels -----------------------------------------------------------
 620 
 621       gdImageStringFT( im, NULL,                       // imagespace, bounding 
 622                        white, FNT, FSZ,                // color, font, fontsize
 623                        0.0, 0.00*XSIZE, 0.95*YSIZE,    // angle, x, y
 624                        " Side View" );                 // the text
 625 
 626       sprintf( command, "%7d thinsats", nvert );
 627 
 628       gdImageStringFT( im, NULL,                       // imagespace, bounding 
 629                        white, FNT, FSZ/2,              // color, font, fontsize
 630                        0.0, 0.00*XSIZE, 0.85*YSIZE,    // angle, x, y
 631                        command      );                 // the text
 632 
 633 
 634       gdImageStringFT( im, NULL,                       // imagespace, bounding
 635                        white, FNT, FSZ/2,              // color, font, fontsize
 636                        0.0, 0.26*XSIZE, 2.2*YSCENT,    // angle, x, y
 637                        "    Top View" );               // the text
 638 
 639       gdImageStringFT( im, NULL,                       // imagespace, bounding
 640                        white, FNT, FSZ/2,              // color, font, fontsize
 641                        0.0, 0.50*XSIZE, 2.2*YSCENT,    // angle, x, y
 642                        " Front View" );                // the text
 643 
 644       gdImageStringFT( im, NULL,                       // imagespace, bounding
 645                        white, FNT, FSZ/2,              // color, font, fontsize
 646                        0.0, 0.71*XSIZE, 2.2*YSCENT,    // angle, x, y
 647                        "   Side View" );               // the text
 648 
 649       //  time code -------------------------------------------------------
 650  
 651       double time0 = ( 4.0 * (double) frame ) / ( (double) NPICS ) ;
 652       char   timestring[64] ;
 653 
 654       sprintf( timestring, "%3.1f/4.0 hours,%6.0fx speedup", time0, SPEEDUP );
 655 
 656       gdImageStringFT( im, NULL,                       // imagespace, bounding 
 657                        white, FNT, 0.6*FSZ,            // color, font, fontsize
 658                        0.0, XSIZE-16.0*FSZ,            // angle, x
 659                        0.95*YSIZE, timestring );       // y, the text 
 660 
 661       //  direction arrows ------------------------------------------------
 662           // large side view, small top view, small side view
 663 
 664       arrow( 0.00*XSIZE, YBCENT, 2.0*XBCENT, YBCENT, gray ); // large side view
 665       arrow( 0.25*XSIZE, YSCENT, 0.52*XSIZE, YSCENT, gray );
 666       arrow( 0.71*XSIZE, YSCENT, 0.95*XSIZE, YSCENT, gray );
 667 
 668       //  compute element locations in X, Y, Z --------------------------
 669       //  this is where we skew and distort the array.    The Z scaling
 670       //  becomes the x distance skew.  The orbit of each element can be
 671       //  elliptical in Y (north-south)  versus Z (radial), while the X
 672       //  distance (along the orbit) gets the skew 
 673 
 674       double zm = SPACE * ZEXP ;   // radial
 675       double ym = SPACE * YEXP ;
 676       double xm = SPACE * XEXP ;
 677       double xt = SPACE * ( 1.0 + ASKEW * ZEXP ) ; // range of x values
 678       double as = ASKEW * zm/ym ;
 679       double dtest = (257.0/1024.0) ;
 680       double test0 = (-1025.0/1024.0);
 681       double test1 = ( 1025.0/1024.0);
 682       double zfactr = ZFACTR/zm ;
 683 
 684       // offsets from small top view to screen
 685       int  txo = (int) ( (0.26+0.125) * xsize );
 686       int  tyo = (int) ( YSCENT );
 687 
 688       // offsets from small front view to screen 
 689       int  fxo = (int) ( (0.48+0.125) * xsize );
 690       int  fyo = (int) ( YSCENT );
 691 
 692       // offsets from small side view to screen
 693       int ssxo = (int) ( (0.71+0.125) * xsize );
 694       int ssyo = (int) ( YSCENT );
 695 
 696       // offsets from large side view to screen
 697       int lsxo =         XBCENT ;
 698       int lsyo = (int) ( YBCENT );
 699 
 700       // compute the points of the rotated, scaled, and skewed array
 701       for( i=0 ; i < nvert ; i++ ) {
 702          zdis[i] = zm*(cos(angle)*vertex[i].z+sin(angle)*vertex[i].y ) ;
 703          ydis[i] = ym*(cos(angle)*vertex[i].y-sin(angle)*vertex[i].z ) ;
 704          xdis[i] = xm*vertex[i].x            + as*ydis[i] ;
 705       }
 706 
 707 
 708       //  layer from back, so that bigger points in front overlap smaller
 709       //  points in back.   This will make a large loop
 710       for( test = test0; test < test1 ; test += dtest ) {
 711 
 712          // compute display slices
 713          double z0 =  test*zm ;
 714          double z1 = dtest*zm + z0 ;
 715          double y0 =  test*ym ;
 716          double y1 = dtest*ym + y0 ;
 717          double x0 =  test*xt ;
 718          double x1 = dtest*xt + x0 ;
 719  
 720          // step through all the vertices miltiple times
 721          for( i=0 ; i < nvert; i++ ) {
 722 
 723             //  draw view from top (Z+) --------------------------------------
 724 
 725             if( ( ydis[i] >= y0 ) && ( ydis[i] < y1 ) ) {
 726 
 727                // small top view
 728 
 729 	       ox = txo + (int) ( 0.25*xdis[i] );
 730                oy = tyo + (int) ( 0.25*zdis[i] );
 731                gdImageSetPixel( im, ox, oy, col[i] );
 732                // gdImageFilledArc(im, ox, oy, s0, s0, 0, 360, col[i], gdArc);
 733 
 734                // orbiting tiny top view
 735 
 736                xu = -ORBSIZE - OSCALE * zdis[i] ;
 737                yu =            OSCALE * xdis[i] ; 
 738 
 739                xi = ex + (int)( xu*cos(angle)+yu*sin(angle)) ;
 740                yi = ey + (int)( yu*cos(angle)-xu*sin(angle)) ;
 741 
 742                if(  ( xi <   ex      ) 
 743                   ||( yi < ( ey-re ) )
 744                   ||( yi > ( ey+re ) ) ) icol = col[      i] ;
 745                else                      icol = col[nvert+i] ;
 746 
 747                gdImageSetPixel( im, xi, yi, icol );
 748                // gdImageFilledArc( im, xi, yi, 3, 3, 0, 360, icol, gdArc );
 749             }
 750 
 751             // side views ---------------------------------------------
 752 
 753             if( ( zdis[i] >= z0 ) && ( zdis[i] < z1 ) ) {
 754 
 755                // BIG side view 
 756                ox = lsxo + (int) (  xdis[i] );
 757                oy = lsyo + (int) ( -ydis[i] );
 758                // gdImageSetPixel( im, ox, oy, col[i] );
 759 
 760                s0 = (int) (     SIZE * (1.00 + zfactr * zdis[i] ) ) ;
 761                gdImageFilledArc(im, ox, oy, s0, s0, 0, 360, col[i], gdArc);
 762 
 763                // small side view  
 764                ox = ssxo + (int) (  0.25*xdis[i] );
 765                oy = ssyo + (int) ( -0.25*ydis[i] );
 766                gdImageSetPixel( im, ox, oy, col[i] );
 767 
 768                // s0 = (int) ( 0.35*SIZE * (1.00 + zfactr * zdis[i] ) ) ;
 769                // gdImageFilledArc(im, ox, oy, s0, s0, 0, 360, col[i], gdArc);
 770             }
 771 
 772             // small front view - apogee to left, closer ---------------
 773 
 774             if( ( xdis[i] >= x0 ) && ( xdis[i] < x1 ) ) {
 775                ox = fxo + (int) (  0.25*zdis[i] );
 776                oy = fyo + (int) ( -0.25*ydis[i] );
 777                gdImageSetPixel( im, ox, oy, col[i] );
 778                // s0 = (int) ( 0.35*SIZE * (1.30 + 2.0* zfactr * zdis[i] ) ) ;
 779                // gdImageFilledArc(im, ox, oy, s0, s0, 0, 360, col[i], gdArc);
 780             }
 781          }
 782       }
 783 
 784       // ---------------------------------------------------------------------
 785       // RADIO COMPUTATION LOOP ----------------------------------------------
 786       // setup
 787       // 1) construct the target vector from the center of the array,
 788       //    length 2pi/wavelength
 789       // 2) For each element, construct the vector from the center of
 790       // the array, and compute the dot product of that vector with 
 791       // the target vector to compute the conjugate phase of the element
 792       //
 793       // For each plotted angular point:
 794       //
 795       // 3) construct the target vector from the center of the array
 796       // to the plotted angular point, length 2pi/wavelength
 797       // 4) Repeat 2 for each element, adding to the phase difference
 798       // 5) compute the sin and cos of the phase difference, add to
 799       //    the running total for the angular point
 800       // 6) compute the magnitude, choose color from palette,
 801       //    plot point
 802       // 7) Draw box and key
 803 
 804       double ox0    = (double) OX0               ; // output display left
 805       double ox1    = (double) OX1               ; // output display right
 806       double oy0    = (double) OY0               ; // output display top
 807       double oy1    = (double) OY1               ; // output display bottom
 808       double dangy  = (double) DANGY             ; // display vert angle offset
 809       double dangx  = (double) DANGX             ; // display horz angle offset
 810       double dangw  = (double) DANGW             ; // display width kilometers
 811       double dangh  = (double) DANGH             ; // display height kilometers
 812       double dangsc = (double) DANGSCALE         ; // kilometer distance
 813       double k2r    = 1.0/dangsc                 ; // radians per kilometer
 814       double dangx0 = k2r*(dangx-0.5*dangw)      ; // left pixel value, radians
 815       double dangx1 = k2r*(dangx+0.5*dangw)      ; // right pix value, radians
 816       double dangxs = k2r*dangw/(ox1-ox0)        ; // pixel X slope radians
 817       double dangxf = dangx0-dangxs*ox0          ; // pixel X offset radians
 818       double dangy0 = k2r*(dangy-0.5*dangh)      ; // top pixel value, radians
 819       double dangy1 = k2r*(dangy+0.5*dangh)      ; // bottom pix value, radians
 820       double dangys = k2r*dangh/(oy1-oy0)        ; // pixel Y slope radians
 821       double dangyf = dangy0-dangys*oy0          ; // pixel Y offset radians
 822       int    oxc    = (OX0+OX1)/2                ; // center of output display
 823       int    oyc    = (OY0+OY1)/2                ;
 824 
 825       // make wavevector 
 826       double kx = -k * sin( d2r*TANGX )                  ;
 827       double ky = -k * cos( d2r*TANGX ) * sin( d2r*TANGY );
 828       double kz = -k * cos( d2r*TANGX ) * cos( d2r*TANGY );
 829 
 830       // conjugate phase for each element
 831 
 832       for( i=0 ; i < nvert; i++ ) {
 833          ap[i] = kx*xdis[i] + ky*ydis[i] + kz*zdis[i] ;
 834       }
 835 
 836       // --------------------------------------------------------------------
 837       // Big Computation Loop -----------------------------------------------
 838       // a bunch of points, a bunch of emitters
 839       // The outer ox and oy loops are relatively small, though I can imagine
 840       // image sizes of 20K by 20K pixels.  
 841       // The inner nvert loop is thousands to millions,
 842       // and parallelizing that would save much time.
 843 
 844       int pct  = -1 ;
 845       int pct0 =  0 ;
 846 
 847       for( ox=OX0 ; ox<OX1 ; ox++ ) {
 848          double  angx = dangxf + dangxs*((double)ox );
 849          double  spx  = sin( angx );
 850          double  cpx  = cos( angx );
 851    
 852          pct0  = pct ;
 853          pct   = (100*(ox-OX0))/(OX1-OX0) ;
 854    
 855          if( pct != pct0 ) {
 856             printf( "%04d/%04d%4d\r", frame, NPICS, pct ) ;
 857             fflush(stdout);
 858          }
 859    
 860          for( oy=OY0 ; oy<OY1 ; oy++ ) {
 861             double angy =  dangyf + dangys*((double)oy );
 862             double spy  = sin( angy );
 863             double cpy  = cos( angy );
 864    
 865             // sum the components
 866             double s2 = 0.0 ;
 867             double c2 = 0.0 ;
 868 
 869             kx = k * spx ;
 870             ky = k * cpx * spy ;
 871             kz = k * cpx * cpy ;
 872             
 873             // CENTRAL COMPUTING LOOP
 874             // element loop - spends huge time in this loop, nvert is 
 875             // typically thousands to millions
 876             for( i=0 ; i < nvert; i++ ) {
 877                double phase = ap[i] + kx*xdis[i] + ky*ydis[i] + kz*zdis[i] ;
 878                s2    += sin( phase ) ;
 879                c2    += cos( phase ) ;
 880             }
 881             // END CENTRAL COMPUTING LOOP
 882            
 883             // power/intensity computation
 884             double mm  = s2*s2 + c2*c2 ;
 885             gdImageSetPixel( im, ox, oy, mcolor(mm) ) ;
 886       }  } // END big computation triple loop
 887 
 888       // draw box  ----------------------------------------------------------
 889       
 890       gdImageLine( im, OX0, OY0, OX1, OY0, white );
 891       gdImageLine( im, OX1, OY1, OX1, OY0, white );
 892       gdImageLine( im, OX1, OY1, OX0, OY1, white );
 893       gdImageLine( im, OX0, OY0, OX0, OY1, white );
 894      
 895       // draw ticmarks on edges ----------------------------------------------
 896       
 897       double dangs = (double) DANGS                  ; // display angle step
 898       int    mx    = (int) (0.001+0.5*(dangw/dangs)) ; // tickmark +/- count
 899       int    my    = (int) (0.001+0.5*(dangh/dangs)) ; // tickmark +/- count
 900       double sx    = (ox1-ox0)*dangs/dangw           ; // tickmark step size
 901       double sy    = (oy0-oy1)*dangs/dangh           ; // tickmark step size
 902    
 903       // draw x tickmarks 
 904 
 905       int  ft     = 16 ;
 906    
 907       for( i = -mx ; i <= mx ; i++ ) {
 908          ox = oxc + (int) ( sx*i ) ;
 909          gdImageLine( im, ox, OY1, ox, OY1-5, white );
 910          gdImageLine( im, ox, OY0, ox, OY0+5, white );
 911    
 912          sprintf( labstring, "%3d", ( (int)( i*dangs ) ) );
 913          gdImageStringFT( im, NULL,        // imagespace, bounding
 914             white, FNT, ft, 0.0,           // color, font, fontsize, angle
 915             ox-1.6*ft, OY1+1.5*ft,         // angle, x, y
 916             labstring );                   // the text
 917       }
 918       gdImageLine( im, oxc, OY1, oxc, OY1-10, white );  // big tics
 919       gdImageLine( im, oxc, OY0, oxc, OY0+10, white );
 920 
 921             // draw y tickmarks
 922       
 923       for( i = -mx ; i <= mx ; i++ ) {
 924          oy = oyc + (int) ( sy*i ) ;
 925          gdImageLine( im, OX1, oy,  OX1-5, oy, white );
 926          gdImageLine( im, OX0, oy,  OX0+5, oy, white );
 927    
 928          sprintf( labstring, "%3d", ( (int)( i*dangs ) ) );
 929          gdImageStringFT( im, NULL,       // imagespace, bounding
 930             white, FNT, ft, 0.0,          // color, font, fontsize, angle
 931             OX0-3*ft, oy+ft/2,            // x, y,
 932             labstring );                  // the text
 933    
 934       }
 935       gdImageLine( im, OX0, oyc, OX0+10, oyc, white );
 936       gdImageLine( im, OX1, oyc, OX1-10, oyc, white );
 937 
 938       gdImageStringFT( im, NULL,         // imagespace, bounding
 939           white, FNT, ft+4, 0.0,         // color, font, fontsize, angle
 940           OX0+1*ft, OY1-1*ft,            // x, y
 941           "Ground spot, km");            // the text
 942    
 943       // Output the frame in PNG format.----------------------------
 944 
 945       gdImagePngEx( im, pngout, 1 );
 946       gdImageDestroy(im);
 947       fclose(pngout);
 948    }
 949 
 950    sprintf( command, PNG2SWF, NAME, RATE, NAME );
 951    system( command );
 952    return 0 ;
 953 }

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] (2013-02-27 21:27:00, 405.6 KB) [[attachment:H0_toroid.png]]
  • [get | view] (2013-03-04 17:15:39, 11.1 KB) [[attachment:ap02.c]]
  • [get | view] (2013-02-27 21:26:39, 4323.3 KB) [[attachment:ap02.swf]]
  • [get | view] (2013-03-14 00:03:25, 22.1 KB) [[attachment:gsr02]]
  • [get | view] (2013-03-05 06:53:04, 35.2 KB) [[attachment:gsr02.c]]
  • [get | view] (2013-03-04 16:48:33, 52232.0 KB) [[attachment:gsr02.swf]]
  • [get | view] (2013-03-14 00:03:41, 22.1 KB) [[attachment:gsr03]]
  • [get | view] (2013-03-05 06:53:14, 35.2 KB) [[attachment:gsr03.c]]
  • [get | view] (2013-03-04 16:50:03, 24452.2 KB) [[attachment:gsr03.swf]]
  • [get | view] (2013-02-27 16:30:39, 7.0 KB) [[attachment:incline1.png]]
 All files | Selected Files: delete move to page

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