1 /*
2 * xnec2c - GTK2-based version of nec2c, the C translation of NEC2
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU Library General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
17 */
18
19 /*
20 * gnuplot.c
21 *
22 * gnuplot routines for xnec2c
23 */
24
25 #include "gnuplot.h"
26 #include "shared.h"
27
28 /*-----------------------------------------------------------------------*/
29
30 /* Save_FreqPlots_Gnuplot_Data()
31 *
32 * Saves frequency plots data to a file for gnuplot
33 */
34 void
Save_FreqPlots_Gnuplot_Data(char * filename)35 Save_FreqPlots_Gnuplot_Data( char *filename )
36 {
37 /* Abort if plot data not available */
38 if( isFlagClear(FREQ_LOOP_DONE) )
39 return;
40
41 /* Used to calculate net gain */
42 double Zr, Zo, Zi;
43
44 /* Open gplot file, abort on error */
45 FILE *fp = NULL;
46 if( !Open_File(&fp, filename, "w") )
47 return;
48
49 /* Plot max gain vs frequency, if possible */
50 if( isFlagSet(PLOT_GMAX) && isFlagSet(ENABLE_RDPAT) )
51 {
52 int nth, nph, idx, pol;
53 gboolean no_fbr;
54
55 double
56 gmax, /* Max gain buffer */
57 netgain, /* Viewer direction net gain buffer */
58 gdir_phi, /* Direction in phi of gain */
59 fbratio; /* Front to back ratio */
60
61 /* Find max gain and direction, F/B ratio */
62 no_fbr = FALSE;
63 netgain = 0;
64
65 /* Polarization type and impedance */
66 pol = calc_data.pol_type;
67 Zo = calc_data.zo;
68
69 /* Save data for all frequency steps that were used */
70 fprintf( fp, _("# Gain and F/B Ratio vs Frequency\n") );
71 for( idx = 0; idx <= calc_data.lastf; idx++ )
72 {
73 double fbdir;
74 int fbidx, mgidx;
75
76 /* Index to gtot buffer where max gain
77 * occurs for given polarization type */
78 mgidx = rad_pattern[idx].max_gain_idx[pol];
79
80 /* Max gain for given polarization type */
81 gmax = rad_pattern[idx].gtot[mgidx] +
82 Polarization_Factor(pol, idx, mgidx);
83
84 /* Net gain if selected */
85 if( isFlagSet(PLOT_NETGAIN) )
86 {
87 Zr = impedance_data.zreal[idx];
88 Zi = impedance_data.zimag[idx];
89 netgain = gmax + 10*log10(4*Zr*Zo/(pow(Zr+Zo,2)+pow(Zi,2)));
90 }
91
92 /* Radiation angle/phi where max gain occurs */
93 gdir_phi = rad_pattern[idx].max_gain_phi[pol];
94
95 /* Find F/B direction in theta */
96 fbdir = 180.0 - rad_pattern[idx].max_gain_tht[pol];
97 if( fpat.dth == 0.0 )
98 nth = 0;
99 else
100 nth = (int)( fbdir/fpat.dth + 0.5 );
101
102 /* If the antenna is modelled over ground, then use the same
103 * theta as the max gain direction, relying on phi alone to
104 * take us to the back. Patch supplied by Rik van Riel AB1KW
105 */
106 if( (nth >= fpat.nth) || (nth < 0) )
107 {
108 fbdir = rad_pattern[idx].max_gain_tht[pol];
109 if( fpat.dth == 0.0 )
110 nth = 0;
111 else
112 nth = (int)( fbdir/fpat.dth + 0.5 );
113 }
114
115 /* Find F/B direction in phi */
116 fbdir = gdir_phi + 180.0;
117 if( fbdir >= 360.0 ) fbdir -= 360.0;
118 nph = (int)( fbdir/fpat.dph + 0.5 );
119
120 /* No F/B calc. possible if no phi step at +180 from max gain */
121 if( (nph >= fpat.nph) || (nph < 0) )
122 no_fbr = TRUE;
123
124 /* Index to gtot buffer for gain in back direction */
125 fbidx = nth + nph*fpat.nth;
126
127 /* Front to back ratio */
128 fbratio = pow( 10.0, gmax / 10.0 );
129 fbratio /= pow( 10.0,
130 (rad_pattern[idx].gtot[fbidx] +
131 Polarization_Factor(pol, idx, fbidx)) / 10.0 );
132 fbratio = 10.0 * log10( fbratio );
133
134 if( no_fbr && isFlagClear(PLOT_NETGAIN) ) /* Plot max gain only */
135 fprintf( fp, "%13.6E %10.3E\n", save.freq[idx], gmax );
136 else if( isFlagSet(PLOT_NETGAIN) ) /* Plot max gain and net gain */
137 fprintf( fp, "%13.6E %10.3E %10.3E\n", save.freq[idx], gmax, netgain );
138 else if( !no_fbr ) /* Plot max gain and F/B ratio */
139 fprintf( fp, "%13.6E %10.3E %10.3E\n", save.freq[idx], gmax, fbratio );
140 } /* for( idx = 0; idx < calc_data.lastf; idx++ ) */
141
142 /* Plot gain direction in phi and theta */
143 if( isFlagSet(PLOT_GAIN_DIR) )
144 {
145 fprintf( fp, "\n\n" );
146 fprintf( fp, _("# Direction of gain in theta and phi\n") );
147 for( idx = 0; idx < calc_data.lastf; idx++ )
148 {
149 double gdir_tht; /* Direction in theta of gain */
150
151 /* Radiation angle/phi where max gain occurs */
152 gdir_tht = 90.0 - rad_pattern[idx].max_gain_tht[pol];
153 gdir_phi = rad_pattern[idx].max_gain_phi[pol];
154 fprintf( fp, "%13.6E %10.3E %10.3E\n", save.freq[idx], gdir_tht, gdir_phi );
155 } /* for( idx = 0; idx < calc_data.lastf; idx++ ) */
156 } /* if( isFlagSet(PLOT_GAIN_DIR) ) */
157
158 fprintf( fp, "\n\n" );
159 } /* if( isFlagSet(PLOT_GMAX) && isFlagSet(ENABLE_RDPAT) ) */
160
161 /* Plot VSWR vs freq */
162 if( isFlagSet(PLOT_VSWR) )
163 {
164 int idx;
165 double vswr, gamma;
166 double zrpro2, zrmro2, zimag2;
167
168 /* Calculate VSWR */
169 fprintf( fp, _("# VSWR vs Frequency\n") );
170 for(idx = 0; idx <= calc_data.lastf; idx++ )
171 {
172 zrpro2 = impedance_data.zreal[idx] + calc_data.zo;
173 zrpro2 *= zrpro2;
174 zrmro2 = impedance_data.zreal[idx] - calc_data.zo;
175 zrmro2 *= zrmro2;
176 zimag2 = impedance_data.zimag[idx] * impedance_data.zimag[idx];
177 gamma = sqrt( (zrmro2 + zimag2)/(zrpro2 + zimag2) );
178 vswr = (1+gamma)/(1-gamma);
179 if( vswr > 10.0 ) vswr = 10.0;
180 fprintf( fp, "%13.6E %10.3E\n", save.freq[idx], vswr );
181 }
182
183 fprintf( fp, "\n\n" );
184 } /* if( isFlagSet(PLOT_VSWR) ) */
185
186 /* Plot z-real and z-imag */
187 if( isFlagSet(PLOT_ZREAL_ZIMAG) )
188 {
189 int idx;
190 fprintf( fp, _("# Z real & Z imaginary vs Frequency\n") );
191 for(idx = 0; idx <= calc_data.lastf; idx++ )
192 fprintf( fp, "%13.6E %10.3E %10.3E\n",
193 save.freq[idx], impedance_data.zreal[idx], impedance_data.zimag[idx] );
194
195 fprintf( fp, "\n\n" );
196 } /* if( isFlagSet(PLOT_ZREAL_ZIMAG) ) */
197
198 /* Plot z-magn and z-phase */
199 if( isFlagSet(PLOT_ZMAG_ZPHASE) )
200 {
201 int idx;
202 fprintf( fp, _("# Z magnitude & Z phase vs Frequency\n") );
203 for(idx = 0; idx <= calc_data.lastf; idx++ )
204 fprintf( fp, "%13.6E %10.3E %10.3E\n",
205 save.freq[idx], impedance_data.zmagn[idx], impedance_data.zphase[idx] );
206 } /* if( isFlagSet(PLOT_ZREAL_ZIMAG) ) */
207
208 fclose(fp);
209 } /* Save_FreqPlots_Gnuplot_Data() */
210
211 /*-----------------------------------------------------------------------*/
212
213 /* Save_RadPattern_Gnuplot_Data()
214 *
215 * Saves radiation pattern data to a file for gnuplot
216 */
217 void
Save_RadPattern_Gnuplot_Data(char * filename)218 Save_RadPattern_Gnuplot_Data( char *filename )
219 {
220 int idx, npts; /* Number of points to plot */
221
222 /* Scale factor ref, for normalizing field strength values */
223 double dr;
224
225 double
226 fx, fy, fz, /* Co-ordinates of "free" end of field lines */
227 fscale; /* Scale factor for equalizing field line segments */
228
229 FILE *fp = NULL;
230
231 /* Draw near field pattern if possible */
232 if( isFlagSet(ENABLE_NEAREH) && near_field.valid )
233 {
234 /* Open gplot file, abort on error */
235 if( !Open_File(&fp, filename, "w") )
236 return;
237
238 /* Reference for scale factor used in
239 * normalizing field strength values */
240 if( fpat.near ) /* Spherical co-ordinates */
241 dr = (double)fpat.dxnr;
242 /* Rectangular co-ordinates */
243 else dr = sqrt(
244 (double)fpat.dxnr * (double)fpat.dxnr +
245 (double)fpat.dynr * (double)fpat.dynr +
246 (double)fpat.dznr * (double)fpat.dznr )/1.75;
247
248 npts = fpat.nrx * fpat.nry * fpat.nrz;
249
250 /*** Draw Near E Field ***/
251 if( isFlagSet(DRAW_EFIELD) && (fpat.nfeh & NEAR_EFIELD) )
252 {
253 fprintf( fp, _("# Near E field\n") );
254 /* Write e-field out to file [DJS] */
255 for( idx = 0; idx < npts; idx++ )
256 {
257 fscale = dr / near_field.er[idx];
258 fx = near_field.px[idx] + near_field.erx[idx] * fscale;
259 fy = near_field.py[idx] + near_field.ery[idx] * fscale;
260 fz = near_field.pz[idx] + near_field.erz[idx] * fscale;
261
262 /* Print as x, y, z, dx, dy, dz for gnuplot */
263 fprintf( fp, "%f %f %f %f %f %f\n",
264 near_field.px[idx],
265 near_field.py[idx],
266 near_field.pz[idx],
267 fx - near_field.px[idx],
268 fy - near_field.py[idx],
269 fz - near_field.pz[idx] );
270 }
271 } /* if( isFlagSet(DRAW_EFIELD) */
272
273 /*** Draw Near H Field ***/
274 if( isFlagSet(DRAW_HFIELD) && (fpat.nfeh & NEAR_HFIELD) )
275 {
276 fprintf( fp, _("# Near H field\n") );
277 /* Write h-field out to file [DJS] */
278 for( idx = 0; idx < npts; idx++ )
279 {
280 fscale = dr / near_field.hr[idx];
281 fx = near_field.px[idx] + near_field.hrx[idx] * fscale;
282 fy = near_field.py[idx] + near_field.hry[idx] * fscale;
283 fz = near_field.pz[idx] + near_field.hrz[idx] * fscale;
284
285 /* Print as x, y, z, dx, dy, dz for gnuplot */
286 fprintf( fp, "%f %f %f %f %f %f\n",
287 near_field.px[idx],
288 near_field.py[idx],
289 near_field.pz[idx],
290 fx - near_field.px[idx],
291 fy - near_field.py[idx],
292 fz - near_field.pz[idx] );
293 }
294 } /* if( isFlagSet(DRAW_HFIELD) && (fpat.nfeh & NEAR_HFIELD) ) */
295
296 /*** Draw Poynting Vector ***/
297 if( isFlagSet(DRAW_POYNTING) &&
298 (fpat.nfeh & NEAR_EFIELD) &&
299 (fpat.nfeh & NEAR_HFIELD) )
300 {
301 int ipv;
302 static size_t mreq = 0;
303
304 /* Co-ordinates of Poynting vectors */
305 static double *pov_x = NULL, *pov_y = NULL;
306 static double *pov_z = NULL, *pov_r = NULL;
307
308 /* Range of Poynting vector values,
309 * its max and min and log of max/min */
310 static double pov_max = 0;
311
312 /* Allocate on new near field matrix size */
313 if( mreq != (size_t)npts * sizeof( double ) )
314 {
315 mreq = (size_t)npts * sizeof( double );
316 mem_realloc( (void **)&pov_x, mreq, "in draw_radiation.c" );
317 mem_realloc( (void **)&pov_y, mreq, "in draw_radiation.c" );
318 mem_realloc( (void **)&pov_z, mreq, "in draw_radiation.c" );
319 mem_realloc( (void **)&pov_r, mreq, "in draw_radiation.c" );
320 }
321
322 /* Calculate Poynting vector and its max and min */
323 fprintf( fp, _("# Poynting Vector\n") );
324 for( idx = 0; idx < npts; idx++ )
325 {
326 pov_max = 0;
327 for( ipv = 0; ipv < npts; ipv++ )
328 {
329 pov_x[ipv] =
330 near_field.ery[ipv] * near_field.hrz[ipv] -
331 near_field.hry[ipv] * near_field.erz[ipv];
332 pov_y[ipv] =
333 near_field.erz[ipv] * near_field.hrx[ipv] -
334 near_field.hrz[ipv] * near_field.erx[ipv];
335 pov_z[ipv] =
336 near_field.erx[ipv] * near_field.hry[ipv] -
337 near_field.hrx[ipv] * near_field.ery[ipv];
338 pov_r[ipv] = sqrt(
339 pov_x[ipv] * pov_x[ipv] +
340 pov_y[ipv] * pov_y[ipv] +
341 pov_z[ipv] * pov_z[ipv] );
342 if( pov_max < pov_r[ipv] )
343 pov_max = pov_r[ipv];
344
345 } /* for( ipv = 0; ipv < npts; ipv++ ) */
346
347 /* Scale factor for each field point, to make
348 * near field direction lines equal-sized */
349 fscale = dr / pov_r[idx];
350
351 /* Scaled field values are used to set one end of a
352 * line segment that represents direction of field.
353 * The other end is set by the field point co-ordinates */
354 fx = near_field.px[idx] + pov_x[idx] * fscale;
355 fy = near_field.py[idx] + pov_y[idx] * fscale;
356 fz = near_field.pz[idx] + pov_z[idx] * fscale;
357
358 /* Print as x, y, z, dx, dy, dz for gnuplot */
359 fprintf( fp, "%f %f %f %f %f %f\n",
360 near_field.px[idx],
361 near_field.py[idx],
362 near_field.pz[idx],
363 fx - near_field.px[idx],
364 fy - near_field.py[idx],
365 fz - near_field.pz[idx] );
366 } /* for( idx = 0; idx < npts; idx++ ) */
367
368 } /* if( isFlagSet(DRAW_POYNTING) ) */
369 } /* if( isFlagSet(ENABLE_NEAREH) && near_field.valid ) */
370
371 /* Save radiation pattern data if possible */
372 if( isFlagSet(ENABLE_RDPAT) && (calc_data.fstep >= 0) )
373 {
374 int
375 nth, /* Theta step count */
376 nph; /* Phi step count */
377
378 /* Frequency step and polarization type */
379 int fstep = calc_data.fstep;
380
381 /* Theta and phi angles defining a rad pattern point
382 * and distance of its projection from xyz origin */
383 double theta, phi, r;
384
385 /* theta and phi step in rads */
386 double dth = (double)fpat.dth * (double)TA;
387 double dph = (double)fpat.dph * (double)TA;
388
389 /* Open gplot file, abort on error */
390 if( !Open_File(&fp, filename, "w") )
391 return;
392 fprintf( fp, _("# Radiation Pattern") );
393
394 /* Distance of rdpattern point nearest to xyz origin */
395 /*** Convert radiation pattern values
396 * to points in 3d space in x,y,z axis ***/
397 phi = (double)fpat.phis * (double)TA; /* In rads */
398
399 /* Step phi angle */
400 idx = 0;
401 for( nph = 0; nph < fpat.nph; nph++ )
402 {
403 theta = (double)fpat.thets * (double)TA; /* In rads */
404
405 /* Step theta angle */
406 for( nth = 0; nth < fpat.nth; nth++ )
407 {
408 double x, y, z;
409
410 /* Distance of pattern point from the xyz origin */
411 r = Scale_Gain(
412 rad_pattern[fstep].gtot[idx], fstep, idx );
413
414 /* Distance of point's projection on xyz axis, from origin */
415 z = r * cos(theta);
416 r *= sin(theta);
417 x = r * cos(phi);
418 y = r * sin(phi);
419
420 /* Print to file */
421 fprintf( fp, "%10.3E %10.3E %10.3E\n", x, y, z );
422
423 /* Step theta in rads */
424 theta += dth;
425 idx++;
426 } /* for( nth = 0; nth < fpat.nth; nth++ ) */
427
428 /* Step phi in rads */
429 phi += dph;
430 } /* for( nph = 0; nph < fpat.nph; nph++ ) */
431
432
433 } /* if( isFlagSet(ENABLE_RDPAT) && (calc_data.fstep >= 0) ) */
434
435 if( fp != NULL ) fclose(fp);
436 } /* Save_RadPattern_Gnuplot_Data() */
437
438 /*-----------------------------------------------------------------------*/
439
440 /* Save_Struct_Gnuplot_Data()
441 *
442 * Saves antenna structure data for gnuplot
443 */
444 void
Save_Struct_Gnuplot_Data(char * filename)445 Save_Struct_Gnuplot_Data( char *filename )
446 {
447 FILE *fp = NULL;
448
449 /* Open gplot file, abort on error */
450 if( !Open_File(&fp, filename, "w") )
451 return;
452
453 /* Output if patch segs and no new input pending */
454 if( data.m && isFlagClear(INPUT_PENDING) )
455 {
456 int idx, m2;
457
458 /* Output segments data */
459 fprintf( fp, _("# structure patch segmenets\n") );
460
461 /* Output first segment outside loop to enable separation of wires */
462 fprintf( fp, "%10.3E %10.3E %10.3E\n%10.3E %10.3E %10.3E\n",
463 (double)data.px1[0], (double)data.py1[0], (double)data.pz1[0],
464 (double)data.px2[0], (double)data.py2[0], (double)data.pz2[0] );
465
466 /* Start from second segment and check for connection of ends */
467 m2 = data.m * 2;
468 for( idx = 1; idx < m2; idx++ )
469 {
470 fprintf( fp, "%10.3E %10.3E %10.3E\n%10.3E %10.3E %10.3E\n",
471 (double)data.px1[idx], (double)data.py1[idx], (double)data.pz1[idx],
472 (double)data.px2[idx], (double)data.py2[idx], (double)data.pz2[idx] );
473 } /* for( idx = 1; idx < m2; idx++ ) */
474
475 fprintf( fp, "\n\n" );
476 } /* if( data.m && isFlagSet(INPUT_PENDING) ) */
477
478 /* Output if wire segs and no new input pending */
479 if( data.n && isFlagClear(INPUT_PENDING) )
480 {
481 int idx;
482
483 /* Output segments data */
484 fprintf( fp, _("# structure wire segmenets\n") );
485
486 /* Output first segment outside loop to enable separation of wires */
487 fprintf( fp, "%10.3E %10.3E %10.3E\n%10.3E %10.3E %10.3E\n",
488 (double)data.x1[0], (double)data.y1[0], (double)data.z1[0],
489 (double)data.x2[0], (double)data.y2[0], (double)data.z2[0] );
490
491 /* Start from second segment and check for connection of ends */
492 for( idx = 1; idx < data.n; idx++ )
493 {
494 /* Leave a 2-line gap to next segment */
495 if( (data.icon1[idx] == 0) || (data.icon1[idx] == (idx+1)) )
496 fprintf( fp, "\n\n" );
497 fprintf( fp, "%10.3E %10.3E %10.3E\n%10.3E %10.3E %10.3E\n",
498 (double)data.x1[idx], (double)data.y1[idx], (double)data.z1[idx],
499 (double)data.x2[idx], (double)data.y2[idx], (double)data.z2[idx] );
500 /* Leave a 2-line gap to next segment */
501 if( (data.icon2[idx] == 0) || (data.icon2[idx] == (idx+1)) )
502 fprintf( fp, "\n\n" );
503
504 } /* for( idx = 1; idx < data.n; idx++ ) */
505 } /* if( data.n && isFlagSet(INPUT_PENDING) ) */
506
507 fclose( fp );
508 } /* Save_Struct_Gnuplot_Data() */
509
510 /*-----------------------------------------------------------------------*/
511
512