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 /* xnec2c.c
20  *
21  * Contains functions that carry out various
22  * operations that were packed spaggetti-fashion
23  * in the original NEC2 main() function
24  */
25 
26 #include "xnec2c.h"
27 #include "shared.h"
28 
29 /* Left-overs from fortran code :-( */
30 static double tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
31 
32 /*-----------------------------------------------------------------------*/
33 
34 /* Frequency_Scale_Geometry()
35  *
36  * Scales geometric parameters to frequency
37  */
38   void
Frequency_Scale_Geometry()39 Frequency_Scale_Geometry()
40 {
41   double fr;
42   int idx;
43 
44   /* Calculate wavelength */
45   data.wlam= CVEL/ calc_data.fmhz;
46 
47   /* frequency scaling of geometric parameters */
48   fr= calc_data.fmhz / CVEL;
49   if( data.n != 0)
50   {
51 	for( idx = 0; idx < data.n; idx++ )
52 	{
53 	  data.x[idx] = save.xtemp[idx] * fr;
54 	  data.y[idx] = save.ytemp[idx] * fr;
55 	  data.z[idx] = save.ztemp[idx] * fr;
56 	  data.si[idx]= save.sitemp[idx]* fr;
57 	  data.bi[idx]= save.bitemp[idx]* fr;
58 	}
59   }
60 
61   if( data.m != 0)
62   {
63 	double fr2= fr* fr;
64 	for( idx = 0; idx < data.m; idx++ )
65 	{
66 	  int j;
67 
68 	  j = idx + data.n;
69 	  data.px[idx] = save.xtemp[j] * fr;
70 	  data.py[idx] = save.ytemp[j] * fr;
71 	  data.pz[idx] = save.ztemp[j] * fr;
72 	  data.pbi[idx]= save.bitemp[j]* fr2;
73 	}
74   }
75 
76 } /* Frequency_Scale_Geometry() */
77 
78 /*-----------------------------------------------------------------------*/
79 
80 /* Struct_Impedance_Loading()
81  *
82  * Calculates structure (segment) impedance loading
83  */
84   void
Structure_Impedance_Loading(void)85 Structure_Impedance_Loading( void )
86 {
87   /* Calculate some loading parameters */
88   if( zload.nload != 0)
89 	load(
90 		calc_data.ldtyp,  calc_data.ldtag,
91 		calc_data.ldtagf, calc_data.ldtagt,
92 		calc_data.zlr,    calc_data.zli,
93 		calc_data.zlc );
94 
95 } /* Struct_Impedance_Loading() */
96 
97 /*-----------------------------------------------------------------------*/
98 
99 /* Ground_Parameters()
100  *
101  * Calculates ground parameters (antenna environment)
102  */
103   void
Ground_Parameters(void)104 Ground_Parameters( void )
105 {
106   complex double epsc;
107 
108   if( gnd.ksymp != 1)
109   {
110 	gnd.frati = CPLX_10;
111 
112 	if( gnd.iperf != 1)
113 	{
114 	  if( save.sig < 0.0 )
115 		save.sig = -save.sig / (59.96 * data.wlam);
116 
117 	  epsc = cmplx( save.epsr, -save.sig * data.wlam * 59.96 );
118 	  gnd.zrati = 1.0 / csqrt( epsc);
119 	  gwav.u = gnd.zrati;
120 	  gwav.u2 = gwav.u * gwav.u;
121 
122 	  if( gnd.nradl > 0 )
123 	  {
124 		gnd.scrwl = save.scrwlt / data.wlam;
125 		gnd.scrwr = save.scrwrt / data.wlam;
126 		gnd.t1 = CPLX_01 * 2367.067/ (double)gnd.nradl;
127 		gnd.t2 = gnd.scrwr * (double)gnd.nradl;
128 	  } /* if( gnd.nradl > 0 ) */
129 
130 	  if( gnd.iperf == 2)
131 	  {
132 		somnec( save.epsr, save.sig, calc_data.fmhz );
133 		gnd.frati =( epsc - 1.0) / ( epsc + 1.0);
134 		if( cabs(( ggrid.epscf - epsc) / epsc) >= 1.0e-3 )
135 		{
136 		  fprintf( stderr,
137 			  "xnec2c: Ground_Parameters(): error in ground parameters\n"
138 			  "complex dielectric constant from file: %12.5E%+12.5Ej\n"
139 			  "                            requested: %12.5E%+12.5Ej\n",
140 			  creal(ggrid.epscf), cimag(ggrid.epscf),
141 			  creal(epsc), cimag(epsc) );
142 		  stop( _("Ground_Parameters():"\
143 				"Error in ground parameters"), ERR_STOP );
144 		}
145 	  } /* if( gnd.iperf != 2) */
146 	} /* if( gnd.iperf != 1) */
147 	else
148 	{
149 	  gnd.scrwl = 0.0;
150 	  gnd.scrwr = 0.0;
151 	  gnd.t1 = 0.0;
152 	  gnd.t2 = 0.0;
153 	}
154   } /* if( gnd.ksymp != 1) */
155 
156   return;
157 } /* Ground_Parameters() */
158 
159 /*-----------------------------------------------------------------------*/
160 
161 /* Set_Interaction_Matrix()
162  *
163  * Sets and factors the interaction matrix
164  */
165   void
Set_Interaction_Matrix(void)166 Set_Interaction_Matrix( void )
167 {
168   /* Memory allocation for symmetry array */
169   smat.nop = netcx.neq/netcx.npeq;
170   size_t mreq = (size_t)(smat.nop * smat.nop) * sizeof( complex double);
171   mem_realloc( (void **)&smat.ssx, mreq, "in xnec2c.c" );
172 
173   /* irngf is not used (NGF function not implemented) */
174   int iresrv = data.np2m * (data.np + 2 * data.mp);
175   if( matpar.imat == 0)
176 	fblock( netcx.npeq, netcx.neq, iresrv, data.ipsym);
177 
178   cmset( netcx.neq, cm, calc_data.rkh, calc_data.iexk );
179   factrs( netcx.npeq, netcx.neq, cm, save.ip );
180   netcx.ntsol = 0;
181 
182 } /* Set_Interaction_Matrix() */
183 
184 /*-----------------------------------------------------------------------*/
185 
186 /* Set_Excitation()
187  *
188  * Sets the excitation part of the matrix
189  */
190   void
Set_Excitation(void)191 Set_Excitation( void )
192 {
193   if( (fpat.ixtyp >= 1) && (fpat.ixtyp <= 4) )
194   {
195 	tmp4= TA* calc_data.xpr4;
196 	tmp5= TA* calc_data.xpr5;
197 
198 	if( fpat.ixtyp == 4)
199 	{
200 	  tmp1= calc_data.xpr1/ data.wlam;
201 	  tmp2= calc_data.xpr2/ data.wlam;
202 	  tmp3= calc_data.xpr3/ data.wlam;
203 	  tmp6= calc_data.xpr6/( data.wlam* data.wlam);
204 	}
205 	else
206 	{
207 	  tmp1= TA* calc_data.xpr1;
208 	  tmp2= TA* calc_data.xpr2;
209 	  tmp3= TA* calc_data.xpr3;
210 	  tmp6= calc_data.xpr6;
211 	} /* if( fpat.ixtyp == 4) */
212 
213   } /* if( (fpat.ixtyp >= 1) && (fpat.ixtyp <= 4) ) */
214 
215   /* fills e field right-hand matrix */
216   etmns( tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, fpat.ixtyp, crnt.cur );
217 
218 } /* Set_Excitation() */
219 
220 /*-----------------------------------------------------------------------*/
221 
222 /* Set_Network_Data()
223  *
224  * Sets up network data and solves for currents
225  */
226   void
Set_Network_Data(void)227 Set_Network_Data( void )
228 {
229   if( netcx.nonet != 0 )
230   {
231 	int i, j, itmp1, itmp2, itmp3;
232 
233 	itmp3=0;
234 	itmp1= netcx.ntyp[0];
235 	for( i = 0; i < 2; i++ )
236 	{
237 	  if( itmp1 == 3) itmp1=2;
238 
239 	  for( j = 0; j < netcx.nonet; j++)
240 	  {
241 		itmp2= netcx.ntyp[j];
242 
243 		if( (itmp2/itmp1) != 1 ) itmp3 = itmp2;
244 		else if( (itmp2 >= 2) && (netcx.x11i[j] <= 0.0) )
245 		{
246 		  double xx, yy, zz;
247 		  int idx4, idx5;
248 
249 		  idx4 = netcx.iseg1[j]-1;
250 		  idx5 = netcx.iseg2[j]-1;
251 		  xx = data.x[idx5]- data.x[idx4];
252 		  yy = data.y[idx5]- data.y[idx4];
253 		  zz = data.z[idx5]- data.z[idx4];
254 		  netcx.x11i[j] = data.wlam* sqrt( xx*xx + yy*yy + zz*zz );
255 		}
256 
257 	  } /* for( j = 0; j < netcx.nonet; j++) */
258 
259 	  if( itmp3 == 0) break;
260 
261 	  itmp1= itmp3;
262 
263 	} /* for( i = 0; i < 2; i++ ) */
264 
265   } /* if( netcx.nonet != 0 ) */
266 
267   /* Set network data */
268   netwk( cm, save.ip, crnt.cur );
269   netcx.ntsol = 1;
270 
271   /* Save impedance data for normalization */
272   if( ((calc_data.nfrq > 1) && isFlagSet(FREQ_LOOP_RUNNING)) || CHILD )
273   {
274 	impedance_data.zreal[calc_data.fstep] = (double)creal( netcx.zped);
275 	impedance_data.zimag[calc_data.fstep] = (double)cimag( netcx.zped);
276 	impedance_data.zmagn[calc_data.fstep] = (double)cabs( netcx.zped);
277 	impedance_data.zphase[calc_data.fstep]= (double)cang( netcx.zped);
278 
279 	if( (calc_data.iped == 1) &&
280 		((double)impedance_data.zmagn[calc_data.fstep] >
281 		 calc_data.zpnorm) )
282 	  calc_data.zpnorm =
283 		(double)impedance_data.zmagn[calc_data.fstep];
284   }
285 
286 } /* Set_Network_Data() */
287 
288 /*-----------------------------------------------------------------------*/
289 
290 /* Power_Loss()
291  *
292  * Calculate power loss due to segment loading
293  */
294   void
Power_Loss(void)295 Power_Loss( void )
296 {
297   int i;
298   double cmg;
299   complex double curi;
300 
301 
302   /* No wire/segments in structure */
303   if( data.n == 0) return;
304 
305   fpat.ploss = 0.0;
306   /* Loop over all wire segs */
307   for( i = 0; i < data.n; i++ )
308   {
309 	/* Calculate segment current (mag/phase) */
310 	curi= crnt.cur[i]* data.wlam;
311 	cmg= cabs( curi);
312 
313 	/* Calculate power loss in segment */
314 	if( (zload.nload != 0) &&
315 		(fabs(creal(zload.zarray[i])) >= 1.0e-20) )
316 	  fpat.ploss += 0.5* cmg* cmg* creal( zload.zarray[i])* data.si[i];
317 
318   } /* for( i = 0; i < n; i++ ) */
319 
320 } /* Power_Loss() */
321 
322 /*-----------------------------------------------------------------------*/
323 
324 /* Radiation_Pattern()
325  *
326  * Calculates far field (radiation) pattern
327  */
328   void
Radiation_Pattern(void)329 Radiation_Pattern( void )
330 {
331   if( (gnd.ifar != 1) && isFlagSet(ENABLE_RDPAT) )
332   {
333 	fpat.pinr= netcx.pin;
334 	fpat.pnlr= netcx.pnls;
335 	rdpat();
336   }
337 
338 } /* Radiation_Pattern() */
339 
340 /*-----------------------------------------------------------------------*/
341 
342 /* Near_Field_Pattern()
343  *
344  * Calculates near field pattern if enabled/needed
345  */
346   void
Near_Field_Pattern(void)347 Near_Field_Pattern( void )
348 {
349   if( near_field.valid ||
350 	  isFlagClear(DRAW_EHFIELD) ||
351 	  isFlagClear(ENABLE_NEAREH) )
352 	return;
353 
354   if( isFlagSet(DRAW_EFIELD) )
355 	nfpat(0);
356 
357   if( isFlagSet(DRAW_HFIELD) )
358 	nfpat(1);
359 
360 } /* Near_Field_Pattern() */
361 
362 /*-----------------------------------------------------------------------*/
363 
364 /* New_Frequency()
365  *
366  * (Re)calculates all frequency-dependent parameters
367  */
368   void
New_Frequency(void)369 New_Frequency( void )
370 {
371   /* Abort if freq has not really changed, as when changing
372    * between current or charge density structure coloring */
373   if( (save.last_freq == calc_data.fmhz) || isFlagClear(ENABLE_EXCITN) )
374 	return;
375   save.last_freq = calc_data.fmhz;
376 
377   /* Frequency scaling of geometric parameters */
378   Frequency_Scale_Geometry();
379 
380   /* Structure segment loading */
381   Structure_Impedance_Loading();
382 
383   /* Calculate ground parameters */
384   Ground_Parameters();
385 
386   /* Fill and factor primary interaction matrix */
387   Set_Interaction_Matrix();
388 
389   /* Fill excitation part of matrix */
390   Set_Excitation();
391 
392   /* Matrix solving (netwk calls solves) */
393   crnt.valid = 0;
394   Set_Network_Data();
395 
396   /* Calculate power loss */
397   Power_Loss();
398 
399   /* Calculate radiation pattern */
400   Radiation_Pattern();
401 
402   /* Near field calculation */
403   near_field.valid = 0;
404   Near_Field_Pattern();
405 
406 } /* New_Frequency()  */
407 
408 /*-----------------------------------------------------------------------*/
409 
410 static gboolean retval;	/* Function's return value */
411 
412 /* Frequency_Loop()
413  *
414  * Loops over frequency if calculations over a frequency range is
415  * requested, dividing the job between child processes if forked
416  */
417   gboolean
Frequency_Loop(gpointer udata)418 Frequency_Loop( gpointer udata )
419 {
420   /* Value of frequency and step num in the loop */
421   static double freq;
422 
423   /* Current freq step, saved steps
424    * index, num of busy processes */
425   static int fstep, num_busy_procs;
426 
427   int idx, job_num = 0;
428   size_t len;
429   char *buff;				/* Used to pass on structure poiners */
430   fd_set read_fds;			/* Read file descriptors for select() */
431 
432 
433   /* (Re) initialize freq loop */
434   if( isFlagSet(FREQ_LOOP_INIT) )
435   {
436 	/* Clear global flags */
437 	ClearFlag( FREQ_LOOP_INIT |	FREQ_LOOP_DONE );
438 
439 	/* (Re)-enable freq loop (back to start freq) */
440 	freq = save.fmhz;
441 
442 	/* Step back frequency and step count since incrementing
443 	 * is done at start of frequency loop calculations */
444 	fstep = -1;
445 	if( calc_data.ifrq == 1)
446 	  freq /= calc_data.delfrq;
447 	else
448 	  freq -= calc_data.delfrq;
449 
450 	/* Clear list of "valid" (processed) loop steps */
451 	for( idx = 0; idx < calc_data.nfrq; idx++ )
452 	  save.fstep[idx] = 0;
453 
454 	/* Clear "last-used-frequency" buffer */
455 	save.last_freq = 0.0;
456 
457 	/* Zero num of busy processes */
458 	num_busy_procs = 0;
459 
460 	/* Signal global freq step "illegal" */
461 	calc_data.fstep = -1;
462 
463 	/* Inherited from NEC2 */
464 	if( calc_data.zpnorm > 0.0 )
465 	  calc_data.iped = 2;
466 
467 	/* Continue gtk_main idle callbacks */
468 	retval = TRUE;
469 	return ( retval );
470 
471   } /* isFlagSet(INIT_FREQ_LOOP) */
472 
473   /* Repeat freq stepping over number of child processes
474    * if forked. calc_data.num_jobs = 1 for non-forked runs.
475    * If not forked (no multi-threading), following block will
476    * execute only once, since only one instance is running */
477   for( idx = 0; idx < calc_data.num_jobs; idx++ )
478   {
479 	/* Up frequency step count */
480 	fstep++;
481 
482 	/* Frequency loop is completed or was paused by user */
483 	if( (fstep >= calc_data.nfrq) || isFlagSet(FREQ_LOOP_STOP) )
484 	{
485 	  /* Re-initialize if loop completed all steps */
486 	  if( fstep >= calc_data.nfrq )
487 		SetFlag( FREQ_LOOP_INIT );
488 
489 	  /* Points to last buffer in rad_pattern filled by loop */
490 	  fstep--;
491 
492 	  /* Last freq step that was processed by children */
493 	  calc_data.lastf = fstep;
494 
495 	  /* Re-enable pausing of freq loop */
496 	  ClearFlag( FREQ_LOOP_STOP );
497 
498 	  /* Cancel idle callbacks on exit */
499 	  retval = FALSE;
500 
501 	  break;
502 	} /* if( (fstep >= calc_data.nfrq) || isFlagSet(FREQ_LOOP_STOP) ) */
503 
504 	/* Increment frequency */
505 	if( calc_data.ifrq == 1)
506 	  freq *= calc_data.delfrq;	/* Multiplicative stepping */
507 	else
508 	  freq += calc_data.delfrq; /* Additive stepping */
509 
510 	/* Save frequencies for plotting */
511 	save.freq[fstep] = (double)freq;
512 
513 	/* Delegate calculations to child processes if forked */
514 	if( FORKED )
515 	{
516 	  /* Look for an idle process */
517 	  for( job_num = 0; job_num < calc_data.num_jobs; job_num++ )
518 	  {
519 		/* If an idle process is found, give it a job and
520 		 * then step the frequency loop by breaking out */
521 		if( ! forked_proc_data[job_num]->busy )
522 		{
523 		  /* Signal and count busy processes */
524 		  forked_proc_data[job_num]->busy  = TRUE;
525 		  forked_proc_data[job_num]->fstep = fstep;
526 		  num_busy_procs++;
527 
528 		  /* Tell process to calculate freq dependent data */
529 		  len = strlen( fork_commands[FRQDATA] );
530 		  Write_Pipe( job_num, fork_commands[FRQDATA], (ssize_t)len, TRUE );
531 
532 		  /* When it responds, give it next frequency */
533 		  buff = (char *)&freq;
534 		  len = sizeof( double );
535 		  Write_Pipe( job_num, buff, (ssize_t)len, TRUE );
536 		  break;
537 		}
538 	  } /* for( job_num = 0; job_num < calc_data.num_jobs; job_num++ ) */
539 
540 	} /* if( FORKED ) */
541 	else /* Calculate freq dependent data (no fork) */
542 	{
543 	  calc_data.fmhz  = freq;
544 	  calc_data.fstep = fstep;
545 	  calc_data.lastf = fstep;
546 	  New_Frequency();
547 	  break;
548 	}
549 
550 	/* All idle processes are given a job */
551 	if( num_busy_procs >= calc_data.num_jobs )
552 	  break;
553 
554   } /* for( idx = 0; idx < calc_data.num_jobs; idx++ ) */
555 
556   /* Receive results from forked children */
557   if( FORKED && num_busy_procs )
558 	do
559 	{
560 	  int n = 0;
561 
562 	  /* Set read fd's to watch for child writes */
563 	  FD_ZERO( &read_fds );
564 	  for( idx = 0; idx < calc_data.num_jobs; idx++ )
565 	  {
566 		FD_SET( forked_proc_data[idx]->child2pnt_pipe[READ], &read_fds );
567 		if( n < forked_proc_data[idx]->child2pnt_pipe[READ] )
568 		  n = forked_proc_data[idx]->child2pnt_pipe[READ];
569 	  }
570 
571 	  /* Wait for data from finished child processes */
572 	  if( select( n+1, &read_fds, NULL, NULL, NULL ) == -1 )
573 	  {
574 		perror( "xnec2c: select()" );
575 		_exit(0);
576 	  }
577 
578 	  /* Check for finished child processes */
579 	  for( idx = 0; idx < num_child_procs; idx++ )
580 	  {
581 		if( FD_ISSET(forked_proc_data[idx]->child2pnt_pipe[READ], &read_fds) )
582 		{
583 		  /* Read data from finished child process */
584 		  Get_Freq_Data( idx, forked_proc_data[idx]->fstep );
585 
586 		  /* Mark freq step in list of processed steps */
587 		  save.fstep[forked_proc_data[idx]->fstep] = 1;
588 
589 		  /* Mark finished child process as ready for next job */
590 		  forked_proc_data[idx]->busy = FALSE;
591 
592 		  /* Count down number of busy processes */
593 		  num_busy_procs--;
594 		}
595 	  } /* for( idx = 0; idx < num_child_procs; idx++ ) */
596 
597 	  /* Find highest freq step that has no steps below it
598 	   * that have not been processed by a child process */
599 	  for( idx = 0; idx < calc_data.nfrq; idx++ )
600 		if( save.fstep[idx] ) calc_data.fstep = idx;
601 		else break;
602 
603 	} /* do */
604 	/* Loop terminated and busy children */
605 	while( !retval && num_busy_procs );
606 
607   /* Return if freq step 0 not ready yet */
608   if( calc_data.fstep < 0 ) return( retval );
609 
610   /* Set frequency and step to global variables */
611   calc_data.lastf = calc_data.fstep;
612   calc_data.fmhz = (double)save.freq[calc_data.fstep];
613 
614   /* Trigger a redraw of plots drawingarea */
615   Plot_Frequency_Data();
616 
617   /* Set frequency spinbuttons */
618   gtk_spin_button_set_value(
619 	  mainwin_frequency, (gdouble)calc_data.fmhz );
620 
621   if( isFlagSet(DRAW_ENABLED) )
622 	gtk_spin_button_set_value(
623 		rdpattern_frequency, (gdouble)calc_data.fmhz );
624 
625   if( isFlagSet(PLOT_ENABLED) )
626   {
627 	char txt[10];
628 	snprintf( txt, sizeof(txt), "%10.3f", (gdouble)calc_data.fmhz );
629 	gtk_entry_set_text( GTK_ENTRY(
630 		  lookup_widget(freqplots_window, "freqplots_fmhz_entry")), txt );
631   }
632 
633   /* Change flags at exit if loop is done */
634   if( !retval && !num_busy_procs )
635   {
636 	ClearFlag( FREQ_LOOP_RUNNING );
637 	SetFlag( FREQ_LOOP_DONE );
638   }
639 
640   return( retval );
641 
642 } /* Frequency_Loop() */
643 
644 /*-----------------------------------------------------------------------*/
645 
646 /* Start_Frequency_Loop()
647  *
648  * Starts frequency loop
649  */
650   gboolean
Start_Frequency_Loop(void)651 Start_Frequency_Loop( void )
652 {
653   if( isFlagClear(FREQ_LOOP_RUNNING) && (calc_data.nfrq > 1) )
654   {
655 	retval = TRUE;
656 	SetFlag(FREQ_LOOP_RUNNING);
657 	floop_tag = g_idle_add( Frequency_Loop, NULL );
658 	return( TRUE );
659   }
660   else return( FALSE );
661 
662 } /* Start_Frequency_Loop() */
663 
664 /*-----------------------------------------------------------------------*/
665 
666 /* Stop_Frequency_Loop()
667  *
668  * Stops and resets freq loop
669  */
670   void
Stop_Frequency_Loop(void)671 Stop_Frequency_Loop( void )
672 {
673   if( floop_tag > 0 )
674   {
675 	g_source_remove( floop_tag );
676 	floop_tag = 0;
677   }
678   ClearFlag( FREQ_LOOP_RUNNING );
679 
680 } /* Stop_Frequency_Loop() */
681 
682 /*-----------------------------------------------------------------------*/
683 
684 /* Incident_Field_Loop()
685  *
686  * Loops over incident field directions if
687  * receiving pattern calculations are requested
688  */
689   void
Incident_Field_Loop(void)690 Incident_Field_Loop( void )
691 {
692   int phi_step, theta_step;
693 
694   /* Frequency scaling of geometric parameters */
695   Frequency_Scale_Geometry();
696 
697   /* Structure segment loading */
698   Structure_Impedance_Loading();
699 
700   /* Calculate ground parameters */
701   Ground_Parameters();
702 
703   /* Fill and factor primary interaction matrix */
704   Set_Interaction_Matrix();
705 
706   /* Loop over incident field angles */
707   netcx.nprint=0;
708   /* Loop over phi */
709   for( phi_step = 0; phi_step < calc_data.nphi; phi_step++ )
710   {
711 	/* Loop over theta */
712 	for( theta_step = 0; theta_step < calc_data.nthi; theta_step++ )
713 	{
714 	  /* Fill excitation part of matrix */
715 	  Set_Excitation();
716 
717 	  /* Matrix solving (netwk calls solves) */
718 	  Set_Network_Data();
719 
720 	  /* Calculate power loss */
721 	  Power_Loss();
722 
723 	  calc_data.xpr1 += calc_data.xpr4;
724 
725 	} /* for( theta_step = 0; theta_step < calc_data.nthi.. */
726 
727 	calc_data.xpr1= calc_data.thetis;
728 	calc_data.xpr2= calc_data.xpr2+ calc_data.xpr5;
729 
730   } /* for( phi_step = 0; phi_step < calc_data.nphi.. */
731 
732   calc_data.xpr2  = calc_data.phiss;
733 
734 } /* Incident_Field_Loop() */
735 
736 /*-----------------------------------------------------------------------*/
737 
738