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 "ground.h"
45 #include "shared.h"
46
47 /*-------------------------------------------------------------------*/
48
49 /* segment to obtain the total field due to ground. the method of */
50 /* variable interval width romberg integration is used. there are 9 */
51 /* field components - the x, y, and z components due to constant, */
52 /* sine, and cosine current distributions. */
53 void
rom2(double a,double b,complex double * sum,double dmin)54 rom2( double a, double b,
55 complex double *sum, double dmin )
56 {
57 int i, ns, nt, flag = TRUE;
58 int nts = 4, nx = 1, n = 9;
59 double ze, ep, zend, dz=0.0, dzot=0.0, tmag1, tmag2, tr, ti;
60 double z, s; /***also global***/
61 double rx = 1.0e-4;
62 complex double t00, t02, t11;
63 static complex double *g1 = NULL, *g2 = NULL;
64 static complex double *g3 = NULL, *g4 = NULL;
65 static complex double *t01 = NULL, *t10 = NULL;
66 static complex double *t20 = NULL, *g5 = NULL;
67 static gboolean first_call = TRUE;
68
69 if( first_call )
70 {
71 first_call = FALSE;
72 size_t mreq = 9 * sizeof(complex double);
73 mem_alloc( (void **)&g1, mreq, "in ground.c");
74 mem_alloc( (void **)&g2, mreq, "in ground.c");
75 mem_alloc( (void **)&g3, mreq, "in ground.c");
76 mem_alloc( (void **)&g4, mreq, "in ground.c");
77 mem_alloc( (void **)&g5, mreq, "in ground.c");
78 mem_alloc( (void **)&t01, mreq, "in ground.c");
79 mem_alloc( (void **)&t10, mreq, "in ground.c");
80 mem_alloc( (void **)&t20, mreq, "in ground.c");
81 }
82
83 z= a;
84 ze= b;
85 s= b- a;
86
87 if( s < 0.0)
88 {
89 fprintf( stderr, "xnec2c: b less than a in rom2\n" );
90 stop( _("rom2(): b less than a"), ERR_STOP );
91 }
92
93 ep= s/(1.0e4* data.npm);
94 zend= ze- ep;
95
96 for( i = 0; i < n; i++ )
97 sum[i]=CPLX_00;
98
99 ns= nx;
100 nt=0;
101 sflds( z, g1);
102
103 while( TRUE )
104 {
105 if( flag )
106 {
107 dz= s/ ns;
108 if( z+ dz > ze)
109 {
110 dz= ze- z;
111 if( dz <= ep) return;
112 }
113
114 dzot= dz*.5;
115 sflds( z+ dzot, g3);
116 sflds( z+ dz, g5);
117
118 } /* if( flag ) */
119
120 tmag1=0.0;
121 tmag2=0.0;
122
123 /* evaluate 3 point romberg result and test convergence. */
124 for( i = 0; i < n; i++ )
125 {
126 t00=( g1[i]+ g5[i])* dzot;
127 t01[i]=( t00+ dz* g3[i])*.5;
128 t10[i]=(4.0* t01[i]- t00)/3.0;
129 if( i > 2)
130 continue;
131
132 tr= creal( t01[i]);
133 ti= cimag( t01[i]);
134 tmag1= tmag1+ tr* tr+ ti* ti;
135 tr= creal( t10[i]);
136 ti= cimag( t10[i]);
137 tmag2= tmag2+ tr* tr+ ti* ti;
138
139 } /* for( i = 0; i < n; i++ ) */
140
141 tmag1= sqrt( tmag1);
142 tmag2= sqrt( tmag2);
143 test( tmag1, tmag2, &tr, 0.0, 0.0, &ti, dmin);
144
145 if( tr <= rx)
146 {
147 for( i = 0; i < n; i++ )
148 sum[i] += t10[i];
149 nt += 2;
150
151 z += dz;
152 if( z > zend)
153 return;
154
155 for( i = 0; i < n; i++ )
156 g1[i]= g5[i];
157
158 if( (nt >= nts) && (ns > nx) )
159 {
160 ns= ns/2;
161 nt=1;
162 }
163 flag = TRUE;
164 continue;
165
166 } /* if( tr <= rx) */
167
168 sflds( z+ dz*.25, g2);
169 sflds( z+ dz*.75, g4);
170 tmag1=0.0;
171 tmag2=0.0;
172
173 /* evaluate 5 point romberg result and test convergence. */
174 for( i = 0; i < n; i++ )
175 {
176 t02=( t01[i]+ dzot*( g2[i]+ g4[i]))*.5;
177 t11=( 4.0 * t02- t01[i] )/3.0;
178 t20[i]=(16.0* t11- t10[i])/15.0;
179 if( i > 2)
180 continue;
181
182 tr= creal( t11);
183 ti= cimag( t11);
184 tmag1= tmag1+ tr* tr+ ti* ti;
185 tr= creal( t20[i]);
186 ti= cimag( t20[i]);
187 tmag2= tmag2+ tr* tr+ ti* ti;
188
189 } /* for( i = 0; i < n; i++ ) */
190
191 tmag1= sqrt( tmag1);
192 tmag2= sqrt( tmag2);
193 test( tmag1, tmag2, &tr, 0.0,0.0, &ti, dmin);
194
195 if( tr > rx)
196 {
197 nt=0;
198 if( ns < data.npm )
199 {
200 ns= ns*2;
201 dz= s/ ns;
202 dzot= dz*.5;
203
204 for( i = 0; i < n; i++ )
205 {
206 g5[i]= g3[i];
207 g3[i]= g2[i];
208 }
209
210 flag=FALSE;
211 continue;
212
213 } /* if( ns < npm) */
214
215 fprintf( stderr,
216 "xnec2c: rom2 -- step size limited at z = %12.5E\n", z );
217
218 } /* if( tr > rx) */
219
220 for( i = 0; i < n; i++ )
221 sum[i]= sum[i]+ t20[i];
222 nt= nt+1;
223
224 z= z+ dz;
225 if( z > zend)
226 return;
227
228 for( i = 0; i < n; i++ )
229 g1[i]= g5[i];
230
231 flag = TRUE;
232 if( (nt < nts) || (ns <= nx) )
233 continue;
234
235 ns= ns/2;
236 nt=1;
237
238 } /* while( TRUE ) */
239
240 }
241
242 /*-----------------------------------------------------------------------*/
243
244 /* sfldx returns the field due to ground for a current element on */
245 /* the source segment at t relative to the segment center. */
246 void
sflds(double t,complex double * e)247 sflds( double t, complex double *e )
248 {
249 double xt, yt, zt, rhx, rhy, rhs, rho, phx, phy;
250 double cph, sph, zphs, r2s, rk, sfac, thet;
251 complex double eph, erv, ezv, erh, ezh;
252 complex double er, et, hrv, hzv, hrh;
253
254 xt= dataj.xj+ t* dataj.cabj;
255 yt= dataj.yj+ t* dataj.sabj;
256 zt= dataj.zj+ t* dataj.salpj;
257 rhx= incom.xo- xt;
258 rhy= incom.yo- yt;
259 rhs= rhx* rhx+ rhy* rhy;
260 rho= sqrt( rhs);
261
262 if( rho <= 0.0)
263 {
264 rhx=1.0;
265 rhy=0.0;
266 phx=0.0;
267 phy=1.0;
268 }
269 else
270 {
271 rhx= rhx/ rho;
272 rhy= rhy/ rho;
273 phx= -rhy;
274 phy= rhx;
275 }
276
277 cph= rhx* incom.xsn+ rhy* incom.ysn;
278 sph= rhy* incom.xsn- rhx* incom.ysn;
279
280 if( fabs( cph) < 1.0e-10)
281 cph=0.0;
282 if( fabs( sph) < 1.0e-10)
283 sph=0.0;
284
285 gwav.zph= incom.zo+ zt;
286 zphs= gwav.zph* gwav.zph;
287 r2s= rhs+ zphs;
288 gwav.r2= sqrt( r2s);
289 rk= gwav.r2* TP;
290 gwav.xx2= cmplx( cos( rk),- sin( rk));
291
292 /* use norton approximation for field due to ground.
293 * current is lumped at segment center with current moment
294 * for constant, sine or cosine distribution. */
295 if( incom.isnor != 1)
296 {
297 gwav.zmh=1.0;
298 gwav.r1=1.0;
299 gwav.xx1=0.0;
300 gwave( &erv, &ezv, &erh, &ezh, &eph);
301
302 et=-CONST1* gnd.frati* gwav.xx2/( r2s* gwav.r2);
303 er=2.0* et* cmplx(1.0, rk);
304 et= et* cmplx(1.0 - rk* rk, rk);
305 hrv=( er+ et)* rho* gwav.zph/ r2s;
306 hzv=( zphs* er- rhs* et)/ r2s;
307 hrh=( rhs* er- zphs* et)/ r2s;
308 erv= erv- hrv;
309 ezv= ezv- hzv;
310 erh= erh+ hrh;
311 ezh= ezh+ hrv;
312 eph= eph+ et;
313 erv= erv* dataj.salpj;
314 ezv= ezv* dataj.salpj;
315 erh= erh* incom.sn* cph;
316 ezh= ezh* incom.sn* cph;
317 eph= eph* incom.sn* sph;
318 erh= erv+ erh;
319 e[0]=( erh* rhx+ eph* phx)* dataj.s;
320 e[1]=( erh* rhy+ eph* phy)* dataj.s;
321 e[2]=( ezv+ ezh)* dataj.s;
322 e[3]=0.0;
323 e[4]=0.0;
324 e[5]=0.0;
325 sfac= PI* dataj.s;
326 sfac= sin( sfac)/ sfac;
327 e[6]= e[0]* sfac;
328 e[7]= e[1]* sfac;
329 e[8]= e[2]* sfac;
330
331 return;
332 } /* if( smat.isnor != 1) */
333
334 /* interpolate in sommerfeld field tables */
335 if( rho >= 1.0e-12)
336 thet= atan( gwav.zph/ rho);
337 else
338 thet= POT;
339
340 /* combine vertical and horizontal components and convert */
341 /* to x,y,z components. multiply by exp(-jkr)/r. */
342 intrp( gwav.r2, thet, &erv, &ezv, &erh, &eph );
343 gwav.xx2= gwav.xx2/ gwav.r2;
344 sfac= incom.sn* cph;
345 erh= gwav.xx2*( dataj.salpj* erv+ sfac* erh);
346 ezh= gwav.xx2*( dataj.salpj* ezv- sfac* erv);
347 /* x,y,z fields for constant current */
348 eph= incom.sn* sph* gwav.xx2* eph;
349 e[0]= erh* rhx+ eph* phx;
350 e[1]= erh* rhy+ eph* phy;
351 e[2]= ezh;
352 /* x,y,z fields for sine current */
353 rk= TP* t;
354 sfac= sin( rk);
355 e[3]= e[0]* sfac;
356 e[4]= e[1]* sfac;
357 /* x,y,z fields for cosine current */
358 e[5]= e[2]* sfac;
359 sfac= cos( rk);
360 e[6]= e[0]* sfac;
361 e[7]= e[1]* sfac;
362 e[8]= e[2]* sfac;
363
364 return;
365 }
366
367 /*-----------------------------------------------------------------------*/
368
369