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.You are not allowed to attach a file to this page.