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