1 /* XNECVIEW - a program for visualizing NEC2 input and output data
2 *
3 * Copyright (C) 2000-2006, Pieter-Tjerk de Boer -- pa3fwm@amsat.org
4 *
5 * Distributed on the conditions of version 2 of the GPL: see the files
6 * README and COPYING, which accompany this source file.
7 *
8 * This module contains code for drawing plots of several quantities
9 * like impedance, SWR and gain as a function of frequency.
10 *
11 */
12
13 #include <stdio.h>
14 #include <math.h>
15 #include <float.h>
16 #include <string.h>
17 #include <gdk/gdk.h>
18
19 #include "xnecview.h"
20
21
22 int win2sizex,win2sizey; /* size of window in pixels */
23
24 int plot2_swr=1; /* show the SWR graph? */
25 int plot2_maxgain=1; /* show the maxgain and front/back graph? */
26 int plot2_vgain=0; /* show the vgain graph? */
27 int plot2_z=0; /* show the impedance graph? */
28 int plot2_z2=0; /* show the phi(z)/abs(z) graph? */
29 int plot2_dir=0; /* show the direction-of-maximum-gain graph? */
30
31 double r0=R0; /* reference impedance for SWR calculation */
32
33
34
35
fixrange1(double * mi,double * ma,int * np)36 void fixrange1(double *mi, double *ma, int *np)
37 /* mi and ma point to minimum and maximum value of a range of values to
38 be plotted, and np to the maximum acceptable number of subdivision.
39 This function tries to modify the minimum and maximum and the number
40 of subdivision such that the resulting grid lines are at "round" numbers.
41 */
42 {
43 double d,e;
44 double a;
45 double newmin,newmax;
46 int i;
47 int n=*np;
48 /*
49 static double acceptable[]={10, 5.0, 2.5, 2.0, 1.0, -1};
50 */
51 static double acceptable[]={10, 5.0, 2.0, 1.0, -1};
52
53 if (*ma==*mi) {
54 if (*mi>0) {*mi=0; *ma=2* *ma;}
55 else if (*mi<0) {*mi=2* *mi; *ma=0;}
56 else {*mi=-10; *ma=10;}
57 }
58 d=(*ma-*mi)/n;
59 e=1.0;
60 while (e<d) e*=10;
61 while (e>d) e/=10;
62 a=d/e;
63 i=0;
64 while (acceptable[i]>a) i++;
65 if (acceptable[i]==-1) i--;
66 i++;
67 do {
68 i--;
69 if (i<0) {
70 e*=10;
71 i=0;
72 while (acceptable[i+1]>0) i++;
73 }
74 a=acceptable[i];
75 d=a*e;
76 newmin = d*floor(*mi/d);
77 newmax = d*ceil(*ma/d);
78 n = (int)((newmax-newmin)/d+0.5);
79 } while (n>*np);
80 *np=n;
81 *mi=newmin;
82 *ma=newmax;
83 }
84
85
fixrange2(double * mi1,double * ma1,double * mi2,double * ma2,int * np)86 void fixrange2(double *mi1, double *ma1, double *mi2, double *ma2, int *np)
87 /* like fixrange2(), but for two (vertical) axes simultaneously */
88 {
89 static double acceptable[]={100.0, 50.0, 25.0, 20.0, 10.0, 5.0, 2.5, 2.0, 1.0, 0.5, 0.25, 0.2, 0.1, 0.05, 0.025, 0.02, 0.01, -1};
90 double a,d,e1,e2,s;
91 int n=*np;
92 int i1,i2;
93 int i,j;
94 int ibest,jbest;
95 int n1[5],n2[5];
96 double x1[4],x2[4];
97 double best;
98
99 if (*ma1==*mi1) {
100 if (*mi1>0) {*mi1=0; *ma1=2* *ma1;}
101 else if (*mi1<0) {*mi1=2* *mi1; *ma1=0;}
102 else {*mi1=-10; *ma1=10;}
103 }
104 d=(*ma1-*mi1)/n; /* d is the ideal, but usually not acceptable, stepsize, for axis 1 */
105 *ma1-=0.00001*d; /* prevent rounding errors from causing a boundary of say 1000 to be seen as slightly larger than say 10 steps of 100 each */
106 *mi1+=0.00001*d; /* idem */
107 d-=0.00001*d;
108 e1=1.0;
109 while (e1<d) e1*=10;
110 while (e1>d) e1/=10; /* e1 is the appropriate power of 10 to scale the steps, for axis 1 */
111 a=d/e1;
112 i1=0;
113 while (acceptable[i1+1]>=a) i1++; /* i1 is the index in the acceptable[] array of the highest acceptable stepsize, for axis 1 */
114 for (i=0;i<4;i++) { /* consider this and the next 3 lower stepsizes: */
115 s = e1*acceptable[i1-i];
116 n1[i] = ceil(*ma1/s) - floor(*mi1/s) ; /* minimum number of acceptable steps */
117 x1[i] = (*ma1-*mi1) / s; /* "usage factor": how many of these steps does the data cover? */
118 }
119
120 /* same calculations for axis 2 */
121 if (*ma2==*mi2) {
122 if (*mi2>0) {*mi2=0; *ma2=2* *ma2;}
123 else if (*mi2<0) {*mi2=2* *mi2; *ma2=0;}
124 else {*mi2=-10; *ma2=10;}
125 }
126 d=(*ma2-*mi2)/n;
127 *ma2-=0.00001*d;
128 *mi2+=0.00001*d;
129 d-=0.00001*d;
130 e2=1.0;
131 while (e2<d) e2*=10;
132 while (e2>d) e2/=10;
133 a=d/e2;
134 i2=0;
135 while (acceptable[i2+1]>=a) i2++;
136 for (i=0;i<4;i++) {
137 s = e2*acceptable[i2-i];
138 n2[i] = ceil(*ma2/s) - floor(*mi2/s) ;
139 x2[i] = (*ma2-*mi2) / s;
140 }
141
142 /* search for best combination: the combination for which the data covers as large a fraction of both axes as possible */
143 best=0;
144 ibest=jbest=0;
145 for (i=0;i<4;i++)
146 for (j=0;j<4;j++) {
147 double x;
148 int n;
149 n = n1[i];
150 if (n2[j]>n) n=n2[j];
151 x = (x1[i]/n) * (x2[j]/n);
152 if (x>best*1.1 || (x>best && n>=*np)) { best=x; ibest=i; jbest=j; *np=n; }
153 }
154
155 n = *np;
156 i1-=ibest;
157 i2-=jbest;
158 s = e1*acceptable[i1];
159 *mi1 = s*floor(*mi1/s);
160 *ma1 = *mi1+n*s;
161 s = e2*acceptable[i2];
162 *mi2 = s*floor(*mi2/s);
163 *ma2 = *mi2+n*s;
164 }
165
166
167 double minf,maxf;
168 int xleft, xright;
169
170 #define idxOK(idx,ne) (idx>=0 && ( !ONLY_IF_RP(idx) || ne->rp ) && ne->d[idx]>-DBL_MAX)
171
freqplot(int idx1,int idx2,int idx1a,int idx2a,char * title1,char * title2,char * title,GdkColor * color1,GdkColor * color2,double ybotf,double ytopf)172 void freqplot(
173 int idx1, /* index in neco[].d[] of quantity for left axis */
174 int idx2, /* index in neco[].d[] of quantity for right axis */
175 int idx1a, /* index in neco[].d[] of second quantity for left axis (dotted line) */
176 int idx2a, /* index in neco[].d[] of second quantity for right axis (dotted line) */
177 char *title1, char *title2, /* titles for left and right */
178 char *title, /* center title */
179 GdkColor *color1, GdkColor *color2, /* colours for both curves */
180 double ybotf, double ytopf /* vertical position; 0...1 = top...bottom of window */
181 )
182 {
183 int ybot, ytop;
184 int i;
185 double min1,max1, min2, max2;
186 NECoutput *ne;
187 int ntx,nty;
188 int xx1,xx2,yy1,yy2;
189 int xx1a,xx2a,yy1a,yy2a;
190
191 /* choose the corner points of the graph area */
192 ybot = ybotf*win2sizey - fontheight;
193 ytop = ytopf*win2sizey + fontheight;
194 xleft = 5*fontheight;
195 xright = win2sizex - 5*fontheight;
196
197 /* find the ranges */
198 minf=maxf=neco[0].f;
199 min1=min2=DBL_MAX;
200 max1=max2=-DBL_MAX;
201 for (i=0, ne=neco; i<numneco; i++, ne++) {
202 if (ne->f < minf) minf=ne->f;
203 if (ne->f > maxf) maxf=ne->f;
204 if (idxOK(idx1,ne)) {
205 if (ne->d[idx1] < min1) min1=ne->d[idx1];
206 if (ne->d[idx1] > max1) max1=ne->d[idx1];
207 }
208 if (idxOK(idx1a,ne)) {
209 if (ne->d[idx1a] < min1) min1=ne->d[idx1a];
210 if (ne->d[idx1a] > max1) max1=ne->d[idx1a];
211 }
212 if (idxOK(idx2,ne)) {
213 if (ne->d[idx2] < min2) min2=ne->d[idx2];
214 if (ne->d[idx2] > max2) max2=ne->d[idx2];
215 }
216 if (idxOK(idx2a,ne)) {
217 if (ne->d[idx2a] < min2) min2=ne->d[idx2a];
218 if (ne->d[idx2a] > max2) max2=ne->d[idx2a];
219 }
220 }
221 if (min1>max1) { idx1=-1; idx1a=-1; }
222 if (min2>max2) { idx2=-1; idx2a=-1; }
223
224 /* extend the ranges to have 'round' numbers at each division */
225 ntx=(xright-xleft)/fontheight/3;
226 fixrange1(&minf,&maxf,&ntx);
227 nty = (ybot-ytop)/fontheight;
228 if (idx1>=0) {
229 if (idx1==neco_zr) {
230 if (max1 > 20*r0) max1 = 20*r0;
231 }
232 }
233 if (idx2>=0) {
234 if (idx2==neco_zi || idx2==neco_zabs) {
235 if (max2 > 20*r0 && min2 < 20*r0) max2 = 20*r0;
236 if (min2 < -20*r0 && max2 > -20*r0) min2 = -20*r0;
237 }
238 }
239 if (idx1==neco_swr) { min1=0; if (max1>10) max1=9; else max1-=1; }
240 if (idx2==neco_swr) { min2=0; if (max2>10) max2=9; else max2-=1; }
241 if (idx1<0 && idx1a<0) {
242 if (idx2<0 && idx2a<0) return;
243 fixrange1(&min2,&max2,&nty);
244 } else {
245 if (idx2<0 && idx2a<0) fixrange1(&min1,&max1,&nty);
246 else fixrange2(&min1,&max1,&min2,&max2,&nty);
247 }
248 if (idx1==neco_swr) { min1+=1; max1+=1; }
249 if (idx2==neco_swr) { min2+=1; max2+=1; }
250
251
252
253 /* macros for converting from "real" values to screen coordinates */
254 #define sx(f) (((f)-minf)/(maxf-minf)*(xright-xleft)+xleft)
255 #define sy1(f) (((f)-min1)/(max1-min1)*(ytop-ybot)+ybot)
256 #define sy2(f) (((f)-min2)/(max2-min2)*(ytop-ybot)+ybot)
257
258 SetLineAttributes(0, GDK_LINE_SOLID, GDK_CAP_ROUND, GDK_JOIN_ROUND);
259 /* vertical grid lines and associated labels */
260 for (i=0; i<=ntx; i++) {
261 double f;
262 int x;
263 char s[20];
264 f=minf+(maxf-minf)*((double)i)/ntx;
265 x=sx(f);
266 if (i>0 && i<ntx) {
267 SetForeground(&c_scale);
268 DrawLine(x,ybot,x,ytop);
269 }
270 SetForeground(&c_axis);
271 sprintf(s,"%g",f);
272 DrawString(x,ybot+1,s,0.5,1);
273 if (i==ntx) {
274 int l;
275 l=strlen(s);
276 sprintf(s," MHz"+(15-(l+1)/2));
277 DrawString(x,ybot+1,s,0,1);
278 }
279 }
280 if (idx1<0) { min1=1; max1=2; }
281 /* horizontal grid lines and associated labels */
282 for (i=0; i<=nty; i++) {
283 double f;
284 int y;
285 char s[20];
286 f=min1+(max1-min1)*((double)i)/nty;
287 if (fabs(f/(max1-min1))<0.1/nty) f=0;
288 y=sy1(f);
289 if (i>0 && i<nty) {
290 SetForeground(&c_scale);
291 DrawLine(xleft,y,xright,y);
292 }
293 if (idx1>=0) {
294 sprintf(s,"%g ",f);
295 SetForeground(color1);
296 DrawString(xleft,y,s,1,0.5);
297 }
298 if (idx2>=0) {
299 f=min2+(max2-min2)*((double)i)/nty;
300 if (fabs(f/(max2-min2))<0.1/nty) f=0;
301 y=sy2(f);
302 sprintf(s," %g",f);
303 SetForeground(color2);
304 DrawString(xright,y,s,0,0.5);
305 }
306 }
307 SetForeground(&c_axis);
308 /* border around the graph */
309 DrawLine(xleft,ybot,xright,ybot);
310 DrawLine(xleft,ytop,xright,ytop);
311 DrawLine(xleft,ybot,xleft,ytop);
312 DrawLine(xright,ybot,xright,ytop);
313
314 /* title(s) */
315 if (title) {
316 SetForeground(&c_axis);
317 DrawString((xleft+xright)/2, ytop-1, title, 0.5,0);
318 }
319 if (title1) {
320 SetForeground(color1);
321 DrawString(xleft-fontheight, ytop-1, title1, 0,0);
322 }
323 if (title2) {
324 SetForeground(color2);
325 DrawString(xright+fontheight, ytop-1, title2, 1,0);
326 }
327
328 /* the actual data points and connecting lines */
329 SetClipRectangle(xleft-2,ytop,xright+2,ybot);
330 xx1=xx2=yy1=yy2=-1;
331 xx1a=xx2a=yy1a=yy2a=-1;
332 for (i=0, ne=neco; i<numneco; i++, ne++) {
333 int x,y;
334 x=sx(ne->f);
335 if (idx1a>=0 || idx2a>=0) SetLineAttributes(0, GDK_LINE_ON_OFF_DASH, GDK_CAP_ROUND, GDK_JOIN_ROUND);
336 if (idxOK(idx1a,ne)) {
337 y = sy1(ne->d[idx1a]);
338 SetForeground(color1);
339 if (numneco < win2sizex / 4) {
340 DrawLine(x-2,y-2,x+2,y-2);
341 DrawLine(x-2,y+2,x+2,y+2);
342 DrawLine(x-2,y-2,x-2,y+2);
343 DrawLine(x+2,y-2,x+2,y+2);
344 }
345 if (xx1a!=-1) DrawLine(xx1a,yy1a,x,y);
346 xx1a=x; yy1a=y;
347 }
348 if (idxOK(idx2a,ne)) {
349 y = sy2(ne->d[idx2a]);
350 SetForeground(color2);
351 if (numneco < win2sizex / 4) {
352 DrawLine(x-3,y,x,y-3);
353 DrawLine(x-3,y,x,y+3);
354 DrawLine(x+3,y,x,y-3);
355 DrawLine(x+3,y,x,y+3);
356 }
357 if (xx2a!=-1) DrawLine(xx2a,yy2a,x,y);
358 xx2a=x; yy2a=y;
359 }
360 if (idx1a>=0 || idx2a>=0) SetLineAttributes(0, GDK_LINE_SOLID, GDK_CAP_ROUND, GDK_JOIN_ROUND);
361 if (idxOK(idx1,ne)) {
362 y = sy1(ne->d[idx1]);
363 SetForeground(color1);
364 if (numneco < win2sizex / 4) {
365 DrawLine(x-2,y-2,x+2,y-2);
366 DrawLine(x-2,y+2,x+2,y+2);
367 DrawLine(x-2,y-2,x-2,y+2);
368 DrawLine(x+2,y-2,x+2,y+2);
369 }
370 if (xx1!=-1) DrawLine(xx1,yy1,x,y);
371 xx1=x; yy1=y;
372 }
373 if (idxOK(idx2,ne)) {
374 y = sy2(ne->d[idx2]);
375 SetForeground(color2);
376 if (numneco < win2sizex / 4) {
377 DrawLine(x-3,y,x,y-3);
378 DrawLine(x-3,y,x,y+3);
379 DrawLine(x+3,y,x,y-3);
380 DrawLine(x+3,y,x,y+3);
381 }
382 if (xx2!=-1) DrawLine(xx2,yy2,x,y);
383 xx2=x; yy2=y;
384 }
385 }
386 SetClipRectangle(0,0,win2sizex,win2sizey);
387 }
388
389
xfreq(int x)390 double xfreq(int x)
391 {
392 return ((double)x-xleft)/(xright-xleft)*(maxf-minf)+minf;
393 }
394
395
freqx(double f)396 int freqx(double f)
397 {
398 return sx(f);
399 }
400
401
freqindex(double f)402 int freqindex(double f)
403 {
404 double d,dbest;
405 int i,ibest;
406
407 ibest=-1;
408 dbest=DBL_MAX;
409 for (i=0;i<numneco;i++) {
410 if ( !neco[i].rp && !neco[i].cu && !neco[i].nf ) continue;
411 d=fabs(f-neco[i].f);
412 if (d<dbest) {
413 ibest=i;
414 dbest=d;
415 }
416 }
417 return ibest;
418 }
419
420
421
draw_all2(int onlyvgain)422 void draw_all2(int onlyvgain)
423 {
424 int n;
425 double size,step,y;
426
427 n = plot2_swr+plot2_z+plot2_z2+plot2_maxgain+plot2_dir+plot2_vgain;
428 size = 1.0/(n+(n-1)*0.05);
429 step = 1.05*size;
430 y = 1.0;
431
432 if (!onlyvgain) ClearWindow();
433 if (plot2_z2) {
434 if (!onlyvgain) freqplot(neco_zphi,neco_zabs,-1,-1,"phi(Z) [deg]","[ohm] |Z|","impedance",&c_exci,&c_load,y,y-size);
435 y-=step;
436 }
437 if (plot2_z) {
438 if (!onlyvgain) freqplot(neco_zr,neco_zi,-1,-1,"real [ohm]","[ohm] imag","impedance",&c_exci,&c_load,y,y-size);
439 y-=step;
440 }
441 if (plot2_swr) {
442 if (!onlyvgain) freqplot(neco_swr,-1,-1,-1,"SWR",NULL,NULL,&c_wire,NULL,y,y-size);
443 y-=step;
444 }
445 if (plot2_dir) {
446 if (!onlyvgain) {
447 if (polarization==POLnone || polarization==POLcolour) freqplot(neco_phi,neco_theta,-1,-1,"phi [deg]","[deg] theta","direction of maximum gain",&c_exci,&c_load,y,y-size);
448 else freqplot(Neco_polphi+Neco_gsize*polarization,
449 Neco_poltheta+Neco_gsize*polarization,
450 neco_phi,neco_theta,"phi","theta","direction of maximum gain",&c_exci,&c_load,y,y-size);
451 }
452 y-=step;
453 }
454 if (plot2_maxgain) {
455 if (!onlyvgain) {
456 if (polarization==POLnone || polarization==POLcolour) freqplot(neco_maxgain,neco_fb,-1,-1,"gain [dB]","[dB] f/b","in direction of maximum gain",&c_gain,&c_surf,y,y-size);
457 else freqplot(Neco_polgain+Neco_gsize*polarization,
458 Neco_polfb1+Neco_gsize*polarization,
459 neco_maxgain,
460 Neco_polfb2+Neco_gsize*polarization,
461 "gain","f/b","in direction of maximum gain",&c_gain,&c_surf,y,y-size);
462 }
463 y-=step;
464 }
465 if (plot2_vgain) {
466 if (onlyvgain) ClearRectangle(0,(y-size)*win2sizey,win2sizex,y*win2sizey);
467 if (polarization==POLnone || polarization==POLcolour) freqplot(neco_vgain,neco_vfb,-1,-1,"gain [dB]","[dB] f/b","in direction toward viewer",&c_gain,&c_surf,y,y-size);
468 else freqplot(neco_vgain2,neco_vfb,neco_vgain,neco_vfb2,"gain [dB]","[dB] f/b","in direction toward viewer",&c_gain,&c_surf,y,y-size);
469 y-=step;
470 }
471 out->Complete();
472 }
473
474
475