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 /******* Translated to the C language by N. Kyriazis 20 Aug 2003 ******
20
21 Program NEC(input,tape5=input,output,tape11,tape12,tape13,tape14,
22 tape15,tape16,tape20,tape21)
23
24 Numerical Electromagnetics Code (NEC2) developed at Lawrence
25 Livermore lab., Livermore, CA. (contact G. Burke at 415-422-8414
26 for problems with the NEC code. For problems with the vax implem-
27 entation, contact J. Breakall at 415-422-8196 or E. Domning at 415
28 422-5936)
29 file created 4/11/80.
30
31 ***********Notice**********
32 This computer code material was prepared as an account of work
33 sponsored by the United States government. Neither the United
34 States nor the United States Department Of Energy, nor any of
35 their employees, nor any of their contractors, subcontractors,
36 or their employees, makes any warranty, express or implied, or
37 assumes any legal liability or responsibility for the accuracy,
38 completeness or usefulness of any information, apparatus, product
39 or process disclosed, or represents that its use would not infringe
40 privately-owned rights.
41
42 ***********************************************************************/
43
44 #include "radiation.h"
45 #include "shared.h"
46
47 /* Radiation pattern data */
48 /*-----------------------------------------------------------------------*/
49
50 /* ffld calculates the far zone radiated electric fields, */
51 /* the factor exp(j*k*r)/(r/lamda) not included */
ffld(double thet,double phi,complex double * eth,complex double * eph)52 void ffld( double thet, double phi,
53 complex double *eth, complex double *eph )
54 {
55 int k, i, ip, jump;
56 double phx, phy, roz, rozs, thx, thy, thz, rox, roy;
57 double tthet=0.0, darg=0.0, omega, el, sill, top, bot, a;
58 double too, boo, b, c, d, rr, ri, arg, dr, rfl, rrz;
59 complex double cix=CPLX_00, ciy=CPLX_00, ciz=CPLX_00, ccx=CPLX_00;
60 complex double ccy=CPLX_00, ccz=CPLX_00, exa, cdp;
61 complex double zrsin, rrv=CPLX_00, rrh=CPLX_00, rrv1=CPLX_00;
62 complex double rrh1=CPLX_00, rrv2=CPLX_00, rrh2=CPLX_00;
63 complex double tix, tiy, tiz, zscrn, ex=CPLX_00;
64 complex double ey=CPLX_00, ez=CPLX_00, gx, gy, gz;
65
66 phx= -sin( phi);
67 phy= cos( phi);
68 roz= cos( thet);
69 rozs= roz;
70 thx= roz* phy;
71 thy= -roz* phx;
72 thz= -sin( thet);
73 rox= -thz* phy;
74 roy= thz* phx;
75
76 jump = FALSE;
77 if( data.n != 0)
78 {
79 /* loop for structure image if any */
80 /* calculation of reflection coeffecients */
81 for( k = 0; k < gnd.ksymp; k++ )
82 {
83 if( k != 0 )
84 {
85 /* for perfect ground */
86 if( gnd.iperf == 1)
87 {
88 rrv=-CPLX_10;
89 rrh=-CPLX_10;
90 }
91 else
92 {
93 /* for infinite planar ground */
94 zrsin= csqrt(1.0- gnd.zrati* gnd.zrati* thz* thz);
95 rrv=-( roz- gnd.zrati* zrsin)/( roz+ gnd.zrati* zrsin);
96 rrh=( gnd.zrati* roz- zrsin)/( gnd.zrati* roz+ zrsin);
97 } /* if( gnd.iperf == 1) */
98
99 /* for the cliff problem, two reflection coefficients calculated */
100 if( gnd.ifar > 1)
101 {
102 rrv1= rrv;
103 rrh1= rrh;
104 tthet= tan( thet);
105
106 if( gnd.ifar != 4)
107 {
108 zrsin= csqrt(1.0- gnd.zrati2* gnd.zrati2* thz* thz);
109 rrv2=-( roz- gnd.zrati2* zrsin)/( roz+ gnd.zrati2* zrsin);
110 rrh2=( gnd.zrati2* roz- zrsin)/( gnd.zrati2* roz+ zrsin);
111 darg= -TP*2.0* gnd.ch* roz;
112 }
113 } /* if( gnd.ifar > 1) */
114
115 roz= -roz;
116 ccx= cix;
117 ccy= ciy;
118 ccz= ciz;
119
120 } /* if( k != 0 ) */
121
122 cix=CPLX_00;
123 ciy=CPLX_00;
124 ciz=CPLX_00;
125
126 /* loop over structure segments */
127 for( i = 0; i < data.n; i++ )
128 {
129 omega=-( rox* data.cab[i] + roy* data.sab[i]+ roz* data.salp[i]);
130 el= PI* data.si[i];
131 sill= omega* el;
132 top= el+ sill;
133 bot= el- sill;
134
135 if( fabs( omega) >= 1.0e-7)
136 a=2.0* sin( sill)/ omega;
137 else
138 a=(2.0- omega* omega* el* el/3.0)* el;
139
140 if( fabs( top) >= 1.0e-7)
141 too= sin( top)/ top;
142 else
143 too=1.0- top* top/6.0;
144
145 if( fabs( bot) >= 1.0e-7)
146 boo= sin( bot)/ bot;
147 else
148 boo=1.0- bot* bot/6.0;
149
150 b= el*( boo- too);
151 c= el*( boo+ too);
152 rr= a* crnt.air[i]+ b* crnt.bii[i]+ c* crnt.cir[i];
153 ri= a* crnt.aii[i]- b* crnt.bir[i]+ c* crnt.cii[i];
154 arg= TP*( data.x[i]* rox+ data.y[i]* roy+ data.z[i]* roz);
155
156 if( (k != 1) || (gnd.ifar < 2) )
157 {
158 /* summation for far field integral */
159 exa= cmplx( cos( arg), sin( arg))* cmplx( rr, ri);
160 cix= cix+ exa* data.cab[i];
161 ciy= ciy+ exa* data.sab[i];
162 ciz= ciz+ exa* data.salp[i];
163 continue;
164 }
165
166 /* calculation of image contribution */
167 /* in cliff and ground screen problems */
168
169 /* specular point distance */
170 dr= data.z[i]* tthet;
171 d= dr* phy+ data.x[i];
172 if( gnd.ifar == 2)
173 {
174 if(( gnd.cl- d) > 0.0)
175 {
176 rrv= rrv1;
177 rrh= rrh1;
178 }
179 else
180 {
181 rrv= rrv2;
182 rrh= rrh2;
183 arg= arg+ darg;
184 }
185 } /* if( gnd.ifar == 2) */
186 else
187 {
188 d= sqrt( d*d + (data.y[i]-dr*phx)*(data.y[i]-dr*phx) );
189 if( gnd.ifar == 3)
190 {
191 if(( gnd.cl- d) > 0.0)
192 {
193 rrv= rrv1;
194 rrh= rrh1;
195 }
196 else
197 {
198 rrv= rrv2;
199 rrh= rrh2;
200 arg= arg+ darg;
201 }
202 } /* if( gnd.ifar == 3) */
203 else
204 {
205 if(( gnd.scrwl- d) >= 0.0)
206 {
207 /* radial wire ground screen reflection coefficient */
208 d= d+ gnd.t2;
209 zscrn= gnd.t1* d* log( d/ gnd.t2);
210 zscrn=( zscrn* gnd.zrati)/( ETA* gnd.zrati+ zscrn);
211 zrsin= csqrt(1.0- zscrn* zscrn* thz* thz);
212 rrv=( roz+ zscrn* zrsin)/(- roz+ zscrn* zrsin);
213 rrh=( zscrn* roz+ zrsin)/( zscrn* roz- zrsin);
214 } /* if(( gnd.scrwl- d) < 0.0) */
215 else
216 {
217 if( gnd.ifar == 4)
218 {
219 rrv= rrv1;
220 rrh= rrh1;
221 } /* if( gnd.ifar == 4) */
222 else
223 {
224 if( gnd.ifar == 5)
225 d= dr* phy+ data.x[i];
226
227 if(( gnd.cl- d) > 0.0)
228 {
229 rrv= rrv1;
230 rrh= rrh1;
231 }
232 else
233 {
234 rrv= rrv2;
235 rrh= rrh2;
236 arg= arg+ darg;
237 } /* if(( gnd.cl- d) > 0.0) */
238
239 } /* if( gnd.ifar == 4) */
240
241 } /* if(( gnd.scrwl- d) < 0.0) */
242
243 } /* if( gnd.ifar == 3) */
244
245 } /* if( gnd.ifar == 2) */
246
247 /* contribution of each image segment modified by */
248 /* reflection coef, for cliff and ground screen problems */
249 exa= cmplx( cos( arg), sin( arg))* cmplx( rr, ri);
250 tix= exa* data.cab[i];
251 tiy= exa* data.sab[i];
252 tiz= exa* data.salp[i];
253 cdp=( tix* phx+ tiy* phy)*( rrh- rrv);
254 cix= cix+ tix* rrv+ cdp* phx;
255 ciy= ciy+ tiy* rrv+ cdp* phy;
256 ciz= ciz- tiz* rrv;
257
258 } /* for( i = 0; i < n; i++ ) */
259
260 if( k == 0 )
261 continue;
262
263 /* calculation of contribution of
264 * structure image for infinite ground */
265 if( gnd.ifar < 2)
266 {
267 cdp=( cix* phx+ ciy* phy)*( rrh- rrv);
268 cix= ccx+ cix* rrv+ cdp* phx;
269 ciy= ccy+ ciy* rrv+ cdp* phy;
270 ciz= ccz- ciz* rrv;
271 }
272 else
273 {
274 cix= cix+ ccx;
275 ciy= ciy+ ccy;
276 ciz= ciz+ ccz;
277 }
278
279 } /* for( k=0; k < gnd.ksymp; k++ ) */
280
281 if( data.m > 0)
282 jump = TRUE;
283 else
284 {
285 *eth = ( cix * thx + ciy * thy + ciz * thz ) * CONST3;
286 *eph = ( cix * phx + ciy * phy ) * CONST3;
287 return;
288 }
289
290 } /* if( n != 0) */
291
292 if( ! jump )
293 {
294 cix=CPLX_00;
295 ciy=CPLX_00;
296 ciz=CPLX_00;
297 }
298
299 /* electric field components */
300 roz= rozs;
301 rfl=-1.0;
302 for( ip = 0; ip < gnd.ksymp; ip++ )
303 {
304 rfl= -rfl;
305 rrz= roz* rfl;
306 fflds( rox, roy, rrz, &crnt.cur[data.n], &gx, &gy, &gz);
307
308 if( ip != 1 )
309 {
310 ex= gx;
311 ey= gy;
312 ez= gz;
313 continue;
314 }
315
316 if( gnd.iperf == 1)
317 {
318 gx= -gx;
319 gy= -gy;
320 gz= -gz;
321 }
322 else
323 {
324 rrv= csqrt(1.0- gnd.zrati* gnd.zrati* thz* thz);
325 rrh= gnd.zrati* roz;
326 rrh=( rrh- rrv)/( rrh+ rrv);
327 rrv= gnd.zrati* rrv;
328 rrv=-( roz- rrv)/( roz+ rrv);
329 *eth=( gx* phx+ gy* phy)*( rrh- rrv);
330 gx= gx* rrv+ *eth* phx;
331 gy= gy* rrv+ *eth* phy;
332 gz= gz* rrv;
333
334 } /* if( gnd.iperf == 1) */
335
336 ex= ex+ gx;
337 ey= ey+ gy;
338 ez= ez- gz;
339
340 } /* for( ip = 0; ip < gnd.ksymp; ip++ ) */
341
342 ex = ex + cix * CONST3;
343 ey = ey + ciy * CONST3;
344 ez = ez + ciz * CONST3;
345
346 *eth = ex * thx + ey * thy + ez * thz;
347 *eph = ex * phx + ey * phy;
348
349 return;
350 }
351
352 /*-----------------------------------------------------------------------*/
353
354 /* calculates the xyz components of the electric */
355 /* field due to surface currents */
fflds(double rox,double roy,double roz,complex double * scur,complex double * ex,complex double * ey,complex double * ez)356 void fflds( double rox, double roy, double roz,
357 complex double *scur, complex double *ex,
358 complex double *ey, complex double *ez )
359 {
360 double *xs, *ys, *zs, *s;
361 int j, i, k;
362 double arg;
363 complex double ct;
364
365 xs = data.px;
366 ys = data.py;
367 zs = data.pz;
368 s = data.pbi;
369
370 *ex=CPLX_00;
371 *ey=CPLX_00;
372 *ez=CPLX_00;
373
374 i= -1;
375 for( j = 0; j < data.m; j++ )
376 {
377 i++;
378 arg= TP*( rox* xs[i]+ roy* ys[i]+ roz* zs[i]);
379 ct= cmplx( cos( arg)* s[i], sin( arg)* s[i]);
380 k=3*j;
381 *ex += scur[k ]* ct;
382 *ey += scur[k+1]* ct;
383 *ez += scur[k+2]* ct;
384 }
385
386 ct= rox* *ex+ roy* *ey+ roz* *ez;
387 *ex= CONST4*( ct* rox- *ex);
388 *ey= CONST4*( ct* roy- *ey);
389 *ez= CONST4*( ct* roz- *ez);
390
391 return;
392 }
393
394 /*-----------------------------------------------------------------------*/
395
396 /* gfld computes the radiated field including ground wave. */
gfld(double rho,double phi,double rz,complex double * eth,complex double * epi,complex double * erd,complex double ux,int ksymp)397 void gfld( double rho, double phi, double rz,
398 complex double *eth, complex double *epi,
399 complex double *erd, complex double ux, int ksymp )
400 {
401 int i, k;
402 double b, r, thet, arg, phx, phy, rx, ry;
403 double dx, dy, dz, rix, riy, rhs, rhp;
404 double rhx, rhy, calp, cbet, sbet, cph;
405 double sph, el, rfl, riz, thx, thy, thz;
406 double rxyz, rnx, rny, rnz, omega, sill;
407 double top, bot, a, too, boo, c, rr, ri;
408 complex double cix, ciy, ciz, exa, erv;
409 complex double ezv, erh, eph, ezh, ex, ey;
410
411 r= sqrt( rho*rho+ rz*rz );
412 if( (ksymp == 1) || (cabs(ux) > .5) || (r > 1.0e5) )
413 {
414 /* computation of space wave only */
415 if( rz >= 1.0e-20)
416 thet= atan( rho/ rz);
417 else
418 thet= PI*.5;
419
420 ffld( thet, phi, eth, epi);
421 arg= -TP* r;
422 exa= cmplx( cos( arg), sin( arg))/ r;
423 *eth= *eth* exa;
424 *epi= *epi* exa;
425 *erd=CPLX_00;
426 return;
427 } /* if( (ksymp == 1) && (cabs(ux) > .5) && (r > 1.0e5) ) */
428
429 /* computation of space and ground waves. */
430 gwav.u= ux;
431 gwav.u2= gwav.u* gwav.u;
432 phx= -sin( phi);
433 phy= cos( phi);
434 rx= rho* phy;
435 ry= -rho* phx;
436 cix=CPLX_00;
437 ciy=CPLX_00;
438 ciz=CPLX_00;
439
440 /* summation of field from individual segments */
441 for( i = 0; i < data.n; i++ )
442 {
443 dx= data.cab[i];
444 dy= data.sab[i];
445 dz= data.salp[i];
446 rix= rx- data.x[i];
447 riy= ry- data.y[i];
448 rhs= rix* rix+ riy* riy;
449 rhp= sqrt( rhs);
450
451 if( rhp >= 1.0e-6)
452 {
453 rhx= rix/ rhp;
454 rhy= riy/ rhp;
455 }
456 else
457 {
458 rhx=1.0;
459 rhy=0.0;
460 }
461
462 calp=1.0- dz* dz;
463 if( calp >= 1.0e-6)
464 {
465 calp= sqrt( calp);
466 cbet= dx/ calp;
467 sbet= dy/ calp;
468 cph= rhx* cbet+ rhy* sbet;
469 sph= rhy* cbet- rhx* sbet;
470 }
471 else
472 {
473 cph= rhx;
474 sph= rhy;
475 }
476
477 el= PI* data.si[i];
478 rfl=-1.0;
479
480 /* integration of (current)*(phase factor)
481 * over segment and image for constant,
482 * sine, and cosine current distributions */
483 for( k = 0; k < 2; k++ )
484 {
485 rfl= -rfl;
486 riz= rz- data.z[i]* rfl;
487 rxyz= sqrt( rix* rix+ riy* riy+ riz* riz);
488 rnx= rix/ rxyz;
489 rny= riy/ rxyz;
490 rnz= riz/ rxyz;
491 omega=-( rnx* dx+ rny* dy+ rnz* dz* rfl);
492 sill= omega* el;
493 top= el+ sill;
494 bot= el- sill;
495
496 if( fabs( omega) >= 1.0e-7)
497 a=2.0* sin( sill)/ omega;
498 else
499 a=(2.0- omega* omega* el* el/3.0)* el;
500
501 if( fabs( top) >= 1.0e-7)
502 too= sin( top)/ top;
503 else
504 too=1.0- top* top/6.0;
505
506 if( fabs( bot) >= 1.0e-7)
507 boo= sin( bot)/ bot;
508 else
509 boo=1.0- bot* bot/6.0;
510
511 b= el*( boo- too);
512 c= el*( boo+ too);
513 rr= a* crnt.air[i] +
514 b* crnt.bii[i]+ c* crnt.cir[i];
515 ri= a* crnt.aii[i] -
516 b* crnt.bir[i]+ c* crnt.cii[i];
517 arg= TP*( data.x[i] *
518 rnx+ data.y[i]* rny+ data.z[i]* rnz* rfl);
519 exa= cmplx( cos( arg), sin( arg))* cmplx( rr, ri)/ TP;
520
521 if( k != 1 )
522 {
523 gwav.xx1= exa;
524 gwav.r1= rxyz;
525 gwav.zmh= riz;
526 continue;
527 }
528
529 gwav.xx2= exa;
530 gwav.r2= rxyz;
531 gwav.zph= riz;
532
533 } /* for( k = 0; k < 2; k++ ) */
534
535 /* call subroutine to compute the field */
536 /* of segment including ground wave. */
537 gwave( &erv, &ezv, &erh, &ezh, &eph);
538 erh= erh* cph* calp+ erv* dz;
539 eph= eph* sph* calp;
540 ezh= ezh* cph* calp+ ezv* dz;
541 ex= erh* rhx- eph* rhy;
542 ey= erh* rhy+ eph* rhx;
543 cix= cix+ ex;
544 ciy= ciy+ ey;
545 ciz= ciz+ ezh;
546
547 } /* for( i = 0; i < n; i++ ) */
548
549 arg= -TP* r;
550 exa= cmplx( cos( arg), sin( arg));
551 cix= cix* exa;
552 ciy= ciy* exa;
553 ciz= ciz* exa;
554 rnx= rx/ r;
555 rny= ry/ r;
556 rnz= rz/ r;
557 thx= rnz* phy;
558 thy= -rnz* phx;
559 thz= -rho/ r;
560 *eth= cix* thx+ ciy* thy+ ciz* thz;
561 *epi= cix* phx+ ciy* phy;
562 *erd= cix* rnx+ ciy* rny+ ciz* rnz;
563
564 return;
565 }
566
567 /*-----------------------------------------------------------------------*/
568
569 /* compute radiation pattern, gain, normalized gain */
rdpat(void)570 void rdpat(void)
571 {
572 int kth, kph, isens;
573 double prad, gcon, gcop;
574 double phi, pha, thet;
575 double tha, ethm2, ethm;
576 double etha, ephm2, ephm, epha, tilta, emajr2, eminr2;
577 double dfaz, axrat, dfaz2, cdfaz, tstor1=0.0, tstor2;
578 double gnmn, stilta, gnmj, gnv, gnh, gtot;
579 complex double eth, eph, erd;
580 int idx, pol; /* Gain buffer and pol type index */
581 double gain;
582
583
584 if( gnd.ifar != 4 )
585 {
586 gnd.cl= fpat.clt/ data.wlam;
587 gnd.ch= fpat.cht/ data.wlam;
588 gnd.zrati2= csqrt(1.0/ cmplx(
589 fpat.epsr2,- fpat.sig2* data.wlam*59.96));
590 }
591
592 /* Calculate radiation pattern data */
593 /*** For applied voltage excitation ***/
594 if( (fpat.ixtyp == 0) || (fpat.ixtyp == 5) )
595 {
596 gcop= data.wlam* data.wlam* 2.0* PI/(376.73* fpat.pinr);
597 prad= fpat.pinr- fpat.ploss- fpat.pnlr;
598 gcon= gcop;
599 if( fpat.ipd != 0)
600 gcon *= fpat.pinr/ prad;
601 }
602 /*** For elementary current source ***/
603 else if( fpat.ixtyp == 4)
604 {
605 fpat.pinr=394.510* calc_data.xpr6*
606 calc_data.xpr6* data.wlam* data.wlam;
607 gcop= data.wlam* data.wlam*2.0* PI/(376.73* fpat.pinr);
608 prad= fpat.pinr- fpat.ploss- fpat.pnlr;
609 gcon= gcop;
610 if( fpat.ipd != 0)
611 gcon= gcon* fpat.pinr/ prad;
612 }
613 /*** Incident field source ***/
614 else gcon=4.0* PI/(1.0+ calc_data.xpr6* calc_data.xpr6);
615
616 phi = fpat.phis - fpat.dph;
617
618 /*** Save radiation pattern data ***/
619 /* Prime max and min gains and index */
620 for( pol = 0; pol < NUM_POL; pol++ )
621 {
622 rad_pattern[calc_data.fstep].max_gain[pol] = -10000.0;
623 rad_pattern[calc_data.fstep].min_gain[pol] = 10000.0;
624 rad_pattern[calc_data.fstep].max_gain_idx[pol] = 0;
625 rad_pattern[calc_data.fstep].min_gain_idx[pol] = 0;
626 rad_pattern[calc_data.fstep].max_gain_tht[pol] = 0;
627 rad_pattern[calc_data.fstep].max_gain_phi[pol] = 0;
628 }
629
630 /* Signal new rad pattern data */
631 SetFlag( DRAW_NEW_RDPAT );
632
633 /* Step over theta and phi angles */
634 idx = 0;
635 for( kph = 1; kph <= fpat.nph; kph++ )
636 {
637 phi += fpat.dph;
638 pha= phi* TA;
639 thet= fpat.thets - fpat.dth;
640
641 for( kth = 1; kth <= fpat.nth; kth++ )
642 {
643 thet += fpat.dth;
644
645 if( (gnd.ksymp == 2) && (thet > 90.01) && (gnd.ifar != 1) )
646 {
647 if( rdpattern_window != NULL )
648 gtk_widget_destroy( rdpattern_window );
649 fprintf( stderr, "xnec2c: rdpat(): Theta > 90 deg with ground specified\n"
650 "Please check RP card data and correct\n" );
651 stop( _("rdpat(): Theta > 90 deg with ground specified\n"\
652 "Please check RP card data and correct"), ERR_STOP );
653 }
654
655 tha= thet* TA;
656 if( gnd.ifar != 1)
657 ffld( tha, pha, ð, &eph);
658 else
659 {
660 gfld( fpat.rfld/data.wlam, pha, thet/data.wlam,
661 ð, &eph, &erd, gnd.zrati, gnd.ksymp);
662 }
663
664 ethm2= creal( eth* conj( eth));
665 ethm= sqrt( ethm2);
666 etha= cang( eth);
667 ephm2= creal( eph* conj( eph));
668 ephm= sqrt( ephm2);
669 epha= cang( eph);
670
671 /* elliptical polarization calc. */
672 if( gnd.ifar != 1)
673 {
674 if( (ethm2 <= 1.0e-20) && (ephm2 <= 1.0e-20) )
675 {
676 tilta=0.0;
677 emajr2=0.0;
678 eminr2=0.0;
679 axrat=0.0;
680 isens= 0;
681 }
682 else
683 {
684 dfaz= epha- etha;
685 if( epha >= 0.0)
686 dfaz2= dfaz-360.0;
687 else
688 dfaz2= dfaz+360.0;
689
690 if( fabs(dfaz) > fabs(dfaz2) )
691 dfaz= dfaz2;
692
693 cdfaz= cos( dfaz* TA);
694 tstor1= ethm2- ephm2;
695 tstor2=2.0* ephm* ethm* cdfaz;
696 tilta=atan2( tstor2, tstor1)/2.0;
697 stilta= sin( tilta);
698 tstor1= tstor1* stilta* stilta;
699 tstor2= tstor2* stilta* cos( tilta);
700 emajr2= -tstor1+ tstor2+ ethm2;
701 eminr2= tstor1- tstor2+ ephm2;
702 if( eminr2 < 0.0) eminr2=0.0;
703
704 axrat= sqrt( eminr2/ emajr2);
705 if( axrat <= 1.0e-5)
706 isens= 1;
707 else if( dfaz <= 0.0)
708 isens= 2;
709 else
710 isens= 3;
711
712 } /* if( (ethm2 <= 1.0e-20) && (ephm2 <= 1.0e-20) ) */
713
714 gnmj= db10( gcon* emajr2);
715 gnmn= db10( gcon* eminr2);
716 gnv = db10( gcon* ethm2);
717 gnh = db10( gcon* ephm2);
718 gtot= db10( gcon* (ethm2+ ephm2) );
719
720 switch( fpat.inor )
721 {
722 case 0:
723 tstor1= gtot;
724 break;
725
726 case 1:
727 tstor1= gnmj;
728 break;
729
730 case 2:
731 tstor1= gnmn;
732 break;
733
734 case 3:
735 tstor1= gnv;
736 break;
737
738 case 4:
739 tstor1= gnh;
740 break;
741
742 case 5:
743 tstor1= gtot;
744 }
745
746 /* Save rad pattern gains */
747 rad_pattern[calc_data.fstep].gtot[idx] = tstor1;
748
749 /* Save axial ratio, tilt and pol sense */
750 if( isens == 2 )
751 rad_pattern[calc_data.fstep].axrt[idx] = -axrat;
752 else
753 rad_pattern[calc_data.fstep].axrt[idx] = axrat;
754 rad_pattern[calc_data.fstep].tilt[idx] = tilta;
755 rad_pattern[calc_data.fstep].sens[idx] = isens;
756
757 /* Find and save max value of gain and direction */
758 for( pol = 0; pol < NUM_POL; pol++ )
759 {
760 gain = rad_pattern[calc_data.fstep].gtot[idx] +
761 Polarization_Factor( pol, calc_data.fstep, idx);
762 if( gain < -999.99 ) gain = -999.99;
763
764 /* Find and save max value of gain and direction */
765 if( rad_pattern[calc_data.fstep].max_gain[pol] < gain )
766 {
767 rad_pattern[calc_data.fstep].max_gain[pol] = gain;
768 rad_pattern[calc_data.fstep].max_gain_tht[pol] = thet;
769 rad_pattern[calc_data.fstep].max_gain_phi[pol] = phi;
770 rad_pattern[calc_data.fstep].max_gain_idx[pol] = idx;
771 }
772
773 /* Find and save min value of gain and buffer idx */
774 if( rad_pattern[calc_data.fstep].min_gain[pol] > gain )
775 {
776 rad_pattern[calc_data.fstep].min_gain[pol] = gain;
777 rad_pattern[calc_data.fstep].min_gain_idx[pol] = idx;
778 }
779
780 } /* for( pol = 0; pol < NUM_POL; pol++ ) */
781
782 idx++;
783 continue;
784 } /* if( gnd.ifar != 1) */
785
786 } /* for( kth = 1; kth <= fpat.nth; kth++ ) */
787 } /* for( kph = 1; kph <= fpat.nph; kph++ ) */
788
789 return;
790
791 } /* void rdpat() */
792
793 /*-----------------------------------------------------------------------*/
794
795