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 "fields.h"
45 #include "shared.h"
46
47 /* common /tmi/ */
48 static tmi_t tmi;
49
50 /*common /tmh/ */
51 static tmh_t tmh;
52
53 /*-------------------------------------------------------------------*/
54
55 /* compute near e fields of a segment with sine, cosine, and */
56 /* constant currents. ground effect included. */
efld(double xi,double yi,double zi,double ai,int ij)57 void efld( double xi, double yi,
58 double zi, double ai, int ij )
59 {
60 #define txk egnd[0]
61 #define tyk egnd[1]
62 #define tzk egnd[2]
63 #define txs egnd[3]
64 #define tys egnd[4]
65 #define tzs egnd[5]
66 #define txc egnd[6]
67 #define tyc egnd[7]
68 #define tzc egnd[8]
69
70 int ip, ijx;
71 double xij, yij, rfl, salpr, zij, zp, rhox;
72 double rhoy, rhoz, rh, r, rmag, cth, px, py;
73 double xymag, xspec, yspec, rhospc, dmin;
74 complex double epx, epy, refs, refps, zrsin, zratx, zscrn;
75 complex double tezs, ters, tezc=0.0, terc=0.0, tezk=0.0, terk=0.0;
76 static complex double *egnd = NULL;
77
78 size_t mreq = 9 * sizeof(complex double);
79 mem_alloc( (void **)&egnd, mreq, "in fields.c");
80
81 xij= xi- dataj.xj;
82 yij= yi- dataj.yj;
83 ijx= ij;
84 rfl=-1.0;
85
86 for( ip = 0; ip < gnd.ksymp; ip++ )
87 {
88 if( ip == 1)
89 ijx=1;
90 rfl= -rfl;
91 salpr= dataj.salpj* rfl;
92 zij= zi- rfl* dataj.zj;
93 zp= xij* dataj.cabj+ yij* dataj.sabj+ zij* salpr;
94 rhox= xij- dataj.cabj* zp;
95 rhoy= yij- dataj.sabj* zp;
96 rhoz= zij- salpr* zp;
97
98 rh= sqrt( rhox* rhox+ rhoy* rhoy+ rhoz* rhoz+ ai* ai);
99 if( rh <= 1.0e-10)
100 {
101 rhox=0.0;
102 rhoy=0.0;
103 rhoz=0.0;
104 }
105 else
106 {
107 rhox= rhox/ rh;
108 rhoy= rhoy/ rh;
109 rhoz= rhoz/ rh;
110 }
111
112 /* lumped current element approx. for large separations */
113 r= sqrt( zp* zp+ rh* rh);
114 if( r >= dataj.rkh)
115 {
116 rmag= TP* r;
117 cth= zp/ r;
118 px= rh/ r;
119 txk= cmplx( cos( rmag),- sin( rmag));
120 py= TP* r* r;
121 tyk= ETA* cth* txk* cmplx(1.0,-1.0/ rmag)/ py;
122 tzk= ETA* px* txk* cmplx(1.0, rmag-1.0/ rmag)/(2.0* py);
123 tezk= tyk* cth- tzk* px;
124 terk= tyk* px+ tzk* cth;
125 rmag= sin( PI* dataj.s)/ PI;
126 tezc= tezk* rmag;
127 terc= terk* rmag;
128 tezk= tezk* dataj.s;
129 terk= terk* dataj.s;
130 txs=CPLX_00;
131 tys=CPLX_00;
132 tzs=CPLX_00;
133
134 } /* if( r >= dataj.rkh) */
135
136 if( r < dataj.rkh)
137 {
138 /* eksc for thin wire approx. or ekscx for extended t.w. approx. */
139 if( dataj.iexk != 1)
140 eksc( dataj.s, zp, rh, TP, ijx, &tezs, &ters,
141 &tezc, &terc, &tezk, &terk );
142 else
143 ekscx( dataj.b, dataj.s, zp, rh, TP, ijx, dataj.ind1,
144 dataj.ind2, &tezs, &ters, &tezc, &terc, &tezk, &terk);
145
146 txs= tezs* dataj.cabj+ ters* rhox;
147 tys= tezs* dataj.sabj+ ters* rhoy;
148 tzs= tezs* salpr+ ters* rhoz;
149
150 } /* if( r < dataj.rkh) */
151
152 txk= tezk* dataj.cabj+ terk* rhox;
153 tyk= tezk* dataj.sabj+ terk* rhoy;
154 tzk= tezk* salpr+ terk* rhoz;
155 txc= tezc* dataj.cabj+ terc* rhox;
156 tyc= tezc* dataj.sabj+ terc* rhoy;
157 tzc= tezc* salpr+ terc* rhoz;
158
159 if( ip == 1)
160 {
161 if( gnd.iperf <= 0)
162 {
163 zratx= gnd.zrati;
164 rmag= r;
165 xymag= sqrt( xij* xij+ yij* yij);
166
167 /* set parameters for radial wire ground screen. */
168 if( gnd.nradl > 0)
169 {
170 xspec=( xi* dataj.zj+ zi* dataj.xj)/( zi+ dataj.zj);
171 yspec=( yi* dataj.zj+ zi* dataj.yj)/( zi+ dataj.zj);
172 rhospc= sqrt( xspec* xspec+ yspec* yspec+ gnd.t2* gnd.t2);
173
174 if( rhospc <= gnd.scrwl)
175 {
176 zscrn= gnd.t1* rhospc* log( rhospc/ gnd.t2);
177 zratx=( zscrn* gnd.zrati)/( ETA* gnd.zrati+ zscrn);
178 }
179 } /* if( gnd.nradl > 0) */
180 else
181 {
182 xspec= 0.0;
183 yspec= 0.0;
184 rhospc= 0.0;
185 zscrn= 0.0;
186 zratx= 0.0;
187 }
188
189 /* calculation of reflection coefficients
190 * when ground is specified. */
191 if( xymag <= 1.0e-6)
192 {
193 px=0.0;
194 py=0.0;
195 cth=1.0;
196 zrsin=CPLX_10;
197 }
198 else
199 {
200 px= -yij/ xymag;
201 py= xij/ xymag;
202 cth= zij/ rmag;
203 zrsin= csqrt(1.0 - zratx*zratx*(1.0 - cth*cth) );
204
205 } /* if( xymag <= 1.0e-6) */
206
207 refs=( cth- zratx* zrsin)/( cth+ zratx* zrsin);
208 refps=-( zratx* cth- zrsin)/( zratx* cth+ zrsin);
209 refps= refps- refs;
210 epy= px* txk+ py* tyk;
211 epx= px* epy;
212 epy= py* epy;
213 txk= refs* txk+ refps* epx;
214 tyk= refs* tyk+ refps* epy;
215 tzk= refs* tzk;
216 epy= px* txs+ py* tys;
217 epx= px* epy;
218 epy= py* epy;
219 txs= refs* txs+ refps* epx;
220 tys= refs* tys+ refps* epy;
221 tzs= refs* tzs;
222 epy= px* txc+ py* tyc;
223 epx= px* epy;
224 epy= py* epy;
225 txc= refs* txc+ refps* epx;
226 tyc= refs* tyc+ refps* epy;
227 tzc= refs* tzc;
228
229 } /* if( gnd.iperf <= 0) */
230
231 dataj.exk= dataj.exk- txk* gnd.frati;
232 dataj.eyk= dataj.eyk- tyk* gnd.frati;
233 dataj.ezk= dataj.ezk- tzk* gnd.frati;
234 dataj.exs= dataj.exs- txs* gnd.frati;
235 dataj.eys= dataj.eys- tys* gnd.frati;
236 dataj.ezs= dataj.ezs- tzs* gnd.frati;
237 dataj.exc= dataj.exc- txc* gnd.frati;
238 dataj.eyc= dataj.eyc- tyc* gnd.frati;
239 dataj.ezc= dataj.ezc- tzc* gnd.frati;
240 continue;
241
242 } /* if( ip == 1) */
243
244 dataj.exk= txk;
245 dataj.eyk= tyk;
246 dataj.ezk= tzk;
247 dataj.exs= txs;
248 dataj.eys= tys;
249 dataj.ezs= tzs;
250 dataj.exc= txc;
251 dataj.eyc= tyc;
252 dataj.ezc= tzc;
253
254 } /* for( ip = 0; ip < gnd.ksymp; ip++ ) */
255
256 if( gnd.iperf != 2)
257 return;
258
259 /* field due to ground using sommerfeld/norton */
260 incom.sn= sqrt( dataj.cabj* dataj.cabj+ dataj.sabj* dataj.sabj);
261 if( incom.sn >= 1.0e-5)
262 {
263 incom.xsn= dataj.cabj/ incom.sn;
264 incom.ysn= dataj.sabj/ incom.sn;
265 }
266 else
267 {
268 incom.sn=0.0;
269 incom.xsn=1.0;
270 incom.ysn=0.0;
271 }
272
273 /* displace observation point for thin wire approximation */
274 zij= zi+ dataj.zj;
275 salpr= -dataj.salpj;
276 rhox= dataj.sabj* zij- salpr* yij;
277 rhoy= salpr* xij- dataj.cabj* zij;
278 rhoz= dataj.cabj* yij- dataj.sabj* xij;
279 rh= rhox* rhox+ rhoy* rhoy+ rhoz* rhoz;
280
281 if( rh <= 1.0e-10)
282 {
283 incom.xo= xi- ai* incom.ysn;
284 incom.yo= yi+ ai* incom.xsn;
285 incom.zo= zi;
286 }
287 else
288 {
289 rh= ai/ sqrt( rh);
290 if( rhoz < 0.0)
291 rh= -rh;
292 incom.xo= xi+ rh* rhox;
293 incom.yo= yi+ rh* rhoy;
294 incom.zo= zi+ rh* rhoz;
295
296 } /* if( rh <= 1.0e-10) */
297
298 r= xij* xij+ yij* yij+ zij* zij;
299 if( r <= .95)
300 {
301 double shaf;
302
303 /* field from interpolation is integrated over segment */
304 incom.isnor=1;
305 dmin=
306 creal(dataj.exk*conj(dataj.exk)) +
307 creal(dataj.eyk*conj(dataj.eyk)) +
308 creal(dataj.ezk*conj(dataj.ezk));
309 dmin=.01* sqrt( dmin);
310 shaf=.5* dataj.s;
311 rom2(- shaf, shaf, egnd, dmin);
312 }
313 else
314 {
315 /* norton field equations and lumped
316 * current element approximation */
317 incom.isnor=2;
318 sflds(0.0, egnd);
319 } /* if( r <= .95) */
320
321 if( r > .95)
322 {
323 zp= xij* dataj.cabj+ yij* dataj.sabj+ zij* salpr;
324 rh= r- zp* zp;
325 if( rh <= 1.0e-10)
326 dmin=0.0;
327 else
328 dmin= sqrt( rh/( rh+ ai* ai));
329
330 if( dmin <= .95)
331 {
332 px=1.0- dmin;
333 terk=( txk* dataj.cabj+ tyk* dataj.sabj+ tzk* salpr)* px;
334 txk= dmin* txk+ terk* dataj.cabj;
335 tyk= dmin* tyk+ terk* dataj.sabj;
336 tzk= dmin* tzk+ terk* salpr;
337 ters=( txs* dataj.cabj+ tys* dataj.sabj+ tzs* salpr)* px;
338 txs= dmin* txs+ ters* dataj.cabj;
339 tys= dmin* tys+ ters* dataj.sabj;
340 tzs= dmin* tzs+ ters* salpr;
341 terc=( txc* dataj.cabj+ tyc* dataj.sabj+ tzc* salpr)* px;
342 txc= dmin* txc+ terc* dataj.cabj;
343 tyc= dmin* tyc+ terc* dataj.sabj;
344 tzc= dmin* tzc+ terc* salpr;
345
346 } /* if( dmin <= .95) */
347
348 } /* if( r > .95) */
349
350 dataj.exk= dataj.exk+ txk;
351 dataj.eyk= dataj.eyk+ tyk;
352 dataj.ezk= dataj.ezk+ tzk;
353 dataj.exs= dataj.exs+ txs;
354 dataj.eys= dataj.eys+ tys;
355 dataj.ezs= dataj.ezs+ tzs;
356 dataj.exc= dataj.exc+ txc;
357 dataj.eyc= dataj.eyc+ tyc;
358 dataj.ezc= dataj.ezc+ tzc;
359
360 return;
361 }
362
363 /*-----------------------------------------------------------------------*/
364
365 /* compute e field of sine, cosine, and constant */
366 /* current filaments by thin wire approximation. */
eksc(double s,double z,double rh,double xk,int ij,complex double * ezs,complex double * ers,complex double * ezc,complex double * erc,complex double * ezk,complex double * erk)367 void eksc( double s, double z, double rh,
368 double xk, int ij, complex double *ezs,
369 complex double *ers, complex double *ezc,
370 complex double *erc, complex double *ezk,
371 complex double *erk )
372 {
373 double rhk, sh, shk, ss, cs, z1a, z2a, cint, sint;
374 complex double gz1, gz2, gp1, gp2, gzp1, gzp2;
375
376 tmi.ij= ij;
377 tmi.zpk= xk* z;
378 rhk= xk* rh;
379 tmi.rkb2= rhk* rhk;
380 sh=.5* s;
381 shk= xk* sh;
382 ss= sin( shk);
383 cs= cos( shk);
384 z2a= sh- z;
385 z1a=-( sh+ z);
386 gx( z1a, rh, xk, &gz1, &gp1);
387 gx( z2a, rh, xk, &gz2, &gp2);
388 gzp1= gp1* z1a;
389 gzp2= gp2* z2a;
390 *ezs= CONST1*(( gz2- gz1)* cs* xk-( gzp2+ gzp1)* ss);
391 *ezc= -CONST1*(( gz2+ gz1)* ss* xk+( gzp2- gzp1)* cs);
392 *erk= CONST1*( gp2- gp1)* rh;
393 intx(- shk, shk, rhk, ij, &cint, &sint);
394 *ezk= -CONST1*( gzp2- gzp1+ xk* xk* cmplx( cint,- sint));
395 gzp1= gzp1* z1a;
396 gzp2= gzp2* z2a;
397
398 if( rh >= 1.0e-10)
399 {
400 *ers= -CONST1*(( gzp2+ gzp1+ gz2+ gz1)*
401 ss-( z2a* gz2- z1a* gz1)* cs*xk)/ rh;
402 *erc= -CONST1*(( gzp2- gzp1+ gz2- gz1)*
403 cs+( z2a* gz2+ z1a* gz1)* ss*xk)/ rh;
404 return;
405 }
406
407 *ers = CPLX_00;
408 *erc = CPLX_00;
409
410 return;
411 }
412
413 /*-----------------------------------------------------------------------*/
414
415 /* compute e field of sine, cosine, and constant current */
416 /* filaments by extended thin wire approximation. */
ekscx(double bx,double s,double z,double rhx,double xk,int ij,int inx1,int inx2,complex double * ezs,complex double * ers,complex double * ezc,complex double * erc,complex double * ezk,complex double * erk)417 void ekscx( double bx, double s, double z,
418 double rhx, double xk, int ij, int inx1,
419 int inx2, complex double *ezs, complex double *ers,
420 complex double *ezc, complex double *erc,
421 complex double *ezk, complex double *erk )
422 {
423 int ira;
424 double b, rh, sh, rhk, shk, ss, cs, z1a;
425 double z2a, a2, bk, bk2, cint, sint;
426 complex double gz1, gz2, gzp1, gzp2, gr1, gr2;
427 complex double grp1, grp2, grk1, grk2, gzz1, gzz2;
428
429 if( rhx >= bx)
430 {
431 rh= rhx;
432 b= bx;
433 ira=0;
434 }
435 else
436 {
437 rh= bx;
438 b= rhx;
439 ira=1;
440 }
441
442 sh=.5* s;
443 tmi.ij= ij;
444 tmi.zpk= xk* z;
445 rhk= xk* rh;
446 tmi.rkb2= rhk* rhk;
447 shk= xk* sh;
448 ss= sin( shk);
449 cs= cos( shk);
450 z2a= sh- z;
451 z1a=-( sh+ z);
452 a2= b* b;
453
454 if( inx1 != 2)
455 gxx( z1a, rh, b, a2, xk, ira, &gz1,
456 &gzp1, &gr1, &grp1, &grk1, &gzz1);
457 else
458 {
459 gx( z1a, rhx, xk, &gz1, &grk1);
460 gzp1= grk1* z1a;
461 gr1= gz1/ rhx;
462 grp1= gzp1/ rhx;
463 grk1= grk1* rhx;
464 gzz1= CPLX_00;
465 }
466
467 if( inx2 != 2)
468 gxx( z2a, rh, b, a2, xk, ira, &gz2,
469 &gzp2, &gr2, &grp2, &grk2, &gzz2);
470 else
471 {
472 gx( z2a, rhx, xk, &gz2, &grk2);
473 gzp2= grk2* z2a;
474 gr2= gz2/ rhx;
475 grp2= gzp2/ rhx;
476 grk2= grk2* rhx;
477 gzz2= CPLX_00;
478 }
479
480 *ezs= CONST1*(( gz2- gz1)* cs* xk-( gzp2+ gzp1)* ss);
481 *ezc= -CONST1*(( gz2+ gz1)* ss* xk+( gzp2- gzp1)* cs);
482 *ers= -CONST1*(( z2a* grp2+ z1a* grp1+ gr2+ gr1)*ss
483 -( z2a* gr2- z1a* gr1)* cs* xk);
484 *erc= -CONST1*(( z2a* grp2- z1a* grp1+ gr2- gr1)*cs
485 +( z2a* gr2+ z1a* gr1)* ss* xk);
486 *erk= CONST1*( grk2- grk1);
487 intx(- shk, shk, rhk, ij, &cint, &sint);
488 bk= b* xk;
489 bk2= bk* bk*.25;
490 *ezk= -CONST1*( gzp2- gzp1+ xk* xk*(1.0- bk2)*
491 cmplx( cint,- sint)-bk2*( gzz2- gzz1));
492
493 return;
494 }
495
496 /*-----------------------------------------------------------------------*/
497
498 /* gf computes the integrand exp(jkr)/(kr) for numerical integration. */
gf(double zk,double * co,double * si)499 void gf( double zk, double *co, double *si )
500 {
501 double zdk, rk, rks;
502
503 zdk= zk- tmi.zpk;
504 rk= sqrt( tmi.rkb2+ zdk* zdk);
505 *si= sin( rk)/ rk;
506
507 if( tmi.ij != 0 )
508 {
509 *co= cos( rk)/ rk;
510 return;
511 }
512
513 if( rk >= .2l)
514 {
515 *co=( cos( rk)-1.0)/ rk;
516 return;
517 }
518
519 rks= rk* rk;
520 *co=((-1.38888889e-3* rks+4.16666667e-2)* rks-.5)* rk;
521
522 return;
523 }
524
525 /*-----------------------------------------------------------------------*/
526
527 /* integrand for h field of a wire */
gh(double zk,double * hr,double * hi)528 void gh( double zk, double *hr, double *hi)
529 {
530 double rs, r, ckr, skr, rr2, rr3;
531
532 rs= zk- tmh.zpka;
533 rs= tmh.rhks+ rs* rs;
534 r= sqrt( rs);
535 ckr= cos( r);
536 skr= sin( r);
537 rr2=1.0/ rs;
538 rr3= rr2/ r;
539 *hr= skr* rr2+ ckr* rr3;
540 *hi= ckr* rr2- skr* rr3;
541
542 return;
543 }
544
545 /*-----------------------------------------------------------------------*/
546
547 /* gwave computes the electric field, including ground wave, of a */
548 /* current element over a ground plane using formulas of k.a. norton */
549 /* (proc. ire, sept., 1937, pp.1203,1236) */
550
gwave(complex double * erv,complex double * ezv,complex double * erh,complex double * ezh,complex double * eph)551 void gwave( complex double *erv, complex double *ezv,
552 complex double *erh, complex double *ezh,
553 complex double *eph )
554 {
555 double sppp, sppp2, cppp2, cppp, spp, spp2, cpp2, cpp;
556 complex double rk1, rk2, t1, t2, t3, t4, p1, rv;
557 complex double omr, w, f, q1, rh, v, g, xr1, xr2;
558 complex double x1, x2, x3, x4, x5, x6, x7;
559
560 sppp= gwav.zmh/ gwav.r1;
561 sppp2= sppp* sppp;
562 cppp2=1.0- sppp2;
563
564 if( cppp2 < 1.0e-20)
565 cppp2=1.0e-20;
566
567 cppp= sqrt( cppp2);
568 spp= gwav.zph/ gwav.r2;
569 spp2= spp* spp;
570 cpp2=1.0- spp2;
571
572 if( cpp2 < 1.0e-20)
573 cpp2=1.0e-20;
574
575 cpp= sqrt( cpp2);
576 rk1= -TPJ* gwav.r1;
577 rk2= -TPJ* gwav.r2;
578 t1=1.0 -gwav.u2* cpp2;
579 t2= csqrt( t1);
580 t3=(1.0 -1.0/ rk1)/ rk1;
581 t4=(1.0 -1.0/ rk2)/ rk2;
582 p1= rk2* gwav.u2* t1/(2.0* cpp2);
583 rv=( spp- gwav.u* t2)/( spp+ gwav.u* t2);
584 omr=1.0- rv;
585 w=1.0/ omr;
586 w=(4.0 + I*0.0)* p1* w* w;
587 fbar( w, &f );
588 q1= rk2* t1/(2.0* gwav.u2* cpp2);
589 rh=( t2- gwav.u* spp)/( t2+ gwav.u* spp);
590 v=1.0/(1.0+ rh);
591 v=(4.0 + I*0.0)* q1* v* v;
592 fbar( v, &g );
593 xr1= gwav.xx1/ gwav.r1;
594 xr2= gwav.xx2/ gwav.r2;
595 x1= cppp2* xr1;
596 x2= rv* cpp2* xr2;
597 x3= omr* cpp2* f* xr2;
598 x4= gwav.u* t2* spp*2.0* xr2/ rk2;
599 x5= xr1* t3*(1.0-3.0* sppp2);
600 x6= xr2* t4*(1.0-3.0* spp2);
601 *ezv=( x1+ x2+ x3- x4- x5- x6)* (-CONST4);
602 x1= sppp* cppp* xr1;
603 x2= rv* spp* cpp* xr2;
604 x3= cpp* omr* gwav.u* t2* f* xr2;
605 x4= spp* cpp* omr* xr2/ rk2;
606 x5=3.0* sppp* cppp* t3* xr1;
607 x6= cpp* gwav.u* t2* omr* xr2/ rk2*.50;
608 x7=3.0* spp* cpp* t4* xr2;
609 *erv=-( x1+ x2- x3+ x4- x5+ x6- x7)* (-CONST4);
610 *ezh=-( x1- x2+ x3- x4- x5- x6+ x7)* (-CONST4);
611 x1= sppp2* xr1;
612 x2= rv* spp2* xr2;
613 x4= gwav.u2* t1* omr* f* xr2;
614 x5= t3*(1.0-3.0* cppp2)* xr1;
615 x6= t4*(1.0-3.0* cpp2)*(1.0- gwav.u2*(1.0+ rv )-
616 gwav.u2* omr* f)* xr2;
617 x7= gwav.u2* cpp2* omr*(1.0-1.0/ rk2)*( f*( gwav.u2* t1 -
618 spp2-1.0/ rk2)+1.0/rk2)* xr2;
619 *erh=( x1- x2- x4- x5+ x6+ x7)* (-CONST4);
620 x1= xr1;
621 x2= rh* xr2;
622 x3=( rh+1.0)* g* xr2;
623 x4= t3* xr1;
624 x5= t4*(1.0- gwav.u2*(1.0+ rv)- gwav.u2* omr* f)* xr2;
625 x6=.5* gwav.u2* omr*( f*( gwav.u2* t1- spp2-1.0/ rk2) +
626 1.0/ rk2)* xr2/ rk2;
627 *eph=-( x1- x2+ x3- x4+ x5+ x6)* (-CONST4);
628
629 return;
630 }
631
632 /*-----------------------------------------------------------------------*/
633
634 /* segment end contributions for thin wire approx. */
gx(double zz,double rh,double xk,complex double * gz,complex double * gzp)635 void gx( double zz, double rh, double xk,
636 complex double *gz, complex double *gzp)
637 {
638 double r, r2, rkz;
639
640 r2= zz* zz+ rh* rh;
641 r= sqrt( r2);
642 rkz= xk* r;
643 *gz= cmplx( cos( rkz),- sin( rkz))/ r;
644 *gzp= -cmplx(1.0, rkz)* *gz/ r2;
645
646 return;
647 }
648
649 /*-----------------------------------------------------------------------*/
650
651 /* segment end contributions for ext. thin wire approx. */
gxx(double zz,double rh,double a,double a2,double xk,int ira,complex double * g1,complex double * g1p,complex double * g2,complex double * g2p,complex double * g3,complex double * gzp)652 void gxx( double zz, double rh, double a,
653 double a2, double xk, int ira, complex double *g1,
654 complex double *g1p, complex double *g2,
655 complex double *g2p, complex double *g3,
656 complex double *gzp )
657 {
658 double r, r2, r4, rk, rk2, rh2, t1, t2;
659 complex double gz, c1, c2, c3;
660
661 r2= zz* zz+ rh* rh;
662 r= sqrt( r2);
663 r4= r2* r2;
664 rk= xk* r;
665 rk2= rk* rk;
666 rh2= rh* rh;
667 t1=.25* a2* rh2/ r4;
668 t2=.5* a2/ r2;
669 c1= cmplx(1.0, rk);
670 c2=3.0* c1- rk2;
671 c3= cmplx(6.0, rk)* rk2-15.0* c1;
672 gz= cmplx( cos( rk),- sin( rk))/ r;
673 *g2= gz*(1.0+ t1* c2);
674 *g1= *g2- t2* c1* gz;
675 gz= gz/ r2;
676 *g2p= gz*( t1* c3- c1);
677 *gzp= t2* c2* gz;
678 *g3= *g2p+ *gzp;
679 *g1p= *g3* zz;
680
681 if( ira != 1)
682 {
683 *g3=( *g3+ *gzp)* rh;
684 *gzp= -zz* c1* gz;
685
686 if( rh <= 1.0e-10)
687 {
688 *g2=0.0;
689 *g2p=0.0;
690 return;
691 }
692
693 *g2= *g2/ rh;
694 *g2p= *g2p* zz/ rh;
695 return;
696
697 } /* if( ira != 1) */
698
699 t2=.5* a;
700 *g2= -t2* c1* gz;
701 *g2p= t2* gz* c2/ r2;
702 *g3= rh2* *g2p- a* gz* c1;
703 *g2p= *g2p* zz;
704 *gzp= -zz* c1* gz;
705
706 return;
707 }
708
709 /*-----------------------------------------------------------------------*/
710
711 /* hfk computes the h field of a uniform current */
712 /* filament by numerical integration */
hfk(double el1,double el2,double rhk,double zpkx,double * sgr,double * sgi)713 void hfk( double el1, double el2, double rhk,
714 double zpkx, double *sgr, double *sgi )
715 {
716 int nx = 1, nma = 65536, nts = 4;
717 int ns, nt;
718 int flag = TRUE;
719 double rx = 1.0e-4;
720 double z, ze, s, ep, zend, dz=0.0, zp, dzot=0.0;
721 double t00r, g1r, g5r=0.0, t00i, g1i, g5i=0.0, t01r;
722 double g3r=0.0, t01i, g3i=0.0, t10r, t10i, te1i;
723 double te1r, g2r, g4r, t02i, g2i, g4i, t11r, t11i;
724 double t20r, t20i, te2i, te2r, t02r;
725
726 tmh.zpka= zpkx;
727 tmh.rhks= rhk* rhk;
728 z= el1;
729 ze= el2;
730 s= ze- z;
731 ep= s/(10.0* nma);
732 zend= ze- ep;
733 *sgr=0.0;
734 *sgi=0.0;
735 ns= nx;
736 nt=0;
737 gh( z, &g1r, &g1i);
738
739 while( TRUE )
740 {
741 if( flag )
742 {
743 dz= s/ ns;
744 zp= z+ dz;
745
746 if( zp > ze )
747 {
748 dz= ze- z;
749 if( fabs(dz) <= ep )
750 {
751 *sgr= *sgr* rhk*.5;
752 *sgi= *sgi* rhk*.5;
753 return;
754 }
755 }
756
757 dzot= dz*.5;
758 zp= z+ dzot;
759 gh( zp, &g3r, &g3i);
760 zp= z+ dz;
761 gh( zp, &g5r, &g5i);
762
763 } /* if( flag ) */
764
765 t00r=( g1r+ g5r)* dzot;
766 t00i=( g1i+ g5i)* dzot;
767 t01r=( t00r+ dz* g3r)*0.5;
768 t01i=( t00i+ dz* g3i)*0.5;
769 t10r=(4.0* t01r- t00r)/3.0;
770 t10i=(4.0* t01i- t00i)/3.0;
771
772 test( t01r, t10r, &te1r, t01i, t10i, &te1i, 0.0);
773 if( (te1i <= rx) && (te1r <= rx) )
774 {
775 *sgr= *sgr+ t10r;
776 *sgi= *sgi+ t10i;
777 nt += 2;
778
779 z += dz;
780 if( z >= zend)
781 {
782 *sgr= *sgr* rhk*.5;
783 *sgi= *sgi* rhk*.5;
784 return;
785 }
786
787 g1r= g5r;
788 g1i= g5i;
789 if( nt >= nts)
790 if( ns > nx)
791 {
792 ns= ns/2;
793 nt=1;
794 }
795 flag = TRUE;
796 continue;
797
798 } /* if( (te1i <= rx) && (te1r <= rx) ) */
799
800 zp= z+ dz*0.25;
801 gh( zp, &g2r, &g2i);
802 zp= z+ dz*0.75;
803 gh( zp, &g4r, &g4i);
804 t02r=( t01r+ dzot*( g2r+ g4r))*0.5;
805 t02i=( t01i+ dzot*( g2i+ g4i))*0.5;
806 t11r=(4.0* t02r- t01r)/3.0;
807 t11i=(4.0* t02i- t01i)/3.0;
808 t20r=(16.0* t11r- t10r)/15.0;
809 t20i=(16.0* t11i- t10i)/15.0;
810
811 test( t11r, t20r, &te2r, t11i, t20i, &te2i, 0.0);
812 if( (te2i > rx) || (te2r > rx) )
813 {
814 nt=0;
815 if( ns >= nma)
816 {
817 fprintf( stderr,
818 "\nxnec2c: step size limited at z= %10.5f", z );
819 }
820 else
821 {
822 ns= ns*2;
823 dz= s/ ns;
824 dzot= dz*0.5;
825 g5r= g3r;
826 g5i= g3i;
827 g3r= g2r;
828 g3i= g2i;
829
830 flag = FALSE;
831 continue;
832 }
833
834 } /* if( (te2i > rx) || (te2r > rx) ) */
835
836 *sgr= *sgr+ t20r;
837 *sgi= *sgi+ t20i;
838 nt++;
839
840 z += dz;
841 if( z >= zend)
842 {
843 *sgr= *sgr* rhk*.5;
844 *sgi= *sgi* rhk*.5;
845 return;
846 }
847
848 g1r= g5r;
849 g1i= g5i;
850 if( nt >= nts)
851 if( ns > nx)
852 {
853 ns= ns/2;
854 nt=1;
855 }
856 flag = TRUE;
857
858 } /* while( TRUE ) */
859
860 }
861
862 /*-----------------------------------------------------------------------*/
863
864 /* hintg computes the h field of a patch current */
hintg(double xi,double yi,double zi)865 void hintg( double xi, double yi, double zi )
866 {
867 int ip;
868 double rx, ry, rfl, xymag, pxx, pyy, cth;
869 double rz, rsq, r, rk, cr, sr, t1zr, t2zr;
870 complex double gam, f1x, f1y, f1z, f2x, f2y, f2z, rrv, rrh;
871
872 rx= xi- dataj.xj;
873 ry= yi- dataj.yj;
874 rfl=-1.0;
875 dataj.exk=CPLX_00;
876 dataj.eyk=CPLX_00;
877 dataj.ezk=CPLX_00;
878 dataj.exs=CPLX_00;
879 dataj.eys=CPLX_00;
880 dataj.ezs=CPLX_00;
881
882 for( ip = 1; ip <= gnd.ksymp; ip++ )
883 {
884 rfl= -rfl;
885 rz= zi- dataj.zj* rfl;
886 rsq= rx* rx+ ry* ry+ rz* rz;
887
888 if( rsq < 1.0e-20)
889 continue;
890
891 r = sqrt( rsq );
892 rk= TP* r;
893 cr= cos( rk);
894 sr= sin( rk);
895 gam=-( cmplx(cr,-sr)+rk*cmplx(sr,cr) )/( FPI*rsq*r )* dataj.s;
896 dataj.exc= gam* rx;
897 dataj.eyc= gam* ry;
898 dataj.ezc= gam* rz;
899 t1zr= dataj.t1zj* rfl;
900 t2zr= dataj.t2zj* rfl;
901 f1x= dataj.eyc* t1zr- dataj.ezc* dataj.t1yj;
902 f1y= dataj.ezc* dataj.t1xj- dataj.exc* t1zr;
903 f1z= dataj.exc* dataj.t1yj- dataj.eyc* dataj.t1xj;
904 f2x= dataj.eyc* t2zr- dataj.ezc* dataj.t2yj;
905 f2y= dataj.ezc* dataj.t2xj- dataj.exc* t2zr;
906 f2z= dataj.exc* dataj.t2yj- dataj.eyc* dataj.t2xj;
907
908 if( ip != 1)
909 {
910 if( gnd.iperf == 1)
911 {
912 f1x= -f1x;
913 f1y= -f1y;
914 f1z= -f1z;
915 f2x= -f2x;
916 f2y= -f2y;
917 f2z= -f2z;
918 }
919 else
920 {
921 xymag= sqrt( rx* rx+ ry* ry);
922 if( xymag <= 1.0e-6)
923 {
924 pxx=0.0;
925 pyy=0.0;
926 cth=1.0;
927 rrv=CPLX_10;
928 }
929 else
930 {
931 pxx= -ry/ xymag;
932 pyy= rx/ xymag;
933 cth= rz/ r;
934 rrv= csqrt(1.0- gnd.zrati* gnd.zrati*(1.0- cth* cth));
935
936 } /* if( xymag <= 1.0e-6) */
937
938 rrh= gnd.zrati* cth;
939 rrh=( rrh- rrv)/( rrh+ rrv);
940 rrv= gnd.zrati* rrv;
941 rrv=-( cth- rrv)/( cth+ rrv);
942 gam=( f1x* pxx+ f1y* pyy)*( rrv- rrh);
943 f1x= f1x* rrh+ gam* pxx;
944 f1y= f1y* rrh+ gam* pyy;
945 f1z= f1z* rrh;
946 gam=( f2x* pxx+ f2y* pyy)*( rrv- rrh);
947 f2x= f2x* rrh+ gam* pxx;
948 f2y= f2y* rrh+ gam* pyy;
949 f2z= f2z* rrh;
950
951 } /* if( gnd.iperf == 1) */
952
953 } /* if( ip != 1) */
954
955 dataj.exk += f1x;
956 dataj.eyk += f1y;
957 dataj.ezk += f1z;
958 dataj.exs += f2x;
959 dataj.eys += f2y;
960 dataj.ezs += f2z;
961
962 } /* for( ip = 1; ip <= gnd.ksymp; ip++ ) */
963
964 return;
965 }
966
967 /*-----------------------------------------------------------------------*/
968
969 /* hsfld computes the h field for constant, sine, and */
970 /* cosine current on a segment including ground effects. */
hsfld(double xi,double yi,double zi,double ai)971 void hsfld( double xi, double yi,
972 double zi, double ai )
973 {
974 int ip;
975 double xij, yij, rfl, salpr, zij, zp, rhox;
976 double rhoy, rhoz, rh, phx, phy, phz, rmag;
977 double xymag, xspec, yspec, rhospc, px, py, cth;
978 complex double hpk, hps, hpc, qx, qy, qz, rrv, rrh, zratx;
979
980 xij= xi- dataj.xj;
981 yij= yi- dataj.yj;
982 rfl=-1.0;
983
984 for( ip = 0; ip < gnd.ksymp; ip++ )
985 {
986 rfl= -rfl;
987 salpr= dataj.salpj* rfl;
988 zij= zi- rfl* dataj.zj;
989 zp= xij* dataj.cabj+ yij* dataj.sabj+ zij* salpr;
990 rhox= xij- dataj.cabj* zp;
991 rhoy= yij- dataj.sabj* zp;
992 rhoz= zij- salpr* zp;
993 rh= sqrt( rhox* rhox+ rhoy* rhoy+ rhoz* rhoz+ ai* ai);
994
995 if( rh <= 1.0e-10)
996 {
997 dataj.exk=0.0;
998 dataj.eyk=0.0;
999 dataj.ezk=0.0;
1000 dataj.exs=0.0;
1001 dataj.eys=0.0;
1002 dataj.ezs=0.0;
1003 dataj.exc=0.0;
1004 dataj.eyc=0.0;
1005 dataj.ezc=0.0;
1006 continue;
1007 }
1008
1009 rhox= rhox/ rh;
1010 rhoy= rhoy/ rh;
1011 rhoz= rhoz/ rh;
1012 phx= dataj.sabj* rhoz- salpr* rhoy;
1013 phy= salpr* rhox- dataj.cabj* rhoz;
1014 phz= dataj.cabj* rhoy- dataj.sabj* rhox;
1015
1016 hsflx( dataj.s, rh, zp, &hpk, &hps, &hpc);
1017
1018 if( ip == 1 )
1019 {
1020 if( gnd.iperf != 1 )
1021 {
1022 zratx= gnd.zrati;
1023 rmag= sqrt( zp* zp+ rh* rh);
1024 xymag= sqrt( xij* xij+ yij* yij);
1025
1026 /* set parameters for radial wire ground screen. */
1027 if( gnd.nradl > 0)
1028 {
1029 xspec=( xi* dataj.zj+ zi* dataj.xj)/( zi+ dataj.zj);
1030 yspec=( yi* dataj.zj+ zi* dataj.yj)/( zi+ dataj.zj);
1031 rhospc= sqrt( xspec* xspec+ yspec* yspec+ gnd.t2* gnd.t2);
1032
1033 if( rhospc <= gnd.scrwl)
1034 {
1035 rrv= gnd.t1* rhospc* log( rhospc/ gnd.t2);
1036 zratx=( rrv* gnd.zrati)/( ETA* gnd.zrati+ rrv);
1037 }
1038 } /* if( gnd.nradl > 0) */
1039
1040 /* calculation of reflection coefficients
1041 * when ground is specified. */
1042 if( xymag <= 1.0e-6)
1043 {
1044 px=0.0;
1045 py=0.0;
1046 cth=1.0;
1047 rrv=CPLX_10;
1048 }
1049 else
1050 {
1051 px= -yij/ xymag;
1052 py= xij/ xymag;
1053 cth= zij/ rmag;
1054 rrv= csqrt(1.0- zratx* zratx*(1.0- cth* cth));
1055 }
1056
1057 rrh= zratx* cth;
1058 rrh=-( rrh- rrv)/( rrh+ rrv);
1059 rrv= zratx* rrv;
1060 rrv=( cth- rrv)/( cth+ rrv);
1061 qy=( phx* px+ phy* py)*( rrv- rrh);
1062 qx= qy* px+ phx* rrh;
1063 qy= qy* py+ phy* rrh;
1064 qz= phz* rrh;
1065 dataj.exk= dataj.exk- hpk* qx;
1066 dataj.eyk= dataj.eyk- hpk* qy;
1067 dataj.ezk= dataj.ezk- hpk* qz;
1068 dataj.exs= dataj.exs- hps* qx;
1069 dataj.eys= dataj.eys- hps* qy;
1070 dataj.ezs= dataj.ezs- hps* qz;
1071 dataj.exc= dataj.exc- hpc* qx;
1072 dataj.eyc= dataj.eyc- hpc* qy;
1073 dataj.ezc= dataj.ezc- hpc* qz;
1074 continue;
1075
1076 } /* if( gnd.iperf != 1 ) */
1077 else
1078 {
1079 xspec= 0.0;
1080 yspec= 0.0;
1081 rhospc= 0.0;
1082 rrv= 0.0;
1083 zratx= 0.0;
1084 }
1085
1086 dataj.exk= dataj.exk- hpk* phx;
1087 dataj.eyk= dataj.eyk- hpk* phy;
1088 dataj.ezk= dataj.ezk- hpk* phz;
1089 dataj.exs= dataj.exs- hps* phx;
1090 dataj.eys= dataj.eys- hps* phy;
1091 dataj.ezs= dataj.ezs- hps* phz;
1092 dataj.exc= dataj.exc- hpc* phx;
1093 dataj.eyc= dataj.eyc- hpc* phy;
1094 dataj.ezc= dataj.ezc- hpc* phz;
1095 continue;
1096
1097 } /* if( ip == 1 ) */
1098
1099 dataj.exk= hpk* phx;
1100 dataj.eyk= hpk* phy;
1101 dataj.ezk= hpk* phz;
1102 dataj.exs= hps* phx;
1103 dataj.eys= hps* phy;
1104 dataj.ezs= hps* phz;
1105 dataj.exc= hpc* phx;
1106 dataj.eyc= hpc* phy;
1107 dataj.ezc= hpc* phz;
1108
1109 } /* for( ip = 0; ip < gnd.ksymp; ip++ ) */
1110
1111 return;
1112 }
1113
1114 /*-----------------------------------------------------------------------*/
1115
1116 /* calculates h field of sine cosine, and constant current of segment */
hsflx(double s,double rh,double zpx,complex double * hpk,complex double * hps,complex double * hpc)1117 void hsflx( double s, double rh, double zpx,
1118 complex double *hpk, complex double *hps,
1119 complex double *hpc )
1120 {
1121 complex double fjk, ekr1, ekr2, t1, t2, cons;
1122
1123 fjk = -TPJ;
1124 if( rh >= 1.0e-10)
1125 {
1126 double zp, z2a, hss, dh, z1;
1127 double rhz, dk, cdk, sdk, hkr, hki;
1128
1129 if( zpx >= 0.0)
1130 {
1131 zp= zpx;
1132 hss=1.0;
1133 }
1134 else
1135 {
1136 zp= -zpx;
1137 hss=-1.0;
1138 }
1139
1140 dh=.5* s;
1141 z1= zp+ dh;
1142 z2a= zp- dh;
1143 if( z2a >= 1.0e-7)
1144 rhz= rh/ z2a;
1145 else
1146 rhz=1.0;
1147
1148 dk= TP* dh;
1149 cdk= cos( dk);
1150 sdk= sin( dk);
1151 hfk(- dk, dk, rh* TP, zp* TP, &hkr, &hki);
1152 *hpk= cmplx( hkr, hki);
1153
1154 if( rhz >= 1.0e-3)
1155 {
1156 double rh2, r1, r2;
1157
1158 rh2= rh* rh;
1159 r1= sqrt( rh2+ z1* z1);
1160 r2= sqrt( rh2+ z2a* z2a);
1161 ekr1= cexp( fjk* r1);
1162 ekr2= cexp( fjk* r2);
1163 t1= z1* ekr1/ r1;
1164 t2= z2a* ekr2/ r2;
1165 *hps=( cdk*( ekr2- ekr1)- CPLX_01* sdk*( t2+ t1))* hss;
1166 *hpc= -sdk*( ekr2+ ekr1)- CPLX_01* cdk*( t2- t1);
1167 cons= -CPLX_01/(2.0* TP* rh);
1168 *hps= cons* *hps;
1169 *hpc= cons* *hpc;
1170 return;
1171
1172 } /* if( rhz >= 1.0e-3) */
1173
1174 ekr1= cmplx( cdk, sdk)/( z2a* z2a);
1175 ekr2= cmplx( cdk,- sdk)/( z1* z1);
1176 t1= TP*(1.0/ z1-1.0/ z2a);
1177 t2= cexp( fjk* zp)* rh/PI8;
1178 *hps= t2*( t1+( ekr1+ ekr2)* sdk)* hss;
1179 *hpc= t2*(- CPLX_01* t1+( ekr1- ekr2)* cdk);
1180 return;
1181
1182 } /* if( rh >= 1.0e-10) */
1183
1184 *hps=CPLX_00;
1185 *hpc=CPLX_00;
1186 *hpk=CPLX_00;
1187
1188 return;
1189 }
1190
1191 /*-----------------------------------------------------------------------*/
1192
1193 /* nefld computes the near field at specified points in space after */
1194 /* the structure currents have been computed. */
nefld(double xob,double yob,double zob,complex double * ex,complex double * ey,complex double * ez)1195 void nefld( double xob, double yob,
1196 double zob, complex double *ex,
1197 complex double *ey, complex double *ez )
1198 {
1199 int i, ix, ipr, iprx, jc, ipa;
1200 double zp, xi, ax;
1201 complex double acx, bcx, ccx;
1202
1203 *ex=CPLX_00;
1204 *ey=CPLX_00;
1205 *ez=CPLX_00;
1206 ax=0.0;
1207
1208 if( data.n != 0)
1209 {
1210 for( i = 0; i < data.n; i++ )
1211 {
1212 dataj.xj= xob- data.x[i];
1213 dataj.yj= yob- data.y[i];
1214 dataj.zj= zob- data.z[i];
1215 zp= data.cab[i]* dataj.xj+ data.sab[i] *
1216 dataj.yj+ data.salp[i]* dataj.zj;
1217
1218 if( fabs( zp) > 0.5001* data.si[i])
1219 continue;
1220
1221 zp= dataj.xj* dataj.xj+ dataj.yj* dataj.yj +
1222 dataj.zj* dataj.zj- zp* zp;
1223 dataj.xj= data.bi[i];
1224
1225 if( zp > 0.9* dataj.xj* dataj.xj)
1226 continue;
1227
1228 ax= dataj.xj;
1229 break;
1230
1231 } /* for( i = 0; i < n; i++ ) */
1232
1233 for( i = 0; i < data.n; i++ )
1234 {
1235 ix = i+1;
1236 dataj.s= data.si[i];
1237 dataj.b= data.bi[i];
1238 dataj.xj= data.x[i];
1239 dataj.yj= data.y[i];
1240 dataj.zj= data.z[i];
1241 dataj.cabj= data.cab[i];
1242 dataj.sabj= data.sab[i];
1243 dataj.salpj= data.salp[i];
1244
1245 if( dataj.iexk != 0)
1246 {
1247 ipr= data.icon1[i];
1248
1249 if (ipr > PCHCON) dataj.ind1 = 2;
1250 else if( ipr < 0 )
1251 {
1252 ipr = -ipr;
1253 iprx = ipr-1;
1254
1255 if( -data.icon1[iprx] != ix )
1256 dataj.ind1=2;
1257 else
1258 {
1259 xi= fabs( dataj.cabj* data.cab[iprx]+ dataj.sabj*
1260 data.sab[iprx]+ dataj.salpj* data.salp[iprx]);
1261 if( (xi < 0.999999) ||
1262 (fabs(data.bi[iprx]/dataj.b-1.0) > 1.0e-6) )
1263 dataj.ind1=2;
1264 else
1265 dataj.ind1=0;
1266 }
1267 } /* if( ipr < 0 ) */
1268 else if( ipr == 0 )
1269 dataj.ind1=1;
1270 else
1271 {
1272 iprx = ipr-1;
1273
1274 if( ipr != ix )
1275 {
1276 if( data.icon2[iprx] != ix )
1277 dataj.ind1=2;
1278 else
1279 {
1280 xi= fabs( dataj.cabj* data.cab[iprx]+ dataj.sabj*
1281 data.sab[iprx]+ dataj.salpj* data.salp[iprx]);
1282 if( (xi < 0.999999) ||
1283 (fabs(data.bi[iprx]/dataj.b-1.0) > 1.0e-6) )
1284 dataj.ind1=2;
1285 else
1286 dataj.ind1=0;
1287 }
1288 } /* if( ipr != ix ) */
1289 else
1290 {
1291 if( dataj.cabj* dataj.cabj +
1292 dataj.sabj* dataj.sabj > 1.0e-8)
1293 dataj.ind1=2;
1294 else
1295 dataj.ind1=0;
1296 }
1297 } /* else */
1298
1299 ipr= data.icon2[i];
1300
1301 if (ipr > PCHCON) dataj.ind2 = 2;
1302 else if( ipr < 0 )
1303 {
1304 ipr = -ipr;
1305 iprx = ipr-1;
1306
1307 if( -data.icon2[iprx] != ix )
1308 dataj.ind1=2;
1309 else
1310 {
1311 xi= fabs( dataj.cabj* data.cab[iprx]+ dataj.sabj*
1312 data.sab[iprx]+ dataj.salpj* data.salp[iprx]);
1313 if( (xi < 0.999999) ||
1314 (fabs(data.bi[iprx]/dataj.b-1.0) > 1.0e-6) )
1315 dataj.ind1=2;
1316 else
1317 dataj.ind1=0;
1318 }
1319 } /* if( ipr < 0 ) */
1320 else if( ipr == 0 )
1321 dataj.ind2=1;
1322 else
1323 {
1324 iprx = ipr-1;
1325
1326 if( ipr != ix )
1327 {
1328 if( data.icon1[iprx] != ix )
1329 dataj.ind2=2;
1330 else
1331 {
1332 xi= fabs( dataj.cabj* data.cab[iprx]+ dataj.sabj*
1333 data.sab[iprx]+ dataj.salpj* data.salp[iprx]);
1334 if( (xi < 0.999999) ||
1335 (fabs(data.bi[iprx]/dataj.b-1.0) > 1.0e-6) )
1336 dataj.ind2=2;
1337 else
1338 dataj.ind2=0;
1339 }
1340 } /* if( ipr != (i+1) ) */
1341 else
1342 {
1343 if( dataj.cabj* dataj.cabj +
1344 dataj.sabj* dataj.sabj > 1.0e-8)
1345 dataj.ind1=2;
1346 else
1347 dataj.ind1=0;
1348 }
1349
1350 } /* else */
1351
1352 } /* if( dataj.iexk != 0) */
1353
1354 efld( xob, yob, zob, ax,1);
1355 acx= cmplx( crnt.air[i], crnt.aii[i]);
1356 bcx= cmplx( crnt.bir[i], crnt.bii[i]);
1357 ccx= cmplx( crnt.cir[i], crnt.cii[i]);
1358 *ex += dataj.exk* acx+ dataj.exs* bcx+ dataj.exc* ccx;
1359 *ey += dataj.eyk* acx+ dataj.eys* bcx+ dataj.eyc* ccx;
1360 *ez += dataj.ezk* acx+ dataj.ezs* bcx+ dataj.ezc* ccx;
1361
1362 } /* for( i = 0; i < n; i++ ) */
1363
1364 if( data.m == 0)
1365 return;
1366
1367 } /* if( n != 0) */
1368
1369 jc= data.n-1;
1370 for( i = 0; i < data.m; i++ )
1371 {
1372 dataj.s= data.pbi[i];
1373 dataj.xj= data.px[i];
1374 dataj.yj= data.py[i];
1375 dataj.zj= data.pz[i];
1376 dataj.t1xj= data.t1x[i];
1377 dataj.t1yj= data.t1y[i];
1378 dataj.t1zj= data.t1z[i];
1379 dataj.t2xj= data.t2x[i];
1380 dataj.t2yj= data.t2y[i];
1381 dataj.t2zj= data.t2z[i];
1382 jc += 3;
1383 acx= dataj.t1xj* crnt.cur[jc-2]+ dataj.t1yj *
1384 crnt.cur[jc-1]+ dataj.t1zj* crnt.cur[jc];
1385 bcx= dataj.t2xj* crnt.cur[jc-2]+ dataj.t2yj *
1386 crnt.cur[jc-1]+ dataj.t2zj* crnt.cur[jc];
1387
1388 for( ipa = 0; ipa < gnd.ksymp; ipa++ )
1389 {
1390 dataj.ipgnd= ipa+1;
1391 unere( xob, yob, zob);
1392 *ex= *ex+ acx* dataj.exk+ bcx* dataj.exs;
1393 *ey= *ey+ acx* dataj.eyk+ bcx* dataj.eys;
1394 *ez= *ez+ acx* dataj.ezk+ bcx* dataj.ezs;
1395 }
1396
1397 } /* for( i = 0; i < m; i++ ) */
1398
1399 return;
1400 }
1401
1402 /*-----------------------------------------------------------------------*/
1403
1404 /* compute near e or h fields over a range of points */
nfpat(int nfeh)1405 void nfpat( int nfeh )
1406 {
1407 int i, j, kk, idx;
1408 double znrt, cth=0.0, sth=0.0, ynrt, cph=0.0, sph=0.0, yob;
1409 double xnrt, xob,zob, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
1410 complex double ex, ey, ez;
1411 double r; /* Distance of field point from xyz origin */
1412
1413 Alloc_Nearfield_Buffers( fpat.nrx, fpat.nry, fpat.nrz );
1414
1415 /* Initialize according to E/H flag */
1416 if( nfeh == 1 )
1417 near_field.max_hr = 0.0;
1418 else
1419 near_field.max_er = 0.0;
1420 near_field.r_max = 0.0;
1421
1422 idx = 0;
1423 znrt= fpat.znr- fpat.dznr;
1424 for( i = 0; i < fpat.nrz; i++ )
1425 {
1426 znrt += fpat.dznr;
1427 if( fpat.near != 0)
1428 {
1429 cth= cos( TA* znrt);
1430 sth= sin( TA* znrt);
1431 }
1432
1433 ynrt= fpat.ynr- fpat.dynr;
1434 for( j = 0; j < fpat.nry; j++ )
1435 {
1436 ynrt += fpat.dynr;
1437 if( fpat.near != 0)
1438 {
1439 cph= cos( TA* ynrt);
1440 sph= sin( TA* ynrt);
1441 }
1442
1443 xnrt= fpat.xnr- fpat.dxnr;
1444 for( kk = 0; kk < fpat.nrx; kk++ )
1445 {
1446 xnrt += fpat.dxnr;
1447 if( fpat.near != 0)
1448 {
1449 xob= xnrt* sth* cph;
1450 yob= xnrt* sth* sph;
1451 zob= xnrt* cth;
1452 }
1453 else
1454 {
1455 xob= xnrt;
1456 yob= ynrt;
1457 zob= znrt;
1458 }
1459
1460 tmp1= xob/ data.wlam;
1461 tmp2= yob/ data.wlam;
1462 tmp3= zob/ data.wlam;
1463
1464 if( nfeh == 1 ) /* Magnetic field */
1465 nhfld( tmp1, tmp2, tmp3, &ex, &ey, &ez);
1466 else /* Electric field */
1467 nefld( tmp1, tmp2, tmp3, &ex, &ey, &ez);
1468
1469 /* Calculate total field vector */
1470 Near_Field_Total( ex, ey, ez, nfeh, idx );
1471
1472 /* Save field point co-ordinates */
1473 near_field.px[idx] = (double)xob;
1474 near_field.py[idx] = (double)yob;
1475 near_field.pz[idx] = (double)zob;
1476
1477 /* Find max distance from xyz origin */
1478 r = sqrt(
1479 near_field.px[idx] * near_field.px[idx] +
1480 near_field.py[idx] * near_field.py[idx] +
1481 near_field.pz[idx] * near_field.pz[idx] );
1482 if( near_field.r_max < r )
1483 near_field.r_max = r;
1484
1485 tmp1= cabs(ex);
1486 tmp2= cang (ex);
1487 tmp3= cabs(ey);
1488 tmp4= cang (ey);
1489 tmp5= cabs(ez);
1490 tmp6= cang (ez);
1491
1492 if( nfeh == 1 ) /* Magnetic field */
1493 {
1494 near_field.hx[idx] = (double)tmp1;
1495 near_field.hy[idx] = (double)tmp3;
1496 near_field.hz[idx] = (double)tmp5;
1497 near_field.fhx[idx] = (double)(tmp2 * TA);
1498 near_field.fhy[idx] = (double)(tmp4 * TA);
1499 near_field.fhz[idx] = (double)(tmp6 * TA);
1500 }
1501 else /* Electric field */
1502 {
1503 near_field.ex[idx] = (double)tmp1;
1504 near_field.ey[idx] = (double)tmp3;
1505 near_field.ez[idx] = (double)tmp5;
1506 near_field.fex[idx] = (double)(tmp2 * TA);
1507 near_field.fey[idx] = (double)(tmp4 * TA);
1508 near_field.fez[idx] = (double)(tmp6 * TA);
1509 }
1510
1511 idx++;
1512
1513 } /* for( kk = 0; kk < fpat.nrx; kk++ ) */
1514
1515 } /* for( j = 0; j < fpat.nry; j++ ) */
1516
1517 } /* for( i = 0; i < fpat.nrz; i++ ) */
1518
1519 /* Signal new valid near field data */
1520 near_field.newer = near_field.valid = 1;
1521
1522 /* Signal new E/H pattern data */
1523 SetFlag( DRAW_NEW_EHFIELD );
1524
1525 return;
1526 }
1527
1528 /*-----------------------------------------------------------------------*/
1529
1530 /* nhfld computes the near field at specified points in space after */
1531 /* the structure currents have been computed. */
1532
nhfld(double xob,double yob,double zob,complex double * hx,complex double * hy,complex double * hz)1533 void nhfld( double xob, double yob,
1534 double zob, complex double *hx,
1535 complex double *hy, complex double *hz )
1536 {
1537 int i, jc;
1538 double ax, zp;
1539 complex double acx, bcx, ccx;
1540
1541 *hx=CPLX_00;
1542 *hy=CPLX_00;
1543 *hz=CPLX_00;
1544 ax=0.0;
1545
1546 if( data.n != 0)
1547 {
1548 for( i = 0; i < data.n; i++ )
1549 {
1550 dataj.xj= xob- data.x[i];
1551 dataj.yj= yob- data.y[i];
1552 dataj.zj= zob- data.z[i];
1553 zp= data.cab[i]* dataj.xj+ data.sab[i] *
1554 dataj.yj+ data.salp[i]* dataj.zj;
1555
1556 if( fabs( zp) > 0.5001* data.si[i])
1557 continue;
1558
1559 zp= dataj.xj* dataj.xj+ dataj.yj* dataj.yj +
1560 dataj.zj* dataj.zj- zp* zp;
1561 dataj.xj= data.bi[i];
1562
1563 if( zp > 0.9* dataj.xj* dataj.xj)
1564 continue;
1565
1566 ax= dataj.xj;
1567 break;
1568 }
1569
1570 for( i = 0; i < data.n; i++ )
1571 {
1572 dataj.s= data.si[i];
1573 dataj.b= data.bi[i];
1574 dataj.xj= data.x[i];
1575 dataj.yj= data.y[i];
1576 dataj.zj= data.z[i];
1577 dataj.cabj= data.cab[i];
1578 dataj.sabj= data.sab[i];
1579 dataj.salpj= data.salp[i];
1580 hsfld( xob, yob, zob, ax);
1581 acx= cmplx( crnt.air[i], crnt.aii[i]);
1582 bcx= cmplx( crnt.bir[i], crnt.bii[i]);
1583 ccx= cmplx( crnt.cir[i], crnt.cii[i]);
1584 *hx += dataj.exk* acx+ dataj.exs* bcx+ dataj.exc* ccx;
1585 *hy += dataj.eyk* acx+ dataj.eys* bcx+ dataj.eyc* ccx;
1586 *hz += dataj.ezk* acx+ dataj.ezs* bcx+ dataj.ezc* ccx;
1587 }
1588
1589 if( data.m == 0)
1590 return;
1591
1592 } /* if( data.n != 0) */
1593
1594 jc= data.n-1;
1595 for( i = 0; i < data.m; i++ )
1596 {
1597 dataj.s= data.pbi[i];
1598 dataj.xj= data.px[i];
1599 dataj.yj= data.py[i];
1600 dataj.zj= data.pz[i];
1601 dataj.t1xj= data.t1x[i];
1602 dataj.t1yj= data.t1y[i];
1603 dataj.t1zj= data.t1z[i];
1604 dataj.t2xj= data.t2x[i];
1605 dataj.t2yj= data.t2y[i];
1606 dataj.t2zj= data.t2z[i];
1607 hintg( xob, yob, zob);
1608 jc += 3;
1609 acx= dataj.t1xj* crnt.cur[jc-2]+ dataj.t1yj *
1610 crnt.cur[jc-1]+ dataj.t1zj* crnt.cur[jc];
1611 bcx= dataj.t2xj* crnt.cur[jc-2]+ dataj.t2yj *
1612 crnt.cur[jc-1]+ dataj.t2zj* crnt.cur[jc];
1613 *hx= *hx+ acx* dataj.exk+ bcx* dataj.exs;
1614 *hy= *hy+ acx* dataj.eyk+ bcx* dataj.eys;
1615 *hz= *hz+ acx* dataj.ezk+ bcx* dataj.ezs;
1616 }
1617
1618 return;
1619 }
1620
1621 /*-----------------------------------------------------------------------*/
1622
1623 /* integrate over patches at wire connection point */
pcint(double xi,double yi,double zi,double cabi,double sabi,double salpi,complex double * e)1624 void pcint( double xi, double yi, double zi,
1625 double cabi, double sabi, double salpi,
1626 complex double *e )
1627 {
1628 int nint, i1, i2;
1629 double d, ds, da, gcon, fcon, xxj, xyj, xzj, xs, s1;
1630 double xss, yss, zss, s2x, s2, g1, g2, g3, g4, f2, f1;
1631 complex double e1, e2, e3, e4, e5, e6, e7, e8, e9;
1632
1633 nint = 10;
1634 d= sqrt( dataj.s)*.5;
1635 ds=4.0* d/ (double) nint;
1636 da= ds* ds;
1637 gcon=1.0/ dataj.s;
1638 fcon=1.0/(2.0* TP* d);
1639 xxj= dataj.xj;
1640 xyj= dataj.yj;
1641 xzj= dataj.zj;
1642 xs= dataj.s;
1643 dataj.s= da;
1644 s1= d+ ds*.5;
1645 xss= dataj.xj+ s1*( dataj.t1xj+ dataj.t2xj);
1646 yss= dataj.yj+ s1*( dataj.t1yj+ dataj.t2yj);
1647 zss= dataj.zj+ s1*( dataj.t1zj+ dataj.t2zj);
1648 s1= s1+ d;
1649 s2x= s1;
1650 e1=CPLX_00;
1651 e2=CPLX_00;
1652 e3=CPLX_00;
1653 e4=CPLX_00;
1654 e5=CPLX_00;
1655 e6=CPLX_00;
1656 e7=CPLX_00;
1657 e8=CPLX_00;
1658 e9=CPLX_00;
1659
1660 for( i1 = 0; i1 < nint; i1++ )
1661 {
1662 s1= s1- ds;
1663 s2= s2x;
1664 xss= xss- ds* dataj.t1xj;
1665 yss= yss- ds* dataj.t1yj;
1666 zss= zss- ds* dataj.t1zj;
1667 dataj.xj= xss;
1668 dataj.yj= yss;
1669 dataj.zj= zss;
1670
1671 for( i2 = 0; i2 < nint; i2++ )
1672 {
1673 s2= s2- ds;
1674 dataj.xj= dataj.xj- ds* dataj.t2xj;
1675 dataj.yj= dataj.yj- ds* dataj.t2yj;
1676 dataj.zj= dataj.zj- ds* dataj.t2zj;
1677 unere( xi, yi, zi);
1678 dataj.exk= dataj.exk* cabi+ dataj.eyk *
1679 sabi+ dataj.ezk* salpi;
1680 dataj.exs= dataj.exs* cabi+ dataj.eys *
1681 sabi+ dataj.ezs* salpi;
1682 g1=( d+ s1)*( d+ s2)* gcon;
1683 g2=( d- s1)*( d+ s2)* gcon;
1684 g3=( d- s1)*( d- s2)* gcon;
1685 g4=( d+ s1)*( d- s2)* gcon;
1686 f2=( s1* s1+ s2* s2)* TP;
1687 f1= s1/ f2-( g1- g2- g3+ g4)* fcon;
1688 f2= s2/ f2-( g1+ g2- g3- g4)* fcon;
1689 e1= e1+ dataj.exk* g1;
1690 e2= e2+ dataj.exk* g2;
1691 e3= e3+ dataj.exk* g3;
1692 e4= e4+ dataj.exk* g4;
1693 e5= e5+ dataj.exs* g1;
1694 e6= e6+ dataj.exs* g2;
1695 e7= e7+ dataj.exs* g3;
1696 e8= e8+ dataj.exs* g4;
1697 e9= e9+ dataj.exk* f1+ dataj.exs* f2;
1698
1699 } /* for( i2 = 0; i2 < nint; i2++ ) */
1700
1701 } /* for( i1 = 0; i1 < nint; i1++ ) */
1702
1703 e[0]= e1;
1704 e[1]= e2;
1705 e[2]= e3;
1706 e[3]= e4;
1707 e[4]= e5;
1708 e[5]= e6;
1709 e[6]= e7;
1710 e[7]= e8;
1711 e[8]= e9;
1712 dataj.xj= xxj;
1713 dataj.yj= xyj;
1714 dataj.zj= xzj;
1715 dataj.s= xs;
1716
1717 return;
1718 }
1719
1720 /*-----------------------------------------------------------------------*/
1721
1722 /* calculates the electric field due to unit current */
1723 /* in the t1 and t2 directions on a patch */
unere(double xob,double yob,double zob)1724 void unere( double xob, double yob, double zob )
1725 {
1726 double zr, t1zr, t2zr, rx, ry, rz, r, tt1;
1727 double tt2, rt, xymag, px, py, cth, r2;
1728 complex double er, q1, q2, rrv, rrh, edp;
1729
1730 zr= dataj.zj;
1731 t1zr= dataj.t1zj;
1732 t2zr= dataj.t2zj;
1733
1734 if( dataj.ipgnd == 2)
1735 {
1736 zr= -zr;
1737 t1zr= -t1zr;
1738 t2zr= -t2zr;
1739 }
1740
1741 rx= xob- dataj.xj;
1742 ry= yob- dataj.yj;
1743 rz= zob- zr;
1744 r2= rx* rx+ ry* ry+ rz* rz;
1745
1746 if( r2 <= 1.0e-20)
1747 {
1748 dataj.exk=CPLX_00;
1749 dataj.eyk=CPLX_00;
1750 dataj.ezk=CPLX_00;
1751 dataj.exs=CPLX_00;
1752 dataj.eys=CPLX_00;
1753 dataj.ezs=CPLX_00;
1754 return;
1755 }
1756
1757 r= sqrt( r2);
1758 tt1= -TP* r;
1759 tt2= tt1* tt1;
1760 rt= r2* r;
1761 er= cmplx( sin( tt1),- cos( tt1))*( CONST2* dataj.s);
1762 q1= cmplx( tt2-1.0, tt1)* er/ rt;
1763 q2= cmplx(3.0- tt2,-3.0* tt1)* er/( rt* r2);
1764 er = q2*( dataj.t1xj* rx+ dataj.t1yj* ry+ t1zr* rz);
1765 dataj.exk= q1* dataj.t1xj+ er* rx;
1766 dataj.eyk= q1* dataj.t1yj+ er* ry;
1767 dataj.ezk= q1* t1zr+ er* rz;
1768 er= q2*( dataj.t2xj* rx+ dataj.t2yj* ry+ t2zr* rz);
1769 dataj.exs= q1* dataj.t2xj+ er* rx;
1770 dataj.eys= q1* dataj.t2yj+ er* ry;
1771 dataj.ezs= q1* t2zr+ er* rz;
1772
1773 if( dataj.ipgnd == 1)
1774 return;
1775
1776 if( gnd.iperf == 1)
1777 {
1778 dataj.exk= -dataj.exk;
1779 dataj.eyk= -dataj.eyk;
1780 dataj.ezk= -dataj.ezk;
1781 dataj.exs= -dataj.exs;
1782 dataj.eys= -dataj.eys;
1783 dataj.ezs= -dataj.ezs;
1784 return;
1785 }
1786
1787 xymag= sqrt( rx* rx+ ry* ry);
1788 if( xymag <= 1.0e-6)
1789 {
1790 px=0.0;
1791 py=0.0;
1792 cth=1.0;
1793 rrv=CPLX_10;
1794 }
1795 else
1796 {
1797 px= -ry/ xymag;
1798 py= rx/ xymag;
1799 cth= rz/ sqrt( xymag* xymag+ rz* rz);
1800 rrv= csqrt(1.0- gnd.zrati* gnd.zrati*(1.0- cth* cth));
1801 }
1802
1803 rrh= gnd.zrati* cth;
1804 rrh=( rrh- rrv)/( rrh+ rrv);
1805 rrv= gnd.zrati* rrv;
1806 rrv=-( cth- rrv)/( cth+ rrv);
1807 edp=( dataj.exk* px+ dataj.eyk* py)*( rrh- rrv);
1808 dataj.exk= dataj.exk* rrv+ edp* px;
1809 dataj.eyk= dataj.eyk* rrv+ edp* py;
1810 dataj.ezk= dataj.ezk* rrv;
1811 edp=( dataj.exs* px+ dataj.eys* py)*( rrh- rrv);
1812 dataj.exs= dataj.exs* rrv+ edp* px;
1813 dataj.eys= dataj.eys* rrv+ edp* py;
1814 dataj.ezs= dataj.ezs* rrv;
1815
1816 return;
1817 }
1818
1819 /*-----------------------------------------------------------------------*/
1820
1821 /* Near_Field_Total()
1822 *
1823 * Calculates the value of Total Near Field vector
1824 */
1825 void
Near_Field_Total(complex double ex,complex double ey,complex double ez,int nfeh,int idx)1826 Near_Field_Total(
1827 complex double ex,
1828 complex double ey,
1829 complex double ez,
1830 int nfeh, int idx )
1831 {
1832 /* Display a time-frozen "snapshot" of near field */
1833 if( isFlagSet(NEAREH_SNAPSHOT) )
1834 {
1835 if( nfeh == 1 ) /* Magnetic field */
1836 {
1837 /* Near magnetic field components */
1838 near_field.hrx[idx] = (double)creal(ex);
1839 near_field.hry[idx] = (double)creal(ey);
1840 near_field.hrz[idx] = (double)creal(ez);
1841
1842 /* Near total magnetic field vector*/
1843 near_field.hr[idx] = sqrt(
1844 near_field.hrx[idx] * near_field.hrx[idx] +
1845 near_field.hry[idx] * near_field.hry[idx] +
1846 near_field.hrz[idx] * near_field.hrz[idx] );
1847 if( near_field.max_hr < near_field.hr[idx] )
1848 near_field.max_hr = near_field.hr[idx];
1849 }
1850 else /* Electric field */
1851 {
1852 /* Near electric field components */
1853 /* Near electric field components */
1854 near_field.erx[idx] = (double)creal(ex);
1855 near_field.ery[idx] = (double)creal(ey);
1856 near_field.erz[idx] = (double)creal(ez);
1857
1858 /* Near total electric field vector */
1859 near_field.er[idx] = sqrt(
1860 near_field.erx[idx] * near_field.erx[idx] +
1861 near_field.ery[idx] * near_field.ery[idx] +
1862 near_field.erz[idx] * near_field.erz[idx] );
1863 if( near_field.max_er < near_field.er[idx] )
1864 near_field.max_er = near_field.er[idx];
1865 } /* if( nfeh == 1 ) */
1866
1867 } /* if( isFlagSet(NEAREH_SNAPSHOT) ) */
1868 else /* Display Total near field vector peak */
1869 {
1870 double
1871 exm, eym, ezm, /* Near field magnitude in x, y, z */
1872 exm2, eym2, ezm2, /* Near field magnitude^2 in x, y, z */
1873 fx, fy, fz, /* Time phase of near field vectors */
1874 fx2, fy2, fz2, /* Time phase of near field vectors*2 */
1875 cp, sp, tp, wt; /* Some values needed in calculations */
1876
1877 exm = (double)cabs(ex);
1878 eym = (double)cabs(ey);
1879 ezm = (double)cabs(ez);
1880 /* Near total electric field vector */
1881 fx = (double)cang(ex)/(double)TD;
1882 fy = (double)cang(ey)/(double)TD;
1883 fz = (double)cang(ez)/(double)TD;
1884
1885 fx2 = fx * 2.0;
1886 fy2 = fy * 2.0;
1887 fz2 = fz * 2.0;
1888
1889 exm2 = exm*exm;
1890 eym2 = eym*eym;
1891 ezm2 = ezm*ezm;
1892
1893 cp = exm2*cos(fx2) + eym2*cos(fy2) + ezm2*cos(fz2);
1894 sp = exm2*sin(fx2) + eym2*sin(fy2) + ezm2*sin(fz2);
1895 tp = sqrt(cp*cp + sp*sp);
1896 wt = atan2(-sp, cp)/2.0;
1897
1898 if( nfeh == 1 ) /* Magnetic field */
1899 {
1900 /* Near magnetic field components */
1901 near_field.hrx[idx] = exm * cos(wt + fx);
1902 near_field.hry[idx] = eym * cos(wt + fy);
1903 near_field.hrz[idx] = ezm * cos(wt + fz);
1904
1905 /* Near total magnetic field vector, peak value */
1906 near_field.hr[idx] = sqrt( (exm2 + eym2 + ezm2 + tp)/2.0 );
1907 if( near_field.max_hr < near_field.hr[idx] )
1908 near_field.max_hr = near_field.hr[idx];
1909 }
1910 else /* Electric field */
1911 {
1912 /* Near electric field components */
1913 near_field.erx[idx] = exm * cos(wt + fx);
1914 near_field.ery[idx] = eym * cos(wt + fy);
1915 near_field.erz[idx] = ezm * cos(wt + fz);
1916
1917 /* Near total electric field vector, peak value */
1918 near_field.er[idx] = sqrt( (exm2 + eym2 + ezm2 + tp)/2.0 );
1919 if( near_field.max_er < near_field.er[idx] )
1920 near_field.max_er = near_field.er[idx];
1921 }
1922 }
1923
1924 } /* Near_Field_Total() */
1925
1926 /*-----------------------------------------------------------------------*/
1927
1928