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