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 "calculations.h"
45 #include "shared.h"
46
47 /*-------------------------------------------------------------------*/
48
49 /* fill incident field array for charge discontinuity voltage source */
qdsrc(int is,complex double v,complex double * e)50 void qdsrc( int is, complex double v, complex double *e )
51 {
52 int i, jx, j, jp1, ipr, ij, i1;
53 double xi, yi, zi, ai, cabi, sabi, salpi, tx, ty, tz;
54 complex double curd, etk, ets, etc;
55
56 is--;
57 i= data.icon1[is];
58 data.icon1[is]=0;
59 tbf( is+1,0);
60 data.icon1[is]= i;
61 dataj.s= data.si[is]*.5;
62 curd= CCJ* v/(( log(2.0 * dataj.s/ data.bi[is])-1.0) *
63 ( segj.bx[segj.jsno-1] * cos( TP* dataj.s) +
64 segj.cx[segj.jsno-1] * sin( TP* dataj.s))* data.wlam);
65 vsorc.vqds[vsorc.nqds]= v;
66 vsorc.iqds[vsorc.nqds]= is+1;
67 vsorc.nqds++;
68
69 for( jx = 0; jx < segj.jsno; jx++ )
70 {
71 j= segj.jco[jx]-1;
72 jp1 = j+1;
73 dataj.s= data.si[j];
74 dataj.b= data.bi[j];
75 dataj.xj= data.x[j];
76 dataj.yj= data.y[j];
77 dataj.zj= data.z[j];
78 dataj.cabj= data.cab[j];
79 dataj.sabj= data.sab[j];
80 dataj.salpj= data.salp[j];
81
82 if( dataj.iexk != 0)
83 {
84 ipr= data.icon1[j];
85
86 if (ipr > PCHCON) dataj.ind1=2;
87 else if( ipr < 0 )
88 {
89 ipr= -ipr;
90 ipr--;
91 if( -data.icon1[ipr-1] != jp1 )
92 dataj.ind1=2;
93 else
94 {
95 xi= fabs( dataj.cabj* data.cab[ipr]+ dataj.sabj*
96 data.sab[ipr]+ dataj.salpj* data.salp[ipr]);
97 if( (xi < 0.999999) ||
98 (fabs(data.bi[ipr]/dataj.b-1.0) > 1.0e-6) )
99 dataj.ind1=2;
100 else
101 dataj.ind1=0;
102 }
103 } /* if( ipr < 0 ) */
104 else if( ipr == 0 )
105 dataj.ind1=1;
106 else /* ipr > 0 */
107 {
108 ipr--;
109 if( ipr != j )
110 {
111 if( data.icon2[ipr] != jp1)
112 dataj.ind1=2;
113 else
114 {
115 xi= fabs( dataj.cabj* data.cab[ipr]+ dataj.sabj*
116 data.sab[ipr]+ dataj.salpj* data.salp[ipr]);
117 if( (xi < 0.999999) ||
118 (fabs(data.bi[ipr]/dataj.b-1.0) > 1.0e-6) )
119 dataj.ind1=2;
120 else
121 dataj.ind1=0;
122 }
123 } /* if( ipr != j ) */
124 else
125 {
126 if( (dataj.cabj*dataj.cabj + dataj.sabj*dataj.sabj) > 1.0e-8)
127 dataj.ind1=2;
128 else
129 dataj.ind1=0;
130 }
131 } /* else */
132
133 ipr= data.icon2[j];
134 if (ipr > PCHCON) dataj.ind2=2;
135 else if( ipr < 0 )
136 {
137 ipr = -ipr;
138 ipr--;
139 if( -data.icon2[ipr] != jp1 )
140 dataj.ind1=2;
141 else
142 {
143 xi= fabs( dataj.cabj* data.cab[ipr]+ dataj.sabj *
144 data.sab[ipr]+ dataj.salpj* data.salp[ipr]);
145 if( (xi < 0.999999) ||
146 (fabs(data.bi[ipr]/dataj.b-1.0) > 1.0e-6) )
147 dataj.ind1=2;
148 else
149 dataj.ind1=0;
150 }
151 } /* if( ipr < 0 ) */
152 else if( ipr == 0 )
153 dataj.ind2=1;
154 else /* ipr > 0 */
155 {
156 ipr--;
157 if( ipr != j )
158 {
159 if( data.icon1[ipr] != jp1)
160 dataj.ind2=2;
161 else
162 {
163 xi= fabs( dataj.cabj* data.cab[ipr]+ dataj.sabj*
164 data.sab[ipr]+ dataj.salpj* data.salp[ipr]);
165 if( (xi < 0.9999990) ||
166 (fabs(data.bi[ipr]/dataj.b-1.0) > 1.0e-6) )
167 dataj.ind2=2;
168 else
169 dataj.ind2=0;
170 }
171 } /* if( ipr != j )*/
172 else
173 {
174 if( (dataj.cabj* dataj.cabj + dataj.sabj* dataj.sabj) > 1.0e-8)
175 dataj.ind1=2;
176 else
177 dataj.ind1=0;
178 }
179 } /* else */
180
181 } /* if( dataj.iexk != 0) */
182
183 for( i = 0; i < data.n; i++ )
184 {
185 ij= i- j;
186 xi= data.x[i];
187 yi= data.y[i];
188 zi= data.z[i];
189 ai= data.bi[i];
190 efld( xi, yi, zi, ai, ij);
191 cabi= data.cab[i];
192 sabi= data.sab[i];
193 salpi= data.salp[i];
194 etk= dataj.exk* cabi+ dataj.eyk* sabi+ dataj.ezk* salpi;
195 ets= dataj.exs* cabi+ dataj.eys* sabi+ dataj.ezs* salpi;
196 etc= dataj.exc* cabi+ dataj.eyc* sabi+ dataj.ezc* salpi;
197 e[i]= e[i]-( etk* segj.ax[jx] +
198 ets* segj.bx[jx]+ etc* segj.cx[jx])* curd;
199 }
200
201 if( data.m != 0)
202 {
203 i1= data.n-1;
204 for( i = 0; i < data.m; i++ )
205 {
206 xi= data.px[i];
207 yi= data.py[i];
208 zi= data.pz[i];
209 hsfld( xi, yi, zi, 0.0);
210 i1++;
211 tx= data.t2x[i];
212 ty= data.t2y[i];
213 tz= data.t2z[i];
214 etk= dataj.exk* tx+ dataj.eyk* ty+ dataj.ezk* tz;
215 ets= dataj.exs* tx+ dataj.eys* ty+ dataj.ezs* tz;
216 etc= dataj.exc* tx+ dataj.eyc* ty+ dataj.ezc* tz;
217 e[i1] += ( etk* segj.ax[jx]+ ets* segj.bx[jx]+
218 etc* segj.cx[jx] )* curd* data.psalp[i];
219 i1++;
220 tx= data.t1x[i];
221 ty= data.t1y[i];
222 tz= data.t1z[i];
223 etk= dataj.exk* tx+ dataj.eyk* ty+ dataj.ezk* tz;
224 ets= dataj.exs* tx+ dataj.eys* ty+ dataj.ezs* tz;
225 etc= dataj.exc* tx+ dataj.eyc* ty+ dataj.ezc* tz;
226 e[i1] += ( etk* segj.ax[jx]+ ets* segj.bx[jx]+
227 etc* segj.cx[jx])* curd* data.psalp[i];
228 }
229
230 } /* if( m != 0) */
231
232 if( zload.nload > 0 )
233 e[j] += zload.zarray[j]* curd*(segj.ax[jx]+ segj.cx[jx]);
234
235 } /* for( jx = 0; jx < segj.jsno; jx++ ) */
236
237 return;
238 }
239
240 /*-----------------------------------------------------------------------*/
241
242 /* cabc computes coefficients of the constant (a), sine (b), and */
243 /* cosine (c) terms in the current interpolation functions for the */
244 /* current vector cur. */
cabc(complex double * curx)245 void cabc( complex double *curx)
246 {
247 int i, is, j, jx, jco1, jco2;
248 double ar, ai, sh;
249 complex double curd, cs1, cs2;
250
251 if( data.n != 0)
252 {
253 for( i = 0; i < data.n; i++ )
254 {
255 crnt.air[i]=0.0;
256 crnt.aii[i]=0.0;
257 crnt.bir[i]=0.0;
258 crnt.bii[i]=0.0;
259 crnt.cir[i]=0.0;
260 crnt.cii[i]=0.0;
261 }
262
263 for( i = 0; i < data.n; i++ )
264 {
265 ar= creal( curx[i]);
266 ai= cimag( curx[i]);
267 tbf( i+1, 1 );
268
269 for( jx = 0; jx < segj.jsno; jx++ )
270 {
271 j= segj.jco[jx]-1;
272 crnt.air[j] += segj.ax[jx]* ar;
273 crnt.aii[j] += segj.ax[jx]* ai;
274 crnt.bir[j] += segj.bx[jx]* ar;
275 crnt.bii[j] += segj.bx[jx]* ai;
276 crnt.cir[j] += segj.cx[jx]* ar;
277 crnt.cii[j] += segj.cx[jx]* ai;
278 }
279
280 } /* for( i = 0; i < n; i++ ) */
281
282 if( vsorc.nqds != 0)
283 {
284 for( is = 0; is < vsorc.nqds; is++ )
285 {
286 i= vsorc.iqds[is]-1;
287 jx= data.icon1[i];
288 data.icon1[i]=0;
289 tbf(i+1,0);
290 data.icon1[i]= jx;
291 sh= data.si[i]*.5;
292 curd= CCJ* vsorc.vqds[is]/( (log(2.0* sh/ data.bi[i])-1.0) *
293 (segj.bx[segj.jsno-1]* cos(TP* sh)+ segj.cx[segj.jsno-1] *
294 sin(TP* sh))* data.wlam );
295 ar= creal( curd);
296 ai= cimag( curd);
297
298 for( jx = 0; jx < segj.jsno; jx++ )
299 {
300 j= segj.jco[jx]-1;
301 crnt.air[j]= crnt.air[j]+ segj.ax[jx]* ar;
302 crnt.aii[j]= crnt.aii[j]+ segj.ax[jx]* ai;
303 crnt.bir[j]= crnt.bir[j]+ segj.bx[jx]* ar;
304 crnt.bii[j]= crnt.bii[j]+ segj.bx[jx]* ai;
305 crnt.cir[j]= crnt.cir[j]+ segj.cx[jx]* ar;
306 crnt.cii[j]= crnt.cii[j]+ segj.cx[jx]* ai;
307 }
308
309 } /* for( is = 0; is < vsorc.nqds; is++ ) */
310
311 } /* if( vsorc.nqds != 0) */
312
313 for( i = 0; i < data.n; i++ )
314 curx[i]= cmplx( crnt.air[i]+crnt.cir[i], crnt.aii[i]+crnt.cii[i] );
315
316 } /* if( n != 0) */
317
318 if( data.m == 0)
319 return;
320
321 /* convert surface currents from */
322 /* t1,t2 components to x,y,z components */
323 jco1= data.np2m;
324 jco2= jco1+ data.m;
325 for( i = 1; i <= data.m; i++ )
326 {
327 jco1 -= 2;
328 jco2 -= 3;
329 cs1= curx[jco1];
330 cs2= curx[jco1+1];
331 curx[jco2] = cs1* data.t1x[data.m-i]+ cs2* data.t2x[data.m-i];
332 curx[jco2+1]= cs1* data.t1y[data.m-i]+ cs2* data.t2y[data.m-i];
333 curx[jco2+2]= cs1* data.t1z[data.m-i]+ cs2* data.t2z[data.m-i];
334 }
335
336 return;
337 }
338
339 /*-----------------------------------------------------------------------*/
340
341 /* function db10 returns db for magnitude (field) */
db10(double x)342 double db10( double x )
343 {
344 if( x < 1.0e-20 )
345 return( -999.99 );
346
347 return( 10.0 * log10(x) );
348 }
349
350 /*-----------------------------------------------------------------------*/
351
352 /* function db20 returns db for mag**2 (power) i */
db20(double x)353 double db20( double x )
354 {
355 if( x < 1.0e-20 )
356 return( -999.99 );
357
358 return( 20.0 * log10(x) );
359 }
360
361 /*-----------------------------------------------------------------------*/
362
363 /* intrp uses bivariate cubic interpolation to obtain */
364 /* the values of 4 functions at the point (x,y). */
intrp(double x,double y,complex double * f1,complex double * f2,complex double * f3,complex double * f4)365 void intrp( double x, double y, complex double *f1,
366 complex double *f2, complex double *f3, complex double *f4 )
367 {
368 static int ix, iy, ixs=-10, iys=-10, igrs=-10, ixeg=0, iyeg=0;
369 static int nxm2, nym2, nxms, nyms, nd, ndp;
370 static int *nda = NULL, *ndpa = NULL;
371 int jump;
372 static double dx = 1.0, dy = 1.0, xs = 0.0, ys = 0.0, xz, yz;
373 double xx, yy;
374 static complex double a[4][4], b[4][4], c[4][4], d[4][4];
375 complex double p1=CPLX_00, p2=CPLX_00, p3=CPLX_00, p4=CPLX_00;
376 complex double fx1, fx2, fx3, fx4;
377 static gboolean first_call = TRUE;
378
379 if( first_call )
380 {
381 first_call = FALSE;
382 size_t mreq = 3*sizeof(int);
383 mem_alloc( (void **)&nda, mreq, "in calculations.c");
384 mem_alloc( (void **)&ndpa, mreq, "in calculations.c");
385 nda[0] = 11; nda[1] = 17; nda[2] = 9;
386 ndpa[0] = 110; ndpa[1] = 85; ndpa[2] = 72;
387 }
388
389 jump = FALSE;
390 if( (x < xs) || (y < ys) )
391 jump = TRUE;
392 else
393 {
394 ix= (int)(( x- xs)/ dx)+1;
395 iy= (int)(( y- ys)/ dy)+1;
396 }
397
398 /* if point lies in same 4 by 4 point region */
399 /* as previous point, old values are reused. */
400 if( (ix < ixeg) ||
401 (iy < iyeg) ||
402 (abs(ix- ixs) >= 2) ||
403 (abs(iy- iys) >= 2) ||
404 jump )
405 {
406 int igr, iadd, iadz, i, k;
407
408 /* determine correct grid and grid region */
409 if( x <= ggrid.xsa[1])
410 igr=0;
411 else
412 {
413 if( y > ggrid.ysa[2])
414 igr=2;
415 else
416 igr=1;
417 }
418
419 if( igr != igrs)
420 {
421 igrs= igr;
422 dx= ggrid.dxa[igrs];
423 dy= ggrid.dya[igrs];
424 xs= ggrid.xsa[igrs];
425 ys= ggrid.ysa[igrs];
426 nxm2= ggrid.nxa[igrs]-2;
427 nym2= ggrid.nya[igrs]-2;
428 nxms=(( nxm2+1)/3)*3+1;
429 nyms=(( nym2+1)/3)*3+1;
430 nd= nda[igrs];
431 ndp= ndpa[igrs];
432 ix= (int)(( x- xs)/ dx)+1;
433 iy= (int)(( y- ys)/ dy)+1;
434
435 } /* if( igr != igrs) */
436
437 ixs=(( ix-1)/3)*3+2;
438 if( ixs < 2)
439 ixs=2;
440 ixeg=-10000;
441
442 if( ixs > nxm2)
443 {
444 ixs= nxm2;
445 ixeg= nxms;
446 }
447
448 iys=(( iy-1)/3)*3+2;
449 if( iys < 2)
450 iys=2;
451 iyeg=-10000;
452
453 if( iys > nym2)
454 {
455 iys= nym2;
456 iyeg= nyms;
457 }
458
459 /* compute coefficients of 4 cubic polynomials in x for */
460 /* the 4 grid values of y for each of the 4 functions */
461 iadz= ixs+( iys-3)* nd- ndp;
462 for( k = 0; k < 4; k++ )
463 {
464 iadz += ndp;
465 iadd = iadz;
466
467 for( i = 0; i < 4; i++ )
468 {
469 iadd += nd;
470
471 switch( igrs )
472 {
473 case 0:
474 p1= ggrid.ar1[iadd-2];
475 p2= ggrid.ar1[iadd-1];
476 p3= ggrid.ar1[iadd];
477 p4= ggrid.ar1[iadd+1];
478 break;
479
480 case 1:
481 p1= ggrid.ar2[iadd-2];
482 p2= ggrid.ar2[iadd-1];
483 p3= ggrid.ar2[iadd];
484 p4= ggrid.ar2[iadd+1];
485 break;
486
487 case 2:
488 p1= ggrid.ar3[iadd-2];
489 p2= ggrid.ar3[iadd-1];
490 p3= ggrid.ar3[iadd];
491 p4= ggrid.ar3[iadd+1];
492
493 } /* switch( igrs ) */
494
495 a[i][k]=( p4- p1+3.0*( p2- p3))*.1666666667;
496 b[i][k]=( p1-2.0* p2+ p3)*.5;
497 c[i][k]= p3-(2.0* p1+3.0* p2+ p4)*.1666666667;
498 d[i][k]= p2;
499
500 } /* for( i = 0; i < 4; i++ ) */
501
502 } /* for( k = 0; k < 4; k++ ) */
503
504 xz=( ixs-1)* dx+ xs;
505 yz=( iys-1)* dy+ ys;
506
507 } /* if( (abs(ix- ixs) >= 2) || */
508
509 /* evaluate polymomials in x and use cubic */
510 /* interpolation in y for each of the 4 functions. */
511 xx=( x- xz)/ dx;
512 yy=( y- yz)/ dy;
513 fx1=(( a[0][0]* xx+ b[0][0])* xx+ c[0][0])* xx+ d[0][0];
514 fx2=(( a[1][0]* xx+ b[1][0])* xx+ c[1][0])* xx+ d[1][0];
515 fx3=(( a[2][0]* xx+ b[2][0])* xx+ c[2][0])* xx+ d[2][0];
516 fx4=(( a[3][0]* xx+ b[3][0])* xx+ c[3][0])* xx+ d[3][0];
517 p1= fx4- fx1+3.0*( fx2- fx3);
518 p2=3.0*( fx1-2.0* fx2+ fx3);
519 p3=6.0* fx3-2.0* fx1-3.0* fx2- fx4;
520 *f1=(( p1* yy+ p2)* yy+ p3)* yy*.1666666667+ fx2;
521 fx1=(( a[0][1]* xx+ b[0][1])* xx+ c[0][1])* xx+ d[0][1];
522 fx2=(( a[1][1]* xx+ b[1][1])* xx+ c[1][1])* xx+ d[1][1];
523 fx3=(( a[2][1]* xx+ b[2][1])* xx+ c[2][1])* xx+ d[2][1];
524 fx4=(( a[3][1]* xx+ b[3][1])* xx+ c[3][1])* xx+ d[3][1];
525 p1= fx4- fx1+3.0*( fx2- fx3);
526 p2=3.0*( fx1-2.0* fx2+ fx3);
527 p3=6.0* fx3-2.0* fx1-3.0* fx2- fx4;
528 *f2=(( p1* yy+ p2)* yy+ p3)* yy*.1666666667+ fx2;
529 fx1=(( a[0][2]* xx+ b[0][2])* xx+ c[0][2])* xx+ d[0][2];
530 fx2=(( a[1][2]* xx+ b[1][2])* xx+ c[1][2])* xx+ d[1][2];
531 fx3=(( a[2][2]* xx+ b[2][2])* xx+ c[2][2])* xx+ d[2][2];
532 fx4=(( a[3][2]* xx+ b[3][2])* xx+ c[3][2])* xx+ d[3][2];
533 p1= fx4- fx1+3.0*( fx2- fx3);
534 p2=3.0*( fx1-2.0* fx2+ fx3);
535 p3=6.0* fx3-2.0* fx1-3.0* fx2- fx4;
536 *f3=(( p1* yy+ p2)* yy+ p3)* yy*.1666666667+ fx2;
537 fx1=(( a[0][3]* xx+ b[0][3])* xx+ c[0][3])* xx+ d[0][3];
538 fx2=(( a[1][3]* xx+ b[1][3])* xx+ c[1][3])* xx+ d[1][3];
539 fx3=(( a[2][3]* xx+ b[2][3])* xx+ c[2][3])* xx+ d[2][3];
540 fx4=(( a[3][3]* xx+ b[3][3])* xx+ c[3][3])* xx+ d[3][3];
541 p1= fx4- fx1+3.0*( fx2- fx3);
542 p2=3.0*( fx1-2.0* fx2+ fx3);
543 p3=6.0* fx3-2.0* fx1-3.0* fx2- fx4;
544 *f4=(( p1* yy+ p2)* yy+ p3)* yy*.1666666667+ fx2;
545
546 return;
547 }
548
549 /*-----------------------------------------------------------------------*/
550
551 /* intx performs numerical integration of exp(jkr)/r by the method of */
552 /* variable interval width romberg integration. the integrand value */
553 /* is supplied by subroutine gf. */
intx(double el1,double el2,double b,int ij,double * sgr,double * sgi)554 void intx( double el1, double el2, double b,
555 int ij, double *sgr, double *sgi )
556 {
557 int ns, nt;
558 int nx = 1, nma = 65536, nts = 4;
559 int flag = TRUE;
560 double z, s, ze, fnm, ep, zend, fns, dz=0.0, zp;
561 double t00r, g1r, g5r=0.0, t00i, g1i, g5i=0.0, t01r, g3r=0.0;
562 double t01i, g3i=0.0, t10r, t10i, te1i, te1r, t02r;
563 double g2r, g4r, t02i, g2i, g4i, t11r, t11i, t20r;
564 double t20i, te2i, te2r, rx = 1.0e-4, dzot=0.0;
565
566 z= el1;
567 ze= el2;
568 if( ij == 0)
569 ze=0.0;
570 s= ze- z;
571 fnm= nma;
572 ep= s/(10.0* fnm);
573 zend= ze- ep;
574 *sgr=0.0;
575 *sgi=0.0;
576 ns= nx;
577 nt=0;
578 gf( z, &g1r, &g1i);
579
580 while( TRUE )
581 {
582 if( flag )
583 {
584 fns= ns;
585 dz= s/ fns;
586 zp= z+ dz;
587
588 if( zp > ze)
589 {
590 dz= ze- z;
591 if( fabs(dz) <= ep)
592 {
593 /* add contribution of near singularity for diagonal term */
594 if(ij == 0)
595 {
596 *sgr=2.0*( *sgr+ log(( sqrt( b* b+ s* s)+ s)/ b));
597 *sgi=2.0* *sgi;
598 }
599 return;
600 }
601
602 } /* if( zp > ze) */
603
604 dzot= dz*.5;
605 zp= z+ dzot;
606 gf( zp, &g3r, &g3i);
607 zp= z+ dz;
608 gf( zp, &g5r, &g5i);
609
610 } /* if( flag ) */
611
612 t00r=( g1r+ g5r)* dzot;
613 t00i=( g1i+ g5i)* dzot;
614 t01r=( t00r+ dz* g3r)*0.5;
615 t01i=( t00i+ dz* g3i)*0.5;
616 t10r=(4.0* t01r- t00r)/3.0;
617 t10i=(4.0* t01i- t00i)/3.0;
618
619 /* test convergence of 3 point romberg result. */
620 test( t01r, t10r, &te1r, t01i, t10i, &te1i, 0.0);
621 if( (te1i <= rx) && (te1r <= rx) )
622 {
623 *sgr= *sgr+ t10r;
624 *sgi= *sgi+ t10i;
625 nt += 2;
626
627 z += dz;
628 if( z >= zend)
629 {
630 /* add contribution of near singularity for diagonal term */
631 if(ij == 0)
632 {
633 *sgr=2.0*( *sgr+ log(( sqrt( b* b+ s* s)+ s)/ b));
634 *sgi=2.0* *sgi;
635 }
636 return;
637 }
638
639 g1r= g5r;
640 g1i= g5i;
641 if( nt >= nts)
642 if( ns > nx)
643 {
644 /* Double step size */
645 ns= ns/2;
646 nt=1;
647 }
648 flag = TRUE;
649 continue;
650
651 } /* if( (te1i <= rx) && (te1r <= rx) ) */
652
653 zp= z+ dz*0.25;
654 gf( zp, &g2r, &g2i);
655 zp= z+ dz*0.75;
656 gf( zp, &g4r, &g4i);
657 t02r=( t01r+ dzot*( g2r+ g4r))*0.5;
658 t02i=( t01i+ dzot*( g2i+ g4i))*0.5;
659 t11r=(4.0* t02r- t01r)/3.0;
660 t11i=(4.0* t02i- t01i)/3.0;
661 t20r=(16.0* t11r- t10r)/15.0;
662 t20i=(16.0* t11i- t10i)/15.0;
663
664 /* test convergence of 5 point romberg result. */
665 test( t11r, t20r, &te2r, t11i, t20i, &te2i, 0.0);
666 if( (te2i > rx) || (te2r > rx) )
667 {
668 nt=0;
669 if( ns >= nma)
670 fprintf( stderr,
671 "xnec2c: step size limited at z= %10.5f\n", z );
672 else
673 {
674 /* halve step size */
675 ns= ns*2;
676 fns= ns;
677 dz= s/ fns;
678 dzot= dz*0.5;
679 g5r= g3r;
680 g5i= g3i;
681 g3r= g2r;
682 g3i= g2i;
683
684 flag = FALSE;
685 continue;
686 }
687
688 } /* if( (te2i > rx) || (te2r > rx) ) */
689
690 *sgr= *sgr+ t20r;
691 *sgi= *sgi+ t20i;
692 nt++;
693
694 z += dz;
695 if( z >= zend)
696 {
697 /* add contribution of near singularity for diagonal term */
698 if(ij == 0)
699 {
700 *sgr=2.0*( *sgr+ log(( sqrt( b* b+ s* s)+ s)/ b));
701 *sgi=2.0* *sgi;
702 }
703 return;
704 }
705
706 g1r= g5r;
707 g1i= g5i;
708 if( nt >= nts)
709 if( ns > nx)
710 {
711 /* Double step size */
712 ns= ns/2;
713 nt=1;
714 }
715 flag = TRUE;
716
717 } /* while( TRUE ) */
718
719 }
720
721 /*-----------------------------------------------------------------------*/
722
723 /* returns smallest of two arguments */
min(int a,int b)724 int min( int a, int b )
725 {
726 if( a < b )
727 return(a);
728 else
729 return(b);
730 }
731
732 /*-----------------------------------------------------------------------*/
733
734 /* test for convergence in numerical integration */
test(double f1r,double f2r,double * tr,double f1i,double f2i,double * ti,double dmin)735 void test( double f1r, double f2r, double *tr,
736 double f1i, double f2i, double *ti, double dmin )
737 {
738 double den;
739
740 den= fabs( f2r);
741 *tr= fabs( f2i);
742
743 if( den < *tr)
744 den= *tr;
745 if( den < dmin)
746 den= dmin;
747
748 if( den < 1.0e-37)
749 {
750 *tr=0.0;
751 *ti=0.0;
752 return;
753 }
754
755 *tr= fabs(( f1r- f2r)/ den);
756 *ti= fabs(( f1i- f2i)/ den);
757
758 return;
759 }
760
761 /*-----------------------------------------------------------------------*/
762
763 /* compute component of basis function i on segment is. */
764 void
sbf(int i,int is,double * aa,double * bb,double * cc)765 sbf( int i, int is, double *aa, double *bb, double *cc )
766 {
767 int ix, jsno, june, jcox, jcoxx, jend, iend, njun1=0, njun2;
768 double d, sig, pp, sdh, cdh, sd, omc;
769 double aj, pm=0, cd, ap, qp, qm, xxi;
770
771 *aa=0.0;
772 *bb=0.0;
773 *cc=0.0;
774 june=0;
775 jsno=0;
776 pp=0.0;
777 ix=i-1;
778
779 jcox= data.icon1[ix];
780 if( jcox > PCHCON)
781 jcox= i;
782
783 jend=-1;
784 iend=-1;
785 sig=-1.0;
786
787 do
788 {
789 if( jcox != 0 )
790 {
791 if( jcox < 0 )
792 jcox= -jcox;
793 else
794 {
795 sig= -sig;
796 jend= -jend;
797 }
798
799 jcoxx = jcox-1;
800 jsno++;
801 d= PI* data.si[jcoxx];
802 sdh= sin( d);
803 cdh= cos( d);
804 sd=2.0* sdh* cdh;
805
806 if( d <= 0.015)
807 {
808 omc=4.0* d* d;
809 omc=((1.3888889e-3* omc -4.1666666667e-2)* omc +.5)* omc;
810 }
811 else omc=1.0- cdh* cdh+ sdh* sdh;
812
813 aj=1.0/( log(1.0/( PI* data.bi[jcoxx]))-.577215664);
814 pp -= omc/ sd* aj;
815
816 if( jcox == is)
817 {
818 *aa= aj/ sd* sig;
819 *bb= aj/(2.0* cdh);
820 *cc= -aj/(2.0* sdh)* sig;
821 june= iend;
822 }
823
824 if( jcox != i )
825 {
826 if( jend != 1)
827 jcox= data.icon1[jcoxx];
828 else
829 jcox= data.icon2[jcoxx];
830
831 if( abs(jcox) != i )
832 {
833 if( jcox == 0 )
834 {
835 fprintf( stderr,
836 "xnec2c: sbf(): segment connection error for segment %d\n", i );
837 stop( _("Segment connection error in sbf()"), ERR_STOP );
838 }
839 else continue;
840 }
841
842 } /* if( jcox != i ) */
843 else if( jcox == is)
844 *bb= -*bb;
845
846 if( iend == 1) break;
847
848 } /* if( jcox != 0 ) */
849
850 pm= -pp;
851 pp=0.0;
852 njun1= jsno;
853
854 jcox= data.icon2[ix];
855 if( jcox > PCHCON)
856 jcox= i;
857
858 jend=1;
859 iend=1;
860 sig=-1.0;
861
862 } /* do */
863 while( jcox != 0 );
864
865 njun2= jsno- njun1;
866 d= PI* data.si[ix];
867 sdh= sin( d);
868 cdh= cos( d);
869 sd=2.0* sdh* cdh;
870 cd= cdh* cdh- sdh* sdh;
871
872 if( d <= 0.015)
873 {
874 omc=4.0* d* d;
875 omc=((1.3888889e-3* omc -4.1666666667e-2)* omc +.5)* omc;
876 }
877 else omc=1.0- cd;
878
879 ap=1.0/( log(1.0/( PI* data.bi[ix])) -.577215664);
880 aj= ap;
881
882 if( njun1 == 0)
883 {
884 if( njun2 == 0)
885 {
886 *aa =-1.0;
887 qp= PI* data.bi[ix];
888 xxi= qp* qp;
889 xxi= qp*(1.0-.5* xxi)/(1.0- xxi);
890 *cc=1.0/( cdh- xxi* sdh);
891 return;
892 }
893
894 qp= PI* data.bi[ix];
895 xxi= qp* qp;
896 xxi= qp*(1.0-.5* xxi)/(1.0- xxi);
897 qp=-( omc+ xxi* sd)/( sd*( ap+ xxi* pp)+ cd*( xxi* ap- pp));
898
899 if( june == 1)
900 {
901 *aa= -*aa* qp;
902 *bb= *bb* qp;
903 *cc= -*cc* qp;
904 if( i != is)
905 return;
906 }
907
908 *aa -= 1.0;
909 d = cd - xxi * sd;
910 *bb += (sdh + ap * qp * (cdh - xxi * sdh)) / d;
911 *cc += (cdh + ap * qp * (sdh + xxi * cdh)) / d;
912 return;
913
914 } /* if( njun1 == 0) */
915
916 if( njun2 == 0)
917 {
918 qm= PI* data.bi[ix];
919 xxi= qm* qm;
920 xxi= qm*(1.0-.5* xxi)/(1.0- xxi);
921 qm=( omc+ xxi* sd)/( sd*( aj- xxi* pm)+ cd*( pm+ xxi* aj));
922
923 if( june == -1)
924 {
925 *aa= *aa* qm;
926 *bb= *bb* qm;
927 *cc= *cc* qm;
928 if( i != is)
929 return;
930 }
931
932 *aa -= 1.0;
933 d= cd- xxi* sd;
934 *bb += ( aj* qm*( cdh- xxi* sdh)- sdh)/ d;
935 *cc += ( cdh- aj* qm*( sdh+ xxi* cdh))/ d;
936 return;
937
938 } /* if( njun2 == 0) */
939
940 qp= sd*( pm* pp+ aj* ap)+ cd*( pm* ap- pp* aj);
941 qm=( ap* omc- pp* sd)/ qp;
942 qp=-( aj* omc+ pm* sd)/ qp;
943
944 if( june != 0 )
945 {
946 if( june < 0 )
947 {
948 *aa= *aa* qm;
949 *bb= *bb* qm;
950 *cc= *cc* qm;
951 }
952 else
953 {
954 *aa= -*aa* qp;
955 *bb= *bb* qp;
956 *cc= -*cc* qp;
957 }
958
959 if( i != is)
960 return;
961
962 } /* if( june != 0 ) */
963
964 *aa -= 1.0;
965 *bb += ( aj* qm+ ap* qp)* sdh/ sd;
966 *cc += ( aj* qm- ap* qp)* cdh/ sd;
967
968 return;
969 }
970
971 /*-----------------------------------------------------------------------*/
972
973 /* compute basis function i */
974 void
tbf(int i,int icap)975 tbf( int i, int icap )
976 {
977 int ix, jcox, jcoxx, jend, iend, njun1=0, njun2, jsnop, jsnox;
978 double pp, sdh, cdh, sd, omc, aj, pm=0, cd, ap, qp, qm, xxi;
979 double d, sig; /*** also global ***/
980
981 segj.jsno=0;
982 pp=0.0;
983 ix = i-1;
984 jcox= data.icon1[ix];
985
986 if( jcox > PCHCON)
987 jcox= i;
988
989 jend=-1;
990 iend=-1;
991 sig=-1.0;
992
993 do
994 {
995 if( jcox != 0 )
996 {
997 if( jcox < 0 )
998 jcox= -jcox;
999 else
1000 {
1001 sig= -sig;
1002 jend= -jend;
1003 }
1004
1005 jcoxx = jcox-1;
1006 segj.jsno++;
1007 jsnox = segj.jsno-1;
1008 segj.jco[jsnox]= jcox;
1009 d= PI* data.si[jcoxx];
1010 sdh= sin( d);
1011 cdh= cos( d);
1012 sd=2.0* sdh* cdh;
1013
1014 if( d <= 0.015)
1015 {
1016 omc=4.0* d* d;
1017 omc=((1.3888889e-3* omc-4.1666666667e-2)* omc+.5)* omc;
1018 }
1019 else omc=1.0- cdh* cdh+ sdh* sdh;
1020
1021 aj=1.0/( log(1.0/( PI* data.bi[jcoxx]))-.577215664);
1022 pp= pp- omc/ sd* aj;
1023 segj.ax[jsnox]= aj/ sd* sig;
1024 segj.bx[jsnox]= aj/(2.0* cdh);
1025 segj.cx[jsnox]= -aj/(2.0* sdh)* sig;
1026
1027 if( jcox != i)
1028 {
1029 if( jend == 1)
1030 jcox= data.icon2[jcoxx];
1031 else
1032 jcox= data.icon1[jcoxx];
1033
1034 if( abs(jcox) != i )
1035 {
1036 if( jcox != 0 )
1037 continue;
1038 else
1039 {
1040 fprintf( stderr,
1041 "xnec2c: tbf(): segment connection error for segment %5d\n", i );
1042 stop( _("Segment connection error in tbf()"), ERR_STOP );
1043 }
1044 }
1045
1046 } /* if( jcox != i) */
1047 else segj.bx[jsnox] = -segj.bx[jsnox];
1048
1049 if( iend == 1) break;
1050
1051 } /* if( jcox != 0 ) */
1052
1053 pm= -pp;
1054 pp=0.0;
1055 njun1= segj.jsno;
1056
1057 jcox= data.icon2[ix];
1058 if( jcox > PCHCON)
1059 jcox= i;
1060
1061 jend=1;
1062 iend=1;
1063 sig=-1.0;
1064
1065 } /* do */
1066 while( jcox != 0 );
1067
1068 njun2= segj.jsno- njun1;
1069 jsnop= segj.jsno;
1070 segj.jco[jsnop]= i;
1071 d= PI* data.si[ix];
1072 sdh= sin( d);
1073 cdh= cos( d);
1074 sd=2.0* sdh* cdh;
1075 cd= cdh* cdh- sdh* sdh;
1076
1077 if( d <= 0.015)
1078 {
1079 omc=4.0* d* d;
1080 omc=((1.3888889e-3* omc-4.1666666667e-2)* omc+.5)* omc;
1081 }
1082 else omc=1.0- cd;
1083
1084 ap=1.0/( log(1.0/( PI* data.bi[ix]))-.577215664);
1085 aj= ap;
1086
1087 if( njun1 == 0)
1088 {
1089 if( njun2 == 0)
1090 {
1091 segj.bx[jsnop]=0.0;
1092
1093 if( icap == 0)
1094 xxi=0.0;
1095 else
1096 {
1097 qp= PI* data.bi[ix];
1098 xxi= qp* qp;
1099 xxi= qp*(1.0-.5* xxi)/(1.0- xxi);
1100 }
1101
1102 segj.cx[jsnop]=1.0/( cdh- xxi* sdh);
1103 segj.jsno= jsnop+1;
1104 segj.ax[jsnop]=-1.0;
1105 return;
1106
1107 } /* if( njun2 == 0) */
1108
1109 if( icap == 0) xxi=0.0;
1110 else
1111 {
1112 qp= PI* data.bi[ix];
1113 xxi= qp* qp;
1114 xxi= qp*(1.0-.5* xxi)/(1.0- xxi);
1115 }
1116
1117 qp=-( omc+ xxi* sd)/( sd*( ap+ xxi* pp)+ cd*( xxi* ap- pp));
1118 d= cd- xxi* sd;
1119 segj.bx[jsnop]=( sdh+ ap* qp*( cdh- xxi* sdh))/ d;
1120 segj.cx[jsnop]=( cdh+ ap* qp*( sdh+ xxi* cdh))/ d;
1121
1122 for( iend = 0; iend < njun2; iend++ )
1123 {
1124 segj.ax[iend]= -segj.ax[iend]* qp;
1125 segj.bx[iend]= segj.bx[iend]* qp;
1126 segj.cx[iend]= -segj.cx[iend]* qp;
1127 }
1128
1129 segj.jsno= jsnop+1;
1130 segj.ax[jsnop]=-1.0;
1131 return;
1132
1133 } /* if( njun1 == 0) */
1134
1135 if( njun2 == 0)
1136 {
1137 if( icap == 0)
1138 xxi=0.0;
1139 else
1140 {
1141 qm= PI* data.bi[ix];
1142 xxi= qm* qm;
1143 xxi= qm*(1.0-.5* xxi)/(1.0- xxi);
1144 }
1145
1146 qm=( omc+ xxi* sd)/( sd*( aj- xxi* pm)+ cd*( pm+ xxi* aj));
1147 d= cd- xxi* sd;
1148 segj.bx[jsnop]=( aj* qm*( cdh- xxi* sdh)- sdh)/ d;
1149 segj.cx[jsnop]=( cdh- aj* qm*( sdh+ xxi* cdh))/ d;
1150
1151 for( iend = 0; iend < njun1; iend++ )
1152 {
1153 segj.ax[iend]= segj.ax[iend]* qm;
1154 segj.bx[iend]= segj.bx[iend]* qm;
1155 segj.cx[iend]= segj.cx[iend]* qm;
1156 }
1157
1158 segj.jsno= jsnop+1;
1159 segj.ax[jsnop]=-1.0;
1160 return;
1161
1162 } /* if( njun2 == 0) */
1163
1164 qp= sd*( pm* pp+ aj* ap)+ cd*( pm* ap- pp* aj);
1165 qm=( ap* omc- pp* sd)/ qp;
1166 qp=-( aj* omc+ pm* sd)/ qp;
1167 segj.bx[jsnop]=( aj* qm+ ap* qp)* sdh/ sd;
1168 segj.cx[jsnop]=( aj* qm- ap* qp)* cdh/ sd;
1169
1170 for( iend = 0; iend < njun1; iend++ )
1171 {
1172 segj.ax[iend]= segj.ax[iend]* qm;
1173 segj.bx[iend]= segj.bx[iend]* qm;
1174 segj.cx[iend]= segj.cx[iend]* qm;
1175 }
1176
1177 jend= njun1;
1178 for( iend = jend; iend < segj.jsno; iend++ )
1179 {
1180 segj.ax[iend]= -segj.ax[iend]* qp;
1181 segj.bx[iend]= segj.bx[iend]* qp;
1182 segj.cx[iend]= -segj.cx[iend]* qp;
1183 }
1184
1185 segj.jsno= jsnop+1;
1186 segj.ax[jsnop]=-1.0;
1187
1188 return;
1189 }
1190
1191 /*-----------------------------------------------------------------------*/
1192
1193 /* compute the components of all basis functions on segment j */
1194 void
trio(int j)1195 trio( int j )
1196 {
1197 int jcox, jcoxx, jsnox, jx, jend=0, iend=0;
1198
1199 segj.jsno=0;
1200 jx = j-1;
1201 jcox= data.icon1[jx];
1202
1203 if( jcox <= PCHCON)
1204 {
1205 jend=-1;
1206 iend=-1;
1207 }
1208
1209 if( (jcox == 0) || (jcox > PCHCON) )
1210 {
1211 jcox= data.icon2[jx];
1212
1213 if( jcox <= PCHCON)
1214 {
1215 jend=1;
1216 iend=1;
1217 }
1218
1219 if( jcox == 0 || (jcox > PCHCON) )
1220 {
1221 jsnox = segj.jsno;
1222 segj.jsno++;
1223
1224 /* Allocate to connections buffers */
1225 if( segj.jsno >= segj.maxcon )
1226 {
1227 segj.maxcon = segj.jsno +1;
1228 size_t mreq = (size_t)segj.maxcon * sizeof(int);
1229 mem_realloc( (void **)&segj.jco, mreq, "in calculations.c" );
1230 mreq = (size_t)segj.maxcon * sizeof(double);
1231 mem_realloc( (void **) &segj.ax, mreq, "in calculations.c" );
1232 mem_realloc( (void **) &segj.bx, mreq, "in calculations.c" );
1233 mem_realloc( (void **) &segj.cx, mreq, "in calculations.c" );
1234 }
1235
1236 sbf( j, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
1237 segj.jco[jsnox]= j;
1238 return;
1239 }
1240
1241 } /* if( (jcox == 0) || (jcox > PCHCON) ) */
1242
1243 do
1244 {
1245 if( jcox < 0 )
1246 jcox= -jcox;
1247 else
1248 jend= -jend;
1249 jcoxx = jcox-1;
1250
1251 if( jcox != j)
1252 {
1253 jsnox = segj.jsno;
1254 segj.jsno++;
1255
1256 /* Allocate to connections buffers */
1257 if( segj.jsno >= segj.maxcon )
1258 {
1259 segj.maxcon = segj.jsno +1;
1260 size_t mreq = (size_t)segj.maxcon * sizeof(int);
1261 mem_realloc( (void **)&segj.jco, mreq, "in calculations.c" );
1262 mreq = (size_t)segj.maxcon * sizeof(double);
1263 mem_realloc( (void **) &segj.ax, mreq, "in calculations.c" );
1264 mem_realloc( (void **) &segj.bx, mreq, "in calculations.c" );
1265 mem_realloc( (void **) &segj.cx, mreq, "in calculations.c" );
1266 }
1267
1268 sbf( jcox, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
1269 segj.jco[jsnox]= jcox;
1270
1271 if( jend != 1)
1272 jcox= data.icon1[jcoxx];
1273 else
1274 jcox= data.icon2[jcoxx];
1275
1276 if( jcox == 0 )
1277 {
1278 fprintf( stderr,
1279 "xnec2c: trio(): segment connention error for segment %5d\n", j );
1280 stop( _("Segment connention error in trio()"), ERR_STOP );
1281 }
1282 else continue;
1283
1284 } /* if( jcox != j) */
1285
1286 if( iend == 1)
1287 break;
1288
1289 jcox= data.icon2[jx];
1290
1291 if( jcox > PCHCON)
1292 break;
1293
1294 jend=1;
1295 iend=1;
1296
1297 } /* do */
1298 while( jcox != 0 );
1299
1300 jsnox = segj.jsno;
1301 segj.jsno++;
1302
1303 /* Allocate to connections buffers */
1304 if( segj.jsno >= segj.maxcon )
1305 {
1306 segj.maxcon = segj.jsno +1;
1307 size_t mreq = (size_t)segj.maxcon * sizeof(int);
1308 mem_realloc( (void **)&segj.jco, mreq, "in calculations.c" );
1309 mreq = (size_t)segj.maxcon * sizeof(double);
1310 mem_realloc( (void **) &segj.ax, mreq, "in calculations.c" );
1311 mem_realloc( (void **) &segj.bx, mreq, "in calculations.c" );
1312 mem_realloc( (void **) &segj.cx, mreq, "in calculations.c" );
1313 }
1314
1315 sbf( j, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
1316 segj.jco[jsnox]= j;
1317
1318 return;
1319 }
1320
1321 /*-----------------------------------------------------------------------*/
1322
1323 /* cang returns the phase angle of a complex number in degrees. */
cang(complex double z)1324 double cang( complex double z )
1325 {
1326 return( carg(z)*TD );
1327 }
1328
1329 /*-----------------------------------------------------------------------*/
1330
1331 /* zint computes the internal impedance of a circular wire */
1332 void
zint(double sigl,double rolam,complex double * zint)1333 zint( double sigl, double rolam, complex double *zint )
1334 {
1335 #define cc1 ( 6.0e-7 + I*1.9e-6)
1336 #define cc2 (-3.4e-6 + I*5.1e-6)
1337 #define cc3 (-2.52e-5 + I*0.0)
1338 #define cc4 (-9.06e-5 - I*9.01e-5)
1339 #define cc5 ( 0.0 - I*9.765e-4)
1340 #define cc6 (.0110486 - I*0.0110485)
1341 #define cc7 ( 0.0 - I*0.3926991)
1342 #define cc8 ( 1.6e-6 - I*3.2e-6)
1343 #define cc9 ( 1.17e-5 - I*2.4e-6)
1344 #define cc10 ( 3.46e-5 + I*3.38e-5)
1345 #define cc11 ( 5.0e-7 + I*2.452e-4)
1346 #define cc12 (-1.3813e-3 + I*1.3811e-3)
1347 #define cc13 (-6.25001e-2 - I*1.0e-7)
1348 #define cc14 (.7071068 + I*0.7071068)
1349 #define cn cc14
1350
1351 #define th(d) ( (((((cc1*(d)+cc2)*(d)+cc3)*(d)+cc4)*(d)+cc5)*(d)+cc6)*(d) + cc7 )
1352 #define ph(d) ( (((((cc8*(d)+cc9)*(d)+cc10)*(d)+cc11)*(d)+cc12)*(d)+cc13)*(d)+cc14 )
1353 #define f(d) ( csqrt(POT/(d))*cexp(-cn*(d)+th(-8.0/x)) )
1354 #define g(d) ( cexp(cn*(d)+th(8.0/x))/csqrt(TP*(d)) )
1355
1356 double x;
1357 double tpcmu = 2.368705e+3;
1358 double cmotp = 60.0;
1359 complex double br1, br2;
1360
1361 x= sqrt( tpcmu* sigl)* rolam;
1362 if( x <= 110.0)
1363 {
1364 if( x <= 8.0)
1365 {
1366 double y, s, ber, bei;
1367
1368 y= x/8.0;
1369 y= y* y;
1370 s= y* y;
1371
1372 ber=((((((-9.01e-6* s+1.22552e-3)* s-.08349609)* s+ 2.6419140)*
1373 s-32.363456)* s+113.77778)* s-64.0)* s+1.0;
1374
1375 bei=((((((1.1346e-4* s-.01103667)* s+.52185615)* s-10.567658)*
1376 s+72.817777)* s-113.77778)* s+16.0)* y;
1377
1378 br1= cmplx( ber, bei);
1379
1380 ber=(((((((-3.94e-6* s+4.5957e-4)* s-.02609253)* s+ .66047849)*
1381 s-6.0681481)* s+14.222222)* s-4.0)* y)* x;
1382
1383 bei=((((((4.609e-5* s-3.79386e-3)* s+.14677204)* s- 2.3116751)*
1384 s+11.377778)* s-10.666667)* s+.5)* x;
1385
1386 br2= cmplx( ber, bei);
1387 br1= br1/ br2;
1388 *zint= CPLX_01* sqrt( cmotp/sigl )* br1/ rolam;
1389
1390 } /* if( x <= 8.0) */
1391 else
1392 {
1393 br2= CPLX_01* f(x)/ PI;
1394 br1= g( x)+ br2;
1395 br2= g( x)* ph(8.0/ x)- br2* ph(-8.0/ x);
1396 br1= br1/ br2;
1397 *zint= CPLX_01* sqrt( cmotp/ sigl)* br1/ rolam;
1398 }
1399
1400 } /* if( x <= 110.0) */
1401 else
1402 {
1403 br1= cmplx(.70710678,-.70710678);
1404 *zint= CPLX_01* sqrt( cmotp/ sigl)* br1/ rolam;
1405 }
1406 }
1407
1408 /*-----------------------------------------------------------------------*/
1409