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 /* last change: pgm 8 nov 2000 1:04 pm
20 program somnec(input,output,tape21)
21
22 program to generate nec interpolation grids for fields due to
23 ground. field components are computed by numerical evaluation
24 of modified sommerfeld integrals.
25
26 somnec2d is a double precision version of somnec for use with
27 nec2d. an alternate version (somnec2sd) is also provided in which
28 computation is in single precision but the output file is written
29 in double precision for use with nec2d. somnec2sd runs about twic
30 as fast as the full double precision somnec2d. the difference
31 between nec2d results using a for021 file from this code rather
32 than from somnec2sd was insignficant in the cases tested.
33
34 changes made by j bergervoet, 31-5-95:
35 parameter 0.0 --> 0.0d0 in calling of routine test
36 status of output files set to 'unknown' */
37
38 #include "somnec.h"
39 #include "shared.h"
40
41 #include "clog.c"
42 #include "cexp.c"
43
44 /* common /evlcom/ */
45 static int jh;
46 static double ck2, ck2sq, tkmag, tsmag, ck1r, zph, rho;
47 static complex double ct1, ct2, ct3, ck1, ck1sq, cksm;
48
49 /* common /cntour/ */
50 static complex double a, b;
51
52 /*-----------------------------------------------------------------------*/
53
54 /* This is the "main" of somnec */
somnec(double epr,double sig,double fmhz)55 void somnec( double epr, double sig, double fmhz )
56 {
57 int k, nth, ith, irs, ir, nr;
58 double wlam, dr, dth=0.0, r, rk, thet, tfac1, tfac2;
59 complex double erv, ezv, erh, eph, cl1, cl2, con;
60
61 static gboolean first_call = TRUE;
62
63 if( first_call )
64 {
65 first_call = FALSE;
66
67 /* Allocate some buffers */
68 size_t mreq = sizeof(complex double) * 11 * 10 * 4;
69 mem_alloc( (void **)&ggrid.ar1, mreq, "in somnec.c");
70 mreq = sizeof(complex double) * 17 * 5 * 4;
71 mem_alloc( (void **)&ggrid.ar2, mreq, "in somnec.c");
72 mreq = sizeof(complex double) * 9 * 8 * 4;
73 mem_alloc( (void **)&ggrid.ar3, mreq, "in somnec.c");
74
75 mreq = sizeof(int) * 3;
76 mem_alloc( (void **)&ggrid.nxa, mreq, "in somnec.c");
77 mem_alloc( (void **)&ggrid.nya, mreq, "in somnec.c");
78
79 mreq = sizeof(double) * 3;
80 mem_alloc( (void **)&ggrid.dxa, mreq, "in somnec.c");
81 mem_alloc( (void **)&ggrid.dya, mreq, "in somnec.c");
82 mem_alloc( (void **)&ggrid.xsa, mreq, "in somnec.c");
83 mem_alloc( (void **)&ggrid.ysa, mreq, "in somnec.c");
84
85 /* Initialize ground grid parameters for somnec */
86 ggrid.nxa[0] = 11;
87 ggrid.nxa[1] = 17;
88 ggrid.nxa[2] = 9;
89
90 ggrid.nya[0] = 10;
91 ggrid.nya[1] = 5;
92 ggrid.nya[2] = 8;
93
94 ggrid.dxa[0] = .02;
95 ggrid.dxa[1] = .05;
96 ggrid.dxa[2] = .1;
97
98 ggrid.dya[0] = .1745329252;
99 ggrid.dya[1] = .0872664626;
100 ggrid.dya[2] = .1745329252;
101
102 ggrid.xsa[0] = 0.0;
103 ggrid.xsa[1] = .2;
104 ggrid.xsa[2] = .2;
105
106 ggrid.ysa[0] = 0.0;
107 ggrid.ysa[1] = 0.0;
108 ggrid.ysa[2] = .3490658504;
109
110 } /* if( first_call ) */
111
112 if(sig >= 0.0)
113 {
114 wlam=CVEL/fmhz;
115 ggrid.epscf=cmplx(epr,-sig*wlam*59.96);
116 }
117 else ggrid.epscf=cmplx(epr,sig);
118
119 ck2=TP;
120 ck2sq=ck2*ck2;
121
122 /* sommerfeld integral evaluation uses exp(-jwt),
123 * nec uses exp(+jwt), hence need conjg(ggrid.epscf).
124 * conjugate of fields occurs in subroutine evlua. */
125
126 ck1sq=ck2sq*conj(ggrid.epscf);
127 ck1=csqrt(ck1sq);
128 ck1r=creal(ck1);
129 tkmag=100.0*cabs(ck1);
130 tsmag=100.0*creal( ck1*conj(ck1) );
131 cksm=ck2sq/(ck1sq+ck2sq);
132 ct1=.5*(ck1sq-ck2sq);
133 erv=ck1sq*ck1sq;
134 ezv=ck2sq*ck2sq;
135 ct2=.125*(erv-ezv);
136 erv *= ck1sq;
137 ezv *= ck2sq;
138 ct3=.0625*(erv-ezv);
139
140 /* loop over 3 grid regions */
141 for( k = 0; k < 3; k++ )
142 {
143 nr=ggrid.nxa[k];
144 nth=ggrid.nya[k];
145 dr=ggrid.dxa[k];
146 dth=ggrid.dya[k];
147 r=ggrid.xsa[k]-dr;
148 irs=1;
149 if(k == 0)
150 {
151 r=ggrid.xsa[k];
152 irs=2;
153 }
154
155 /* loop over r. (r=sqrt(rho**2 + (z+h)**2)) */
156 for( ir = irs-1; ir < nr; ir++ )
157 {
158 r += dr;
159 thet = ggrid.ysa[k]-dth;
160
161 /* loop over theta. (theta=atan((z+h)/rho)) */
162 for( ith = 0; ith < nth; ith++ )
163 {
164 thet += dth;
165 rho=r*cos(thet);
166 zph=r*sin(thet);
167 if(rho < 1.0e-7)
168 rho=1.0e-8;
169 if(zph < 1.0e-7)
170 zph=0.0;
171
172 evlua( &erv, &ezv, &erh, &eph );
173
174 rk=ck2*r;
175 con=-CONST1*r/cmplx(cos(rk),-sin(rk));
176
177 switch( k )
178 {
179 case 0:
180 ggrid.ar1[ir+ith*11+ 0]=erv*con;
181 ggrid.ar1[ir+ith*11+110]=ezv*con;
182 ggrid.ar1[ir+ith*11+220]=erh*con;
183 ggrid.ar1[ir+ith*11+330]=eph*con;
184 break;
185
186 case 1:
187 ggrid.ar2[ir+ith*17+ 0]=erv*con;
188 ggrid.ar2[ir+ith*17+ 85]=ezv*con;
189 ggrid.ar2[ir+ith*17+170]=erh*con;
190 ggrid.ar2[ir+ith*17+255]=eph*con;
191 break;
192
193 case 2:
194 ggrid.ar3[ir+ith*9+ 0]=erv*con;
195 ggrid.ar3[ir+ith*9+ 72]=ezv*con;
196 ggrid.ar3[ir+ith*9+144]=erh*con;
197 ggrid.ar3[ir+ith*9+216]=eph*con;
198
199 } /* switch( k ) */
200
201 } /* for( ith = 0; ith < nth; ith++ ) */
202
203 } /* for( ir = irs-1; ir < nr; ir++; ) */
204
205 } /* for( k = 0; k < 3; k++; ) */
206
207 /* fill grid 1 for r equal to zero. */
208 cl2=-CONST4*(ggrid.epscf-1.0)/(ggrid.epscf+1.0);
209 cl1=cl2/(ggrid.epscf+1.0);
210 ezv=ggrid.epscf*cl1;
211 thet=-dth;
212 nth=ggrid.nya[0];
213
214 for( ith = 0; ith < nth; ith++ )
215 {
216 thet += dth;
217 if( (ith+1) != nth )
218 {
219 tfac2=cos(thet);
220 tfac1=(1.0-sin(thet))/tfac2;
221 tfac2=tfac1/tfac2;
222 erv=ggrid.epscf*cl1*tfac1;
223 erh=cl1*(tfac2-1.0)+cl2;
224 eph=cl1*tfac2-cl2;
225 }
226 else
227 {
228 erv=0.0;
229 erh=cl2-.5*cl1;
230 eph=-erh;
231 }
232
233 ggrid.ar1[0+ith*11+ 0]=erv;
234 ggrid.ar1[0+ith*11+110]=ezv;
235 ggrid.ar1[0+ith*11+220]=erh;
236 ggrid.ar1[0+ith*11+330]=eph;
237 }
238
239 return;
240 }
241
242 /*-----------------------------------------------------------------------*/
243
244 /* bessel evaluates the zero-order bessel function */
245 /* and its derivative for complex argument z. */
bessel(complex double z,complex double * j0,complex double * j0p)246 void bessel( complex double z,
247 complex double *j0, complex double *j0p )
248 {
249 int k, ib;
250 static int *m = NULL, init = FALSE;
251 static double *a1 = NULL, *a2 = NULL;
252 double zms;
253 complex double p0z, p1z, q0z, q1z, zi, zi2, zk, cz, sz;
254 complex double j0x=CPLX_00, j0px=CPLX_00;
255
256 /* initialization of constants */
257 if( !init )
258 {
259 int i;
260
261 size_t mreq = 101 * sizeof(double);
262 mem_alloc( (void **)&m, mreq, "in somnec.c");
263 mem_alloc( (void **)&a1, mreq, "in somnec.c");
264 mem_alloc( (void **)&a2, mreq, "in somnec.c");
265
266 for( k = 1; k <= 25; k++ )
267 {
268 i = k-1;
269 a1[i]=-.25/(k*k);
270 a2[i]=1.0/(k+1.0);
271 }
272
273 for( i = 1; i <= 101; i++ )
274 {
275 double tst=1.0;
276 for( k = 0; k < 24; k++ )
277 {
278 init = k;
279 tst *= -i*a1[k];
280 if( tst < 1.0e-6 )
281 break;
282 }
283
284 m[i-1] = init+1;
285 } /* for( i = 1; i<= 101; i++ ) */
286
287 init = TRUE;
288 } /* if(init == 0) */
289
290 zms=creal( z*conj(z) );
291 if(zms <= 1.0e-12)
292 {
293 *j0=CPLX_10;
294 *j0p=-.5*z;
295 return;
296 }
297
298 ib=0;
299 if(zms <= 37.21)
300 {
301 int iz, miz;
302
303 if(zms > 36.0)
304 ib=1;
305
306 /* series expansion */
307 iz=(int)zms;
308 miz=m[iz];
309 *j0=CPLX_10;
310 *j0p=*j0;
311 zk=*j0;
312 zi=z*z;
313
314 for( k = 0; k < miz; k++ )
315 {
316 zk *= a1[k]*zi;
317 *j0 += zk;
318 *j0p += a2[k]*zk;
319 }
320 *j0p *= -.5*z;
321
322 if(ib == 0)
323 return;
324
325 j0x=*j0;
326 j0px=*j0p;
327 }
328
329 /* asymptotic expansion */
330 zi=1.0/z;
331 zi2=zi*zi;
332 p0z=1.0+(P20*zi2-P10)*zi2;
333 p1z=1.0+(P11-P21*zi2)*zi2;
334 q0z=(Q20*zi2-Q10)*zi;
335 q1z=(Q11-Q21*zi2)*zi;
336 zk=cexp(CPLX_01*(z-POF));
337 zi2=1.9/zk;
338 cz=.5*(zk+zi2);
339 sz=CPLX_01*.5*(zi2-zk);
340 zk=C3*csqrt(zi);
341 *j0=zk*(p0z*cz-q0z*sz);
342 *j0p=-zk*(p1z*sz+q1z*cz);
343
344 if(ib == 0)
345 return;
346
347 zms=cos((sqrt(zms)-6.0)*PI10);
348 *j0=.5*(j0x*(1.0+zms)+ *j0*(1.0-zms));
349 *j0p=.5*(j0px*(1.0+zms)+ *j0p*(1.0-zms));
350
351 return;
352 }
353
354 /*-----------------------------------------------------------------------*/
355
356 /* evlua controls the integration contour in the complex */
357 /* lambda plane for evaluation of the sommerfeld integrals */
358 void
evlua(complex double * erv,complex double * ezv,complex double * erh,complex double * eph)359 evlua( complex double *erv, complex double *ezv,
360 complex double *erh, complex double *eph )
361 {
362 int i, jump;
363 static double del, slope, rmis;
364 static complex double cp1, cp2, cp3, bk, delta, delta2;
365 static complex double *sum = NULL, *ans = NULL;
366 static gboolean first_call = TRUE;
367
368 if( first_call )
369 {
370 first_call = FALSE;
371 size_t mreq = 6 * sizeof(complex double);
372 mem_alloc( (void **)&sum, mreq, "in somnec.c");
373 mem_alloc( (void **)&ans, mreq, "in somnec.c");
374 }
375
376 del=zph;
377 if( rho > del )
378 del=rho;
379
380 if(zph >= 2.0*rho)
381 {
382 /* bessel function form of sommerfeld integrals */
383 jh=0;
384 a=CPLX_00;
385 del=1.0/del;
386
387 if( del > tkmag)
388 {
389 b=cmplx(0.1*tkmag,-0.1*tkmag);
390 rom1(6,sum,2);
391 a=b;
392 b=cmplx(del,-del);
393 rom1 (6,ans,2);
394 for( i = 0; i < 6; i++ )
395 sum[i] += ans[i];
396 }
397 else
398 {
399 b=cmplx(del,-del);
400 rom1(6,sum,2);
401 }
402
403 delta=PTP*del;
404 gshank(b,delta,ans,6,sum,0,b,b);
405 ans[5] *= ck1;
406
407 /* conjugate since nec uses exp(+jwt) */
408 *erv=conj(ck1sq*ans[2]);
409 *ezv=conj(ck1sq*(ans[1]+ck2sq*ans[4]));
410 *erh=conj(ck2sq*(ans[0]+ans[5]));
411 *eph=-conj(ck2sq*(ans[3]+ans[5]));
412
413 return;
414 } /* if(zph >= 2.0*rho) */
415
416 /* hankel function form of sommerfeld integrals */
417 jh=1;
418 cp1=cmplx(0.0, 0.4*ck2);
419 cp2=cmplx(0.6*ck2, -0.2*ck2);
420 cp3=cmplx(1.02*ck2,-0.2*ck2);
421 a=cp1;
422 b=cp2;
423 rom1(6,sum,2);
424 a=cp2;
425 b=cp3;
426 rom1(6,ans,2);
427
428 for( i = 0; i < 6; i++ )
429 sum[i]=-(sum[i]+ans[i]);
430
431 /* path from imaginary axis to -infinity */
432 if(zph > .001*rho)
433 slope=rho/zph;
434 else
435 slope=1000.0;
436
437 del=PTP/del;
438 delta=cmplx(-1.0,slope)*del/sqrt(1.0+slope*slope);
439 delta2=-conj(delta);
440 gshank(cp1,delta,ans,6,sum,0,bk,bk);
441 rmis=rho*(creal(ck1)-ck2);
442
443 jump = FALSE;
444 if( (rmis >= 2.0*ck2) && (rho >= 1.0e-10) )
445 {
446 if(zph >= 1.0e-10)
447 {
448 bk=cmplx(-zph,rho)*(ck1-cp3);
449 rmis=-creal(bk)/fabs(cimag(bk));
450 if(rmis > 4.0*rho/zph)
451 jump = TRUE;
452 }
453
454 if( ! jump )
455 {
456 /* integrate up between branch cuts, then to + infinity */
457 cp1=ck1-(0.1+I*0.2);
458 cp2=cp1+0.2;
459 bk=cmplx(0.0,del);
460 gshank(cp1,bk,sum,6,ans,0,bk,bk);
461 a=cp1;
462 b=cp2;
463 rom1(6,ans,1);
464 for( i = 0; i < 6; i++ )
465 ans[i] -= sum[i];
466
467 gshank(cp3,bk,sum,6,ans,0,bk,bk);
468 gshank(cp2,delta2,ans,6,sum,0,bk,bk);
469 }
470
471 jump = TRUE;
472
473 } /* if( (rmis >= 2.0*ck2) || (rho >= 1.0e-10) ) */
474 else jump = FALSE;
475
476 if( ! jump )
477 {
478 /* integrate below branch points, then to + infinity */
479 for( i = 0; i < 6; i++ )
480 sum[i]=-ans[i];
481
482 rmis=creal(ck1)*1.01;
483 if( (ck2+1.0) > rmis )
484 rmis=ck2+1.0;
485
486 bk=cmplx(rmis,0.99*cimag(ck1));
487 delta=bk-cp3;
488 delta *= del/cabs(delta);
489 gshank(cp3,delta,ans,6,sum,1,bk,delta2);
490
491 } /* if( ! jump ) */
492
493 ans[5] *= ck1;
494
495 /* conjugate since nec uses exp(+jwt) */
496 *erv=conj(ck1sq*ans[2]);
497 *ezv=conj(ck1sq*(ans[1]+ck2sq*ans[4]));
498 *erh=conj(ck2sq*(ans[0]+ans[5]));
499 *eph=-conj(ck2sq*(ans[3]+ans[5]));
500
501 return;
502 }
503
504 /*-----------------------------------------------------------------------*/
505
506 /* fbar is sommerfeld attenuation function for numerical distance p */
507 void
fbar(complex double p,complex double * fbar)508 fbar( complex double p, complex double *fbar )
509 {
510 int i, minus;
511 double tms, sms;
512 complex double z, zs, sum, pow, term;
513
514 z= CPLX_01* csqrt( p);
515 if( cabs( z) <= 3.0)
516 {
517 /* series expansion */
518 zs= z* z;
519 sum= z;
520 pow= z;
521
522 for( i = 1; i <= 100; i++ )
523 {
524 pow= -pow* zs/ (double)i;
525 term= pow/(2.0* i+1.0);
526 sum= sum+ term;
527 tms= creal( term* conj( term));
528 sms= creal( sum* conj( sum));
529 if( tms/sms < ACCS)
530 break;
531 }
532
533 *fbar=1.0-(1.0- sum* TOSP)* z* cexp( zs)* SP;
534
535 } /* if( cabs( z) <= 3.0) */
536
537 /* asymptotic expansion */
538 if( creal( z) < 0.0)
539 {
540 minus=1;
541 z= -z;
542 }
543 else minus=0;
544
545 zs=0.5/( z* z);
546 sum=CPLX_00;
547 term=CPLX_10;
548
549 for( i = 1; i <= 6; i++ )
550 {
551 term = -term*(2.0*i -1.0)* zs;
552 sum += term;
553 }
554
555 if( minus == 1)
556 sum -= 2.0* SP* z* cexp( z* z);
557 *fbar= -sum;
558
559 }
560
561 /*-----------------------------------------------------------------------*/
562
563 /* gshank integrates the 6 sommerfeld integrals from start to */
564 /* infinity (until convergence) in lambda. at the break point, bk, */
565 /* the step increment may be changed from dela to delb. shank's */
566 /* algorithm to accelerate convergence of a slowly converging series */
567 /* is used */
568 void
gshank(complex double start,complex double dela,complex double * sum,int nans,complex double * seed,int ibk,complex double bk,complex double delb)569 gshank( complex double start, complex double dela,
570 complex double *sum, int nans, complex double *seed,
571 int ibk, complex double bk, complex double delb )
572 {
573 int ibx, j, i, jm, intx, inx, brk=0, idx;
574 static double rbk, amg, den, denm;
575 complex double a1, a2, as1, as2, del, aa;
576 static complex double *q1 = NULL, *q2 = NULL;
577 static complex double *ans1 = NULL, *ans2 = NULL;
578 static gboolean first_call = TRUE;
579
580 if( first_call )
581 {
582 first_call = FALSE;
583 size_t mreq = 6 * 20 * sizeof(complex double);
584 mem_alloc( (void **)&q1, mreq, "in somnec.c");
585 mem_alloc( (void **)&q2, mreq, "in somnec.c");
586 mreq = 6 * sizeof(complex double);
587 mem_alloc( (void **)&ans1, mreq, "in somnec.c");
588 mem_alloc( (void **)&ans2, mreq, "in somnec.c");
589 }
590
591 rbk=creal(bk);
592 del=dela;
593 if(ibk == 0)
594 ibx=1;
595 else
596 ibx=0;
597
598 for( i = 0; i < nans; i++ )
599 ans2[i]=seed[i];
600
601 b=start;
602 for( intx = 1; intx <= MAXH; intx++ )
603 {
604 inx=intx-1;
605 a=b;
606 b += del;
607
608 if( (ibx == 0) && (creal(b) >= rbk) )
609 {
610 /* hit break point. reset seed and start over. */
611 ibx=1;
612 b=bk;
613 del=delb;
614 rom1(nans,sum,2);
615 for( i = 0; i < nans; i++ )
616 ans2[i] += sum[i];
617 intx = 0;
618 continue;
619 } /* if( (ibx == 0) && (creal(b) >= rbk) ) */
620
621 rom1(nans,sum,2);
622 for( i = 0; i < nans; i++ )
623 ans1[i] = ans2[i]+sum[i];
624 a=b;
625 b += del;
626
627 if( (ibx == 0) && (creal(b) >= rbk) )
628 {
629 /* hit break point. reset seed and start over. */
630 ibx=2;
631 b=bk;
632 del=delb;
633 rom1(nans,sum,2);
634 for( i = 0; i < nans; i++ )
635 ans2[i] = ans1[i]+sum[i];
636 intx = 0;
637 continue;
638
639 } /* if( (ibx == 0) && (creal(b) >= rbk) ) */
640
641 rom1(nans,sum,2);
642 for( i = 0; i < nans; i++ )
643 ans2[i]=ans1[i]+sum[i];
644
645 den=0.0;
646 for( i = 0; i < nans; i++ )
647 {
648 as1=ans1[i];
649 as2=ans2[i];
650
651 if(intx >= 2)
652 {
653 for( j = 1; j < intx; j++ )
654 {
655 jm=j-1;
656 idx = i + 6*jm;
657 aa=q2[idx];
658 a1=q1[idx]+as1-2.0*aa;
659
660 if( (creal(a1) != 0.0) || (cimag(a1) != 0.0) )
661 {
662 a2=aa-q1[idx];
663 a1=q1[idx]-a2*a2/a1;
664 }
665 else a1=q1[idx];
666
667 a2=aa+as2-2.0*as1;
668 if( (creal(a2) != 0.0) || (cimag(a2) != 0.0) )
669 a2=aa-(as1-aa)*(as1-aa)/a2;
670 else
671 a2=aa;
672
673 q1[idx]=as1;
674 q2[idx]=as2;
675 as1=a1;
676 as2=a2;
677
678 } /* for( j = 1; i < intx; i++ ) */
679
680 } /* if(intx >= 2) */
681
682 idx = i + 6*(intx-1);
683 q1[idx]=as1;
684 q2[idx]=as2;
685 amg=fabs(creal(as2))+fabs(cimag(as2));
686 if(amg > den)
687 den=amg;
688
689 } /* for( i = 0; i < nans; i++ ) */
690
691 denm=1.0e-3*den*CRIT;
692 jm=intx-3;
693 if(jm < 1)
694 jm=1;
695
696 for( j = jm-1; j < intx; j++ )
697 {
698 brk = FALSE;
699 for( i = 0; i < nans; i++ )
700 {
701 idx = i +6*j;
702 a1=q2[idx];
703 den=(fabs(creal(a1))+fabs(cimag(a1)))*CRIT;
704 if(den < denm)
705 den=denm;
706 a1=q1[idx]-a1;
707 amg=fabs(creal(a1)+fabs(cimag(a1)));
708 if(amg > den)
709 {
710 brk = TRUE;
711 break;
712 }
713
714 } /* for( i = 0; i < nans; i++ ) */
715
716 if( brk ) break;
717
718 } /* for( j = jm-1; j < intx; j++ ) */
719
720 if( ! brk )
721 {
722 for( i = 0; i < nans; i++ )
723 {
724 idx = i + 6*inx;
725 sum[i]=0.5*(q1[idx]+q2[idx]);
726 }
727 return;
728 }
729
730 } /* for( intx = 1; intx <= maxh; intx++ ) */
731
732 /* No convergence */
733 fprintf( stderr, "xnec2c: gshank(): No convergence\n" );
734 stop( _("gshank(): No convergencn"), ERR_STOP );
735 return;
736 }
737
738 /*-----------------------------------------------------------------------*/
739
740 /* hankel evaluates hankel function of the first kind, */
741 /* order zero, and its derivative for complex argument z */
742 void
hankel(complex double z,complex double * h0,complex double * h0p)743 hankel( complex double z,
744 complex double *h0, complex double *h0p )
745 {
746 int k, ib;
747 static int *m = NULL, init = FALSE;
748 static double *a1 = NULL, *a2 = NULL, *a3 = NULL;
749 static double *a4, psi, tst, zms;
750 complex double clogz, j0, j0p, p0z, p1z, q0z, q1z;
751 complex double y0 = CPLX_00, y0p = CPLX_00, zi, zi2, zk;
752 static gboolean first_call = TRUE;
753
754 if( first_call )
755 {
756 first_call = FALSE;
757 size_t mreq = 101 * sizeof(int);
758 mem_alloc( (void **)&m, mreq, "in somnec.c");
759 mreq = 25 * sizeof(double);
760 mem_alloc( (void **)&a1, mreq, "in somnec.c");
761 mem_alloc( (void **)&a2, mreq, "in somnec.c");
762 mem_alloc( (void **)&a3, mreq, "in somnec.c");
763 mem_alloc( (void **)&a4, mreq, "in somnec.c");
764 }
765
766 /* initialization of constants */
767 if( ! init )
768 {
769 int i;
770
771 psi=-GAMMA;
772 for( k = 1; k <= 25; k++ )
773 {
774 i = k-1;
775 a1[i]=-0.25/(k*k);
776 a2[i]=1.0/(k+1.0);
777 psi += 1.0/k;
778 a3[i]=psi+psi;
779 a4[i]=(psi+psi+1.0/(k+1.0))/(k+1.0);
780 }
781
782 for( i = 1; i <= 101; i++ )
783 {
784 tst=1.0;
785 for( k = 0; k < 24; k++ )
786 {
787 init = k;
788 tst *= -i*a1[k];
789 if(tst*a3[k] < 1.0e-6)
790 break;
791 }
792 m[i-1]=init+1;
793 }
794
795 init = TRUE;
796
797 } /* if( ! init ) */
798
799 zms=creal( z*conj(z) );
800 if(zms == 0.0)
801 {
802 fprintf( stderr, "xnec2c: hankel(): Hankel not valid for z = 0\n" );
803 stop( _("hankel(): Hankel not valid for z = 0"), ERR_STOP );
804 }
805
806 ib=0;
807 if(zms <= 16.81)
808 {
809 int iz, miz;
810
811 if(zms > 16.0)
812 ib=1;
813
814 /* series expansion */
815 iz=(int)zms;
816 miz=m[iz];
817 j0=CPLX_10;
818 j0p=j0;
819 y0=CPLX_00;
820 y0p=y0;
821 zk=j0;
822 zi=z*z;
823
824 for( k = 0; k < miz; k++ )
825 {
826 zk *= a1[k]*zi;
827 j0 += zk;
828 j0p += a2[k]*zk;
829 y0 += a3[k]*zk;
830 y0p += a4[k]*zk;
831 }
832
833 j0p *= -0.5*z;
834 clogz=clog(0.5*z);
835 y0=(2.0*j0*clogz-y0)/PI+C2;
836 y0p=(2.0/z+2.0*j0p*clogz+0.5*y0p*z)/PI+C1*z;
837 *h0=j0+CPLX_01*y0;
838 *h0p=j0p+CPLX_01*y0p;
839
840 if(ib == 0) return;
841
842 y0=*h0;
843 y0p=*h0p;
844 } /* if(zms <= 16.81) */
845
846 /* asymptotic expansion */
847 zi=1.0/z;
848 zi2=zi*zi;
849 p0z=1.0+(P20*zi2-P10)*zi2;
850 p1z=1.0+(P11-P21*zi2)*zi2;
851 q0z=(Q20*zi2-Q10)*zi;
852 q1z=(Q11-Q21*zi2)*zi;
853 zk=cexp(CPLX_01*(z-POF))*csqrt(zi)*C3;
854 *h0=zk*(p0z+CPLX_01*q0z);
855 *h0p=CPLX_01*zk*(p1z+CPLX_01*q1z);
856
857 if(ib == 0) return;
858
859 zms=cos((sqrt(zms)-4.0)*31.41592654);
860 *h0=0.5*(y0*(1.0+zms)+ *h0*(1.0-zms));
861 *h0p=0.5*(y0p*(1.0+zms)+ *h0p*(1.0-zms));
862
863 return;
864 }
865
866 /*-----------------------------------------------------------------------*/
867
868 /* compute integration parameter xlam=lambda from parameter t. */
lambda(double t,complex double * xlam,complex double * dxlam)869 void lambda( double t,
870 complex double *xlam, complex double *dxlam )
871 {
872 *dxlam=b-a;
873 *xlam=a+*dxlam*t;
874 return;
875 }
876
877 /*-----------------------------------------------------------------------*/
878
879 /* rom1 integrates the 6 sommerfeld integrals from a to b in lambda. */
880 /* the method of variable interval width romberg integration is used. */
rom1(int n,complex double * sum,int nx)881 void rom1( int n, complex double *sum, int nx )
882 {
883 int jump, lstep, nogo, i, ns, nt;
884 static double z, ze, s, ep, zend, dz=0.0, dzot=0.0, tr, ti;
885 static complex double t00, t11, t02;
886 static complex double *g1 = NULL, *g2 = NULL;
887 static complex double *g3 = NULL, *g4 = NULL;
888 static complex double *g5 = NULL, *t01 = NULL;
889 static complex double *t10 = NULL, *t20 = NULL;
890 static gboolean first_call = TRUE;
891
892 if( first_call )
893 {
894 first_call = FALSE;
895 size_t mreq = 6 * sizeof(complex double);
896 mem_alloc( (void **)&g1, mreq, "in somnec.c");
897 mem_alloc( (void **)&g2, mreq, "in somnec.c");
898 mem_alloc( (void **)&g3, mreq, "in somnec.c");
899 mem_alloc( (void **)&g4, mreq, "in somnec.c");
900 mem_alloc( (void **)&g5, mreq, "in somnec.c");
901 mem_alloc( (void **)&t01, mreq, "in somnec.c");
902 mem_alloc( (void **)&t10, mreq, "in somnec.c");
903 mem_alloc( (void **)&t20, mreq, "in somnec.c");
904 }
905
906 lstep=0;
907 z=0.0;
908 ze=1.0;
909 s=1.0;
910 ep=s/(1.0e4*NM);
911 zend=ze-ep;
912 for( i = 0; i < n; i++ )
913 sum[i]=CPLX_00;
914 ns=nx;
915 nt=0;
916 saoa(z,g1);
917
918 jump = FALSE;
919 while( TRUE )
920 {
921 if( ! jump )
922 {
923 dz=s/ns;
924 if( (z+dz) > ze )
925 {
926 dz=ze-z;
927 if( dz <= ep ) return;
928 }
929
930 dzot=dz*.5;
931 saoa(z+dzot,g3);
932 saoa(z+dz,g5);
933
934 } /* if( ! jump ) */
935
936 nogo=FALSE;
937 for( i = 0; i < n; i++ )
938 {
939 t00=(g1[i]+g5[i])*dzot;
940 t01[i]=(t00+dz*g3[i])*.5;
941 t10[i]=(4.0*t01[i]-t00)/3.0;
942
943 /* test convergence of 3 point romberg result */
944 test( creal(t01[i]), creal(t10[i]), &tr,
945 cimag(t01[i]), cimag(t10[i]), &ti, 0.0 );
946 if( (tr > CRIT) || (ti > CRIT) )
947 nogo = TRUE;
948 }
949
950 if( ! nogo )
951 {
952 for( i = 0; i < n; i++ )
953 sum[i] += t10[i];
954
955 nt += 2;
956 z += dz;
957 if(z > zend) return;
958
959 for( i = 0; i < n; i++ )
960 g1[i]=g5[i];
961
962 if( (nt >= NTS) && (ns > nx) )
963 {
964 ns=ns/2;
965 nt=1;
966 }
967
968 jump = FALSE;
969 continue;
970
971 } /* if( ! nogo ) */
972
973 saoa(z+dz*.250,g2);
974 saoa(z+dz*.75,g4);
975 nogo=FALSE;
976 for( i = 0; i < n; i++ )
977 {
978 t02=(t01[i]+dzot*(g2[i]+g4[i]))*0.50;
979 t11=(4.0*t02-t01[i])/3.0;
980 t20[i]=(16.0*t11-t10[i])/15.0;
981
982 /* test convergence of 5 point romberg result */
983 test( creal(t11), creal(t20[i]), &tr,
984 cimag(t11), cimag(t20[i]), &ti, 0.0 );
985 if( (tr > CRIT) || (ti > CRIT) )
986 nogo = TRUE;
987 }
988
989 if( ! nogo )
990 {
991 for( i = 0; i < n; i++ )
992 sum[i] += t20[i];
993
994 nt++;
995 z += dz;
996 if(z > zend) return;
997
998 for( i = 0; i < n; i++ )
999 g1[i]=g5[i];
1000
1001 if( (nt >= NTS) && (ns > nx) )
1002 {
1003 ns=ns/2;
1004 nt=1;
1005 }
1006
1007 jump = FALSE;
1008 continue;
1009
1010 } /* if( ! nogo ) */
1011
1012 nt=0;
1013 if(ns < NM)
1014 {
1015 ns *= 2;
1016 dz=s/ns;
1017 dzot=dz*0.50;
1018
1019 for( i = 0; i < n; i++ )
1020 {
1021 g5[i]=g3[i];
1022 g3[i]=g2[i];
1023 }
1024
1025 jump = TRUE;
1026 continue;
1027
1028 } /* if(ns < nm) */
1029
1030 if( ! lstep )
1031 {
1032 lstep = TRUE;
1033 lambda( z, &t00, &t11 );
1034 }
1035
1036 for( i = 0; i < n; i++ )
1037 sum[i] += t20[i];
1038
1039 nt++;
1040 z += dz;
1041 if(z > zend) return;
1042
1043 for( i = 0; i < n; i++ )
1044 g1[i]=g5[i];
1045
1046 jump = FALSE;
1047
1048 } /* while( TRUE ) */
1049
1050 }
1051
1052 /*-----------------------------------------------------------------------*/
1053
1054 /* saoa computes the integrand for each of the 6 sommerfeld */
1055 /* integrals for source and observer above ground */
saoa(double t,complex double * ans)1056 void saoa( double t, complex double *ans)
1057 {
1058 double xlr;
1059 static complex double xl, dxl, cgam1, cgam2, b0;
1060 static complex double b0p, com, dgam, den1, den2;
1061
1062 lambda(t, &xl, &dxl);
1063 if( jh == 0 )
1064 {
1065 /* bessel function form */
1066 bessel(xl*rho, &b0, &b0p);
1067 b0 *=2.0;
1068 b0p *=2.0;
1069 cgam1=csqrt(xl*xl-ck1sq);
1070 cgam2=csqrt(xl*xl-ck2sq);
1071 if(creal(cgam1) == 0.0)
1072 cgam1=cmplx(0.0,-fabs(cimag(cgam1)));
1073 if(creal(cgam2) == 0.0)
1074 cgam2=cmplx(0.0,-fabs(cimag(cgam2)));
1075 }
1076 else
1077 {
1078 /* hankel function form */
1079 hankel(xl*rho, &b0, &b0p);
1080 com=xl-ck1;
1081 cgam1=csqrt(xl+ck1)*csqrt(com);
1082 if(creal(com) < 0.0 && cimag(com) >= 0.0)
1083 cgam1=-cgam1;
1084 com=xl-ck2;
1085 cgam2=csqrt(xl+ck2)*csqrt(com);
1086 if(creal(com) < 0.0 && cimag(com) >= 0.0)
1087 cgam2=-cgam2;
1088 }
1089
1090 xlr=creal( xl*conj(xl) );
1091 if(xlr >= tsmag)
1092 {
1093 double sign;
1094 if(cimag(xl) >= 0.0)
1095 {
1096 xlr=creal(xl);
1097 if(xlr >= ck2)
1098 {
1099 if(xlr <= ck1r)
1100 dgam=cgam2-cgam1;
1101 else
1102 {
1103 sign=1.0;
1104 dgam=1.0/(xl*xl);
1105 dgam=sign*((ct3*dgam+ct2)*dgam+ct1)/xl;
1106 }
1107 }
1108 else
1109 {
1110 sign=-1.0;
1111 dgam=1.0/(xl*xl);
1112 dgam=sign*((ct3*dgam+ct2)*dgam+ct1)/xl;
1113 } /* if(xlr >= ck2) */
1114
1115 } /* if(cimag(xl) >= 0.0) */
1116 else
1117 {
1118 sign=1.0;
1119 dgam=1.0/(xl*xl);
1120 dgam=sign*((ct3*dgam+ct2)*dgam+ct1)/xl;
1121 }
1122
1123 } /* if(xlr < tsmag) */
1124 else dgam=cgam2-cgam1;
1125
1126 den2=cksm*dgam/(cgam2*(ck1sq*cgam2+ck2sq*cgam1));
1127 den1=1.0/(cgam1+cgam2)-cksm/cgam2;
1128 com=dxl*xl*cexp(-cgam2*zph);
1129 ans[5]=com*b0*den1/ck1;
1130 com *= den2;
1131
1132 if(rho != 0.0)
1133 {
1134 b0p=b0p/rho;
1135 ans[0]=-com*xl*(b0p+b0*xl);
1136 ans[3]=com*xl*b0p;
1137 }
1138 else
1139 {
1140 ans[0]=-com*xl*xl*0.5;
1141 ans[3]=ans[0];
1142 }
1143
1144 ans[1]=com*cgam2*cgam2*b0;
1145 ans[2]=-ans[3]*cgam2*rho;
1146 ans[4]=com*b0;
1147
1148 return;
1149 }
1150
1151 /*-----------------------------------------------------------------------*/
1152
1153