Attachment 'ar06.c'
Download 1 // radar energy pattern
2 //
3 // compile with: gcc -o ar06 ar05.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 #define IMAGE 1
14 // #undef IMAGE
15
16 #include "gd.h"
17 #include "math.h"
18 #include "string.h"
19 #include <stdio.h>
20 #include <stdlib.h>
21
22 #define PNGFMT "ar06.png"
23
24 /*
25 #define SCANX 1000.0 // scan width in meters
26 #define LABX 200.0 // label spacing
27 #define EXAG 5.0 // pixel exaggeration
28 #define UNITS "meters"
29 #define USCALE 1.0
30 #define TUNITS 1000.0
31 #define TSCALE "ms"
32 #define OX0 700 // first output array left
33 #define OX1 2000 // first output array right
34
35 #define SCANX 0.1 // scan width in meters
36 #define LABX 0.02 // label spacing
37 #define EXAG 1.0 // pixel exaggeration
38 #define UNITS "centimeters"
39 #define USCALE 100.0
40 #define TUNITS 1000000.0
41 #define TSCALE "μs"
42 #define OX0 700 // first output array left
43 #define OX1 2000 // first output array right
44 */
45
46 // SCANX LABX EXAG USCALE TUNITS OX0 OX1 UNITS TSCALE
47 #define D0 "1000.0 200.0 5.0 1.0 1000.0 1200 3000 meters ms"
48 #define D1 "0.1 0.02 1.0 100.0 1000000.0 3200 5000 centimeters μs"
49
50 #define XSIZE 5120 // display window in pixels
51 #define OY0 100 // output array top
52 #define OY1 1200 // output array bottom
53 #define GGAIN 250
54 #define GY (OY1+130+GGAIN) // plot centerline
55 #define GT (GY+GGAIN) // time text bottom
56 #define YSIZE (GT+20)
57
58 #define EYC 1100 // earth orbit display vertical middle
59 #define EX0 70 // earth orbit display left
60 #define EX1 940 // earth orbit display right
61
62 #define LNAR 3 // line pixel width
63 #define LWID 5 // line pixel width
64
65 #define OSLOPE -0.392 // orbit slope through field
66
67 #define EINC -50.0 // inclination of target orbit
68 #define EAOP 35.0 // argument of perigee target orbit
69
70 // ------------------------------------------------------------
71
72 #define VORB 3000.0 // orbit velocity
73 #define LAMBDA 0.005 // radio wavelength in meters
74 #define TARGETA 400.0 // altitude of target
75 #define TARGETN 45 // target latitude north
76
77 #define S_WIDTH 200.0 // array width in meters
78 #define S_NUM 181 // thinsats in a row, odd
79
80 #define A_SPACE 12 // array spacing in degrees
81 #define A_NUM 7 // array spacing in degrees
82
83 #define NSAT (S_NUM*A_NUM) // total number of thinsats
84
85 #define RE 6378.0 // earth radius kilometers
86 #define R288 12789.0 // server sky radius kilometers
87
88 #define KM 1000.0 // kilometer
89
90 #define FNT "DejaVuMonoSans"
91 #define CMAX 256 // gray scale
92
93 // ==========================================================================
94
95 int main() {
96
97 FILE *pngout ; // file handle for png output frame
98 char labstring[80] ; // used for labelling
99 int icolor ; // index for computing colors
100 int i ; // general counter
101 int fs ; // large font size
102 double satx[NSAT] ; // satellite x position, east positive
103 double saty[NSAT] ; // satellite y position, out positive
104 double satp[NSAT] ; // satellite phase in radians
105 double sata[NSAT] ; // satellite amplitude
106
107 double arrx[A_NUM] ; // array center position
108 double arry[A_NUM] ; // array center position
109
110 double pi2 = 8.0*atan(1.0) ; // 2*pi, 6.283
111 double d2r = pi2/360.0 ; // degrees to radians
112 double r2d = 360.0/pi2 ; // radians to degrees
113
114 // -------------------
115 // coordinates - z, north south, out of plane of drawing
116 // x, orbit tangential distance from center thinsat array
117
118 double targetr = (RE+TARGETA)*KM ; // target radius in meters
119 double targetn = d2r*TARGETN ; // target latitude in radians north
120 double targetx = 0.0 ; // target x coordinate
121 double targety = targetr*cos(targetn); // target y coordinate
122 double targetz = targetr*sin(targetn); // target z coordinate
123
124 int idx ; //
125 double cscale = (double) CMAX ; //
126
127 // -------------------
128 double k = pi2/LAMBDA ; // wavenumber
129 double re = RE * KM ; // earth radius in meters
130 double r288 = R288 * KM ; // m288 orbit radius in meters
131 int smid = ( S_NUM - 1 ) / 2 ;
132 int amid = ( A_NUM - 1 ) / 2 ;
133 double s_num = (double) S_NUM ;
134 double a_sp = d2r* A_SPACE ; // array spacing in radians
135 double s_sp = S_WIDTH /((S_NUM-1.0)*r288) ; // thinsat spacing rad.
136
137 // -------------------
138 // compute position and phase of each thinsat
139
140 double asum = 0.0 ; // normalize amplitude
141 int acnt, scnt ;
142
143 for( acnt=0 ; acnt < A_NUM ; acnt++ ) {
144 double a_c = (acnt - amid)*a_sp ; //
145 double a_0 = a_c - smid*s_sp ; // array starting angle
146 int aidx = acnt*S_NUM ;
147 arrx[acnt] = r288 * sin(a_c) ;
148 arry[acnt] = r288 * cos(a_c) ;
149
150 for( scnt=0 ; scnt < S_NUM ; scnt++ ) {
151 double s_an = a_0 + scnt * s_sp ; // orbit angle for thinsat
152 int idx = aidx + scnt ; // array index for thinsat
153
154 satx[idx] = r288 * sin(s_an) ; //
155 saty[idx] = r288 * cos(s_an) ; //
156
157 double dx = satx[idx]-targetx ; //
158 double dy = saty[idx]-targety ; //
159 double dz = targetz ; //
160 double d = sqrt( dx*dx + dy*dy + dz*dz ) ; // distance
161 satp[idx] = k*d ; // phase
162 sata[idx] = 1.0 ; // uniform amplitude for now
163 asum += sata[idx]/d ; // summed amplitude
164 }
165 }
166
167 double anorm = 1.0 / asum ;
168
169 // set up array and define colors -----------------------------------------
170
171 gdImagePtr im = gdImageCreateTrueColor(XSIZE, YSIZE );
172
173 // allocate standard colors
174 int black = gdImageColorAllocate(im, 0, 0, 0);
175 int white = gdImageColorAllocate(im, 255, 255, 255);
176 int sun1 = gdImageColorAllocate(im, 51, 51, 102);
177 int red = gdImageColorAllocate(im, 255, 0, 0);
178 int green = gdImageColorAllocate(im, 0, 255, 0);
179 int dgreen= gdImageColorAllocate(im, 0, 128, 0);
180 int blue = gdImageColorAllocate(im, 0, 0, 255);
181 int gray = gdImageColorAllocate(im, 128, 128, 128);
182 int dgray = gdImageColorAllocate(im, 48, 48, 48);
183 int trans = gdImageColorAllocate(im, 1, 1, 1);
184 int gcolor[CMAX] ; // color map numbers
185
186 // allocate gray scale array
187 for( icolor=0 ; icolor<CMAX ; icolor++ ) {
188 int col = 255-icolor ;
189 gcolor[icolor] = gdImageColorAllocate(im, col, col, col);
190 }
191
192 // white background
193
194 gdImageFilledRectangle( im, 0, 0, XSIZE-1, YSIZE-1, white );
195 gdImageSetThickness( im, LWID );
196
197 // slide label -----------------------------------------------------------
198
199 fs= 80 ;
200 gdImageStringFT( im, NULL, // imagespace, bounding box
201 black , FNT, fs, 0.0, // color, font, fontsize, angle
202 fs , fs, // x, y
203 "Phase-synced\nmulti-beam\nradar array" ); // text
204 /*
205 gdImageStringFT( im, NULL, // imagespace, bounding box
206 black , FNT, fs, 0.0, // color, font, fontsize, angle
207 fs , 6*fs/2, // x, y
208 "multi-beam" ); // text
209
210 gdImageStringFT( im, NULL, // imagespace, bounding box
211 black , FNT, fs, 0.0, // color, font, fontsize, angle
212 fs , 9*fs/2, // x, y
213 "radar array" ); // text
214 */
215
216 // draw earth and orbit segments -----------------------------------------
217
218 double ex0 = (double) EX0 ;
219 double ex1 = (double) EX1 ;
220 double eyc = (double) EYC ;
221 double exc = 0.5*(ex0+ex1) ;
222 int exci = (int) exc ;
223 double earc = 0.5*a_sp*(double)A_NUM ;
224 int eari = (int) exc ;
225 double ewide= 2.0*r288*sin(earc);
226 double escal= (ex1-ex0)/ewide ;
227
228 // draw earth
229
230 int ew = (int)( 2.0*re*escal );
231
232 gdImageArc(im, exci, EYC, ew, ew, 0, 360, dgreen );
233
234 // draw orbit fragment
235 int edeg = (int)( r2d*earc );
236 int eow = (int)( escal*r288*2.0 );
237 gdImageArc(im, exci, EYC, eow, eow, 270-edeg, 270+edeg, blue );
238
239 // target position
240 int etx = exc + targetx * escal ;
241 int ety = eyc - targety * escal ;
242
243 // all arrays and beams
244 for( i = 0 ; i < A_NUM ; i++ ) {
245 int eax = (int)( exc + escal*arrx[i] );
246 int eay = (int)( eyc - escal*arry[i] );
247 gdImageLine( im, eax, eay, etx, ety, black ) ;
248 }
249
250 // these will be used for drawing expansion lines
251 int px0 = exci ;
252 int py0 = ety ;
253
254 // draw a tilted orbital ellipse representing a circular orbit
255
256 double einc = d2r*EINC ;
257 double eaop = d2r*EAOP ;
258 double eang ;
259
260 for( eang = 0.0; eang <= pi2; eang += pi2 / 2048.0 ) {
261 double ex00 = targetr * sin( eang ) ;
262 double ey00 = targetr * cos( eang ) ;
263 double ez00 = 0.0 ;
264
265 double ex01 = ex00 ;
266 double ey01 = ey00*cos(einc) - ez00*sin(einc);
267 double ez01 = ey00*sin(einc) + ez00*cos(einc);
268
269 double ex02 = ex01*cos(eaop) - ey01*sin(eaop);
270 double ey02 = ex01*sin(eaop) + ey01*cos(eaop);
271 double ez02 = ez01 ;
272
273 double er02 = sqrt( ex02*ex02 + ey02*ey02 );
274
275 if( ( er02 > re ) || ( ez02 > 0.0 ) ) {
276 int ex = (int)( exc + ex02 * escal ) ;
277 int ey = (int)( eyc + ey02 * escal ) ;
278
279 gdImageFilledEllipse( im, ex, ey, LWID, LWID, red );
280 }
281 }
282
283 // HUGE KLUDGE ---------------------------------------------------------
284 // fake loop to change what were defines
285 int define, OX0, OX1 ;
286 double SCANX, LABX, EXAG, USCALE, TUNITS ;
287 char UNITS[16], TSCALE[16];
288 char str[128] ;
289
290
291 // SCANX LABX EXAG USCALE TUNITS OX0 OX1 UNITS TSCALE
292 for( define = 0 ; define <= 1 ; define++ ) {
293 if( define == 0 ) { strcpy( str, D0 ) ; }
294 if( define == 1 ) { strcpy( str, D1 ) ; }
295
296 sscanf( str, "%lf %lf %lf %lf %lf %d %d %s %s", &SCANX,
297 &LABX, &EXAG, &USCALE, &TUNITS, &OX0, &OX1, UNITS, TSCALE );
298
299 printf( "\n%s\n", str );
300 printf( "SCANX=%9.3e ", SCANX );
301 printf( "LABX=%9.3e ", LABX );
302 printf( "EXAG=%9.3e ", EXAG );
303 printf( "USCALE=%9.3e\n", USCALE);
304 printf( "TUNITS=%9.3e ", USCALE);
305 printf( "OX0=%6d ", OX0 );
306 printf( "OX1=%6d ", OX1 );
307 printf( "UNITS=%s ", UNITS );
308 printf( "TSCALE=%s\n", TSCALE );
309
310 // ------------------------------------------------------------------------
311
312 // plot scaling - compute x and y from pixel value
313
314 double ox0 = (double) OX0 ; // output display left
315 double ox1 = (double) OX1 ; // output display right
316 double oy0 = (double) OY0 ; // output display top
317 double oy1 = (double) OY1 ; // output display bottom
318 int ox, oy ; // output grid indices
319
320 double scany = SCANX * (oy1-oy0)/(ox1-ox0);
321
322 double x_d = SCANX / (ox1-ox0) ;
323 double x_0 = targetx - ( x_d * ox0 + 0.5*SCANX ) ;
324 double y_d = -x_d ; // pixels increase, y gets smaller
325 double y_0 = targety - ( y_d * oy0 - 0.5*scany ) ;
326
327 fs= 60 ;
328
329 gdImageStringFT( im, NULL, // imagespace, bounding box
330 black , FNT, fs, 0.0, // color, font, fontsize, angle
331 OX0 , oy1+3*fs, // x, y
332 UNITS ); // text
333
334 // now compute pixel array
335
336 double dz2 = targetz * targetz ;
337
338 #ifdef IMAGE
339 for( ox=ox0 ; ox<ox1 ; ox++ ) {
340 double grx = x_d*ox + x_0 ;
341
342 for( oy=oy0 ; oy<oy1 ; oy++ ) {
343 double gry = y_d*oy + y_0 ;
344
345 // compute z distance and phase of each source to plane
346 double asum = 0.0 ;
347
348 for( idx = 0; idx < NSAT ; idx ++ ) {
349 double dx = satx[idx] - grx ;
350 double dy = saty[idx] - gry ;
351 double dist = sqrt( dx*dx + dy*dy + dz2 );
352 asum += sata[idx]* cos( k*dist - satp[idx] ) / dist ;
353 }
354 asum *= anorm ; // normalize
355
356 // power/intensity computation
357 int mag = (int)( EXAG*cscale*asum*asum ) ;
358 if( mag > cscale ) { mag = cscale ; }
359
360 gdImageSetPixel( im, ox, oy, gcolor[mag] ) ;
361 }
362 }
363 #endif
364
365 // now compute path through field
366
367 gdImageLine( im, OX0-10, GY, OX1+10, GY, black );
368
369 double gox = -1.0 ;
370 double goy = -1.0 ;
371
372 for( ox=ox0 ; ox<ox1 ; ox++ ) {
373 double grx = x_d*ox + x_0 ;
374 double gry = OSLOPE*grx + targety ;
375 oy = (int)( ( gry-y_0)/ y_d ) ;
376
377 gdImageFilledEllipse( im, ox, oy, 3, 3, red );
378
379 // compute z distance and phase of each source to plane
380 double asum = 0.0 ;
381
382 for( idx = 0; idx < NSAT ; idx ++ ) {
383 double dx = satx[idx] - grx ;
384 double dy = saty[idx] - gry ;
385 double dist = sqrt( dx*dx + dy*dy + dz2 );
386 asum += cos( k*dist - satp[idx] ) / dist ;
387 }
388 asum *= anorm ; // normalize
389 int gpy = GY - (int)( ((double)GGAIN)*asum );
390 if( ox == ox0 ) {
391 gox = ox ;
392 goy = gpy ;
393 }
394 gdImageLine( im, gox, goy, ox, gpy, red ) ;
395 gox = ox ;
396 goy = gpy ;
397 }
398
399 fs = 60;
400 gdImageStringFT( im, NULL, // imagespace, bounding
401 red, FNT, fs, 0.0, // color, font, fontsize, angle
402 0, GY, // x, y
403 " reflected signal" ); // the text
404
405 // draw box ------------------------------------------------------------
406
407 gdImageStringFT( im, NULL, // imagespace, bounding
408 black, FNT, fs, 0.0, // color, font, fontsize, angle
409 OX0+1.0*fs, OY0-fs/2, // x, y
410 "Standing Wave Intensity" ); // the text
411
412 gdImageLine( im, OX0, OY0, OX1, OY0, black );
413 gdImageLine( im, OX1, OY1, OX1, OY0, black );
414 gdImageLine( im, OX1, OY1, OX0, OY1, black );
415 gdImageLine( im, OX0, OY0, OX0, OY1, black );
416
417 // draw ticmarks on edges ----------------------------------------------
418
419 int tx = (int) (0.001+0.5*(SCANX/LABX)) ; // tickmark +/- count
420 int ty = (int) (0.001+0.5*(scany/LABX)) ; // tickmark +/- count
421 double ts = (ox1-ox0)*LABX/SCANX ; // tickmark step size
422
423 // draw x tickmarks
424
425 int oxc = (OX0+OX1)/2 ; // center of output display
426 int oyc = (OY0+OY1)/2 ; // center of output display
427
428 double labcm = LABX*USCALE ;
429
430 for( i = -tx ; i <= tx ; i++ ) {
431 ox = oxc - (int) ( ts*i ) ;
432 gdImageLine( im, ox, OY1, ox, OY1-5, black );
433 gdImageLine( im, ox, OY0, ox, OY0+5, black );
434
435 double lval = -labcm * (double) i ;;
436 sprintf( labstring, "%4.0f", lval );
437 gdImageStringFT( im, NULL, // imagespace, bounding
438 black, FNT, fs, 0.0, // color, font, fontsize, angle
439 ox-(5*fs)/2, OY1+(3*fs)/2, // x, y
440 labstring ); // the text
441 }
442 gdImageLine( im, oxc, OY1, oxc, OY1-10, black );
443 gdImageLine( im, oxc, OY0, oxc, OY0+10, black );
444
445 // draw y tickmarks
446
447 for( i = -ty; i <= ty ; i++ ) {
448 oy = oyc - (int) ( ts*i ) ;
449 gdImageLine( im, OX1, oy, OX1-5, oy, black );
450 gdImageLine( im, OX0, oy, OX0+5, oy, black );
451
452 double lval = labcm * (double) i ;;
453 sprintf( labstring, "%4.0f", lval );
454 gdImageStringFT( im, NULL, // imagespace, bounding
455 black, FNT, fs, 0.0, // color, font, fontsize, angle
456 OX0-4*fs, oy+fs/2, // x, y,
457 labstring ); // the text
458 }
459
460 // time estimate
461
462 double torb = TUNITS * SCANX * cos(einc) * cos(eaop) / VORB ;
463
464 sprintf( labstring, "%5.0f %s at%5.0f m/s", torb, TSCALE, VORB );
465
466 gdImageStringFT( im, NULL, // imagespace, bounding
467 black, FNT, fs, 0.0, // color, font, fontsize, angle
468 OX0+2*fs, GT, // x, y
469 labstring ); // the text
470
471 // draw expansion lines
472 gdImageSetThickness( im, LNAR );
473 gdImageLine( im, px0, py0, OX0, OY0, black );
474 gdImageLine( im, px0, py0, OX0, OY1, black );
475 px0 = oxc ;
476 py0 = oyc ;
477 gdImageSetThickness( im, LWID );
478 }
479
480
481 // output the frame ------------------------------------------------------
482
483 pngout = fopen( PNGFMT, "wb");
484 gdImagePngEx( im, pngout, 1 );
485 gdImageDestroy(im);
486 fclose(pngout);
487 }
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.