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 "geometry.h"
45 #include "shared.h"
46
47 /*-------------------------------------------------------------------*/
48
49 /* arc generates segment geometry data for an arc of ns segments */
50 gboolean
arc(int itg,int ns,double rada,double ang1,double ang2,double rad)51 arc( int itg, int ns, double rada,
52 double ang1, double ang2, double rad )
53 {
54 int ist;
55
56 ist= data.n;
57 data.n += ns;
58 data.np= data.n;
59 data.mp= data.m;
60 data.ipsym=0;
61
62 if( ns < 1)
63 {
64 fprintf( stderr,
65 "xnec2c: arc(): number of segments less than 1 (ns < 1)\n");
66 stop( _("arc(): Number of segments < 1"), ERR_OK );
67 return( FALSE );
68 }
69
70 if( fabs( ang2- ang1) < 360.00001)
71 {
72 int i;
73 double ang, dang, xs1, xs2, zs1, zs2;
74 size_t mreq;
75
76 /* Reallocate tags buffer */
77 mreq = (size_t)(data.n + data.m) * sizeof(int);
78 mem_realloc( (void **)&data.itag, mreq, "in geometry.c" );
79
80 /* Reallocate wire buffers */
81 mreq = (size_t)data.n * sizeof(double);
82 mem_realloc( (void **)&data.x1, mreq, "in geometry.c" );
83 mem_realloc( (void **)&data.y1, mreq, "in geometry.c" );
84 mem_realloc( (void **)&data.z1, mreq, "in geometry.c" );
85 mem_realloc( (void **)&data.x2, mreq, "in geometry.c" );
86 mem_realloc( (void **)&data.y2, mreq, "in geometry.c" );
87 mem_realloc( (void **)&data.z2, mreq, "in geometry.c" );
88 mem_realloc( (void **)&data.bi, mreq, "in geometry.c" );
89
90 ang= ang1* TA;
91 dang=( ang2- ang1)* TA/ ns;
92 xs1= rada* cos( ang);
93 zs1= rada* sin( ang);
94
95 for( i = ist; i < data.n; i++ )
96 {
97 ang += dang;
98 xs2= rada* cos( ang);
99 zs2= rada* sin( ang);
100 data.x1[i]= xs1;
101
102 data.y1[i]=0.0;
103 data.z1[i]= zs1;
104 data.x2[i]= xs2;
105 data.y2[i]=0.0;
106 data.z2[i]= zs2;
107 xs1= xs2;
108 zs1= zs2;
109 data.bi[i]= rad;
110 data.itag[i]= itg;
111
112 } /* for( i = ist; i < data.n; i++ ) */
113
114 } /* if( fabs( ang2- ang1) < 360.00001) */
115 else
116 {
117 fprintf( stderr,
118 "xnec2c: arc(): arc angle exceeds 360 degrees\n");
119 stop( _("arc(): Arc angle exceeds 360 degrees"), ERR_OK );
120 return( FALSE );
121 }
122
123 return( TRUE );
124 }
125
126 /*-----------------------------------------------------------------------*/
127
128 /* connect sets up segment connection data in arrays icon1 and */
129 /* icon2 by searching for segment ends that are in contact. */
130 gboolean
conect(int ignd)131 conect( int ignd )
132 {
133 int i, iz, ic, j, jx, ix, ixx, iseg, iend, jend, jump;
134 double sep=0.0, xi1, yi1, zi1, xi2, yi2, zi2;
135 double slen, xa, ya, za, xs, ys, zs;
136 size_t mreq;
137
138 segj.maxcon = 1;
139
140 if( ignd != 0)
141 {
142 if( data.ipsym == 2)
143 {
144 data.np=2* data.np;
145 data.mp=2* data.mp;
146 }
147
148 if( abs(data.ipsym) > 2 )
149 {
150 data.np= data.n;
151 data.mp= data.m;
152 }
153
154 /*** possibly should be error condition?? **/
155 if( data.np > data.n)
156 {
157 fprintf( stderr, "xnec2c: conect(): np > n\n" );
158 stop( _("Error in conect(): np > n"), ERR_OK );
159 return( FALSE );
160 }
161
162 if( (data.np == data.n) && (data.mp == data.m) )
163 data.ipsym=0;
164
165 } /* if( ignd != 0) */
166
167 if( data.n != 0)
168 {
169 /* Allocate memory to connections */
170 mreq = (size_t)(data.n + data.m) * sizeof(int);
171 mem_realloc( (void **)&data.icon1, mreq, "in geometry.c" );
172 mem_realloc( (void **)&data.icon2, mreq, "in geometry.c" );
173
174 for( i = 0; i < data.n; i++ )
175 {
176 data.icon1[i] = data.icon2[i] = 0;
177 iz = i+1;
178 xi1= data.x1[i];
179 yi1= data.y1[i];
180 zi1= data.z1[i];
181 xi2= data.x2[i];
182 yi2= data.y2[i];
183 zi2= data.z2[i];
184 slen= sqrt( (xi2- xi1)*(xi2- xi1) + (yi2- yi1) *
185 (yi2- yi1) + (zi2- zi1)*(zi2- zi1) ) * SMIN;
186
187 /* determine connection data for end 1 of segment. */
188 jump = FALSE;
189 if( ignd > 0)
190 {
191 if( zi1 <= -slen)
192 {
193 fprintf( stderr,
194 "xnec2c: conect(): geometry data error\n"
195 "Segment %d extends below ground\n", iz );
196 stop( _("conect(): Geometry data error\n"\
197 "Segment extends below ground"), ERR_OK );
198 return( FALSE );
199 }
200
201 if( zi1 <= slen)
202 {
203 data.icon1[i]= iz;
204 data.z1[i]=0.0;
205 jump = TRUE;
206
207 } /* if( zi1 <= slen) */
208
209 } /* if( ignd > 0) */
210
211 if( ! jump )
212 {
213 ic= i;
214 for( j = 1; j < data.n; j++)
215 {
216 ic++;
217 if( ic >= data.n)
218 ic=0;
219
220 sep= fabs( xi1- data.x1[ic]) +
221 fabs(yi1- data.y1[ic])+ fabs(zi1- data.z1[ic]);
222 if( sep <= slen)
223 {
224 data.icon1[i]= -(ic+1);
225 break;
226 }
227
228 sep= fabs( xi1- data.x2[ic]) +
229 fabs(yi1- data.y2[ic])+ fabs(zi1- data.z2[ic]);
230 if( sep <= slen)
231 {
232 data.icon1[i]= (ic+1);
233 break;
234 }
235
236 } /* for( j = 1; j < data.n; j++) */
237
238 } /* if( ! jump ) */
239
240 /* determine connection data for end 2 of segment. */
241 if( (ignd > 0) || jump )
242 {
243 if( zi2 <= -slen)
244 {
245 fprintf( stderr,
246 "xnec2c: conect(): geometry data error\n"
247 "segment %d extends below ground\n", iz );
248 stop( _("conect(): Geometry data error\n"\
249 "Segment extends below ground"), ERR_OK );
250 return( FALSE );
251 }
252
253 if( zi2 <= slen)
254 {
255 if( data.icon1[i] == iz )
256 {
257 fprintf( stderr,
258 "xnec2c: conect(): geometry data error\n"
259 "segment %d lies in ground plane\n", iz );
260 stop( _("conect(): Geometry data error\n"\
261 "Segment lies in ground plane"), ERR_OK );
262 return( FALSE );
263 }
264
265 data.icon2[i]= iz;
266 data.z2[i]=0.0;
267 continue;
268
269 } /* if( zi2 <= slen) */
270
271 } /* if( ignd > 0) */
272
273 ic= i;
274 for( j = 1; j < data.n; j++ )
275 {
276 ic++;
277 if( ic >= data.n)
278 ic=0;
279
280 sep= fabs(xi2- data.x1[ic]) +
281 fabs(yi2- data.y1[ic])+ fabs(zi2- data.z1[ic]);
282 if( sep <= slen)
283 {
284 data.icon2[i]= (ic+1);
285 break;
286 }
287
288 sep= fabs(xi2- data.x2[ic]) +
289 fabs(yi2- data.y2[ic])+ fabs(zi2- data.z2[ic]);
290 if( sep <= slen)
291 {
292 data.icon2[i]= -(ic+1);
293 break;
294 }
295
296 } /* for( j = 1; j < data.n; j++ ) */
297
298 } /* for( i = 0; i < data.n; i++ ) */
299
300 /* find wire-surface connections for new patches */
301 if( data.m != 0)
302 {
303 ix = -1;
304 i = 0;
305 while( ++i <= data.m )
306 {
307 ix++;
308 xs= data.px[ix];
309 ys= data.py[ix];
310 zs= data.pz[ix];
311
312 for( iseg = 0; iseg < data.n; iseg++ )
313 {
314 xi1= data.x1[iseg];
315 yi1= data.y1[iseg];
316 zi1= data.z1[iseg];
317 xi2= data.x2[iseg];
318 yi2= data.y2[iseg];
319 zi2= data.z2[iseg];
320
321 /* for first end of segment */
322 slen=( fabs(xi2- xi1) +
323 fabs(yi2- yi1)+ fabs(zi2- zi1))* SMIN;
324 sep= fabs(xi1- xs)+ fabs(yi1- ys)+ fabs(zi1- zs);
325
326 /* connection - divide patch into 4
327 * patches at present array loc. */
328 if( sep <= slen)
329 {
330 data.icon1[iseg]=PCHCON+ i;
331 ic=0;
332 subph( i, ic );
333 break;
334 }
335
336 sep= fabs(xi2- xs)+ fabs(yi2- ys)+ fabs(zi2- zs);
337 if( sep <= slen)
338 {
339 data.icon2[iseg]=PCHCON+ i;
340 ic=0;
341 subph( i, ic );
342 break;
343 }
344
345 } /* for( iseg = 0; iseg < data.n; iseg++ ) */
346
347 } /* while( ++i <= data.m ) */
348
349 } /* if( data.m != 0) */
350
351 } /* if( data.n != 0) */
352
353 iseg=( data.n+ data.m)/( data.np+ data.mp);
354 if( iseg != 1)
355 {
356 /*** may be error condition?? ***/
357 if( data.ipsym == 0 )
358 {
359 fprintf( stderr, "xnec2c: conect(): ipsym = 0\n" );
360 stop( _("conect(): Error: ipsym = 0"), ERR_OK );
361 return( FALSE );
362 }
363 } /* if( iseg != 1) */
364
365 /* No wire segments */
366 if( data.n <= 0) return( TRUE );
367
368 /* Allocate to connection buffers */
369 mreq = (size_t)segj.maxcon * sizeof(int);
370 mem_realloc( (void **)&segj.jco, mreq, "in geometry.c" );
371
372 /* Adjust connected segment ends to exactly coincide.
373 * Also find old seg. connecting to new segment */
374 iseg = 0;
375 for( j = 0; j < data.n; j++ )
376 {
377 jx = j+1;
378 iend=-1;
379 jend=-1;
380 ix= data.icon1[j];
381 ic=1;
382 segj.jco[0]= -jx;
383 xa= data.x1[j];
384 ya= data.y1[j];
385 za= data.z1[j];
386
387 while( TRUE )
388 {
389 if( (ix != 0) && (ix != (j+1)) && (ix <= PCHCON) )
390 {
391 do
392 {
393 if( ix < 0 )
394 ix= -ix;
395 else
396 jend= -jend;
397
398 jump = FALSE;
399
400 if( ix == jx )
401 break;
402
403 if( ix < jx )
404 {
405 jump = TRUE;
406 break;
407 }
408
409 /* Record max. no. of connections */
410 ic++;
411 if( ic >= segj.maxcon )
412 {
413 segj.maxcon = ic+1;
414 mreq = (size_t)segj.maxcon * sizeof(int);
415 mem_realloc( (void **)&segj.jco, mreq, "in geometry.c" );
416 }
417 segj.jco[ic-1]= ix* jend;
418
419 ixx = ix-1;
420 if( jend != 1)
421 {
422 xa= xa+ data.x1[ixx];
423 ya= ya+ data.y1[ixx];
424 za= za+ data.z1[ixx];
425 ix= data.icon1[ixx];
426 continue;
427 }
428
429 xa= xa+ data.x2[ixx];
430 ya= ya+ data.y2[ixx];
431 za= za+ data.z2[ixx];
432 ix= data.icon2[ixx];
433
434 } /* do */
435 while( ix != 0 );
436
437 if( jump && (iend == 1) ) break;
438 else if( jump )
439 {
440 iend=1;
441 jend=1;
442 ix= data.icon2[j];
443 ic=1;
444 segj.jco[0]= jx;
445 xa= data.x2[j];
446 ya= data.y2[j];
447 za= data.z2[j];
448 continue;
449 }
450
451 sep= (double)ic;
452 xa= xa/ sep;
453 ya= ya/ sep;
454 za= za/ sep;
455
456 for( i = 0; i < ic; i++ )
457 {
458 ix= segj.jco[i];
459 if( ix <= 0)
460 {
461 ix= -ix;
462 ixx = ix-1;
463 data.x1[ixx]= xa;
464 data.y1[ixx]= ya;
465 data.z1[ixx]= za;
466 continue;
467 }
468
469 ixx = ix-1;
470 data.x2[ixx]= xa;
471 data.y2[ixx]= ya;
472 data.z2[ixx]= za;
473
474 } /* for( i = 0; i < ic; i++ ) */
475
476 if( ic >= 3)
477 iseg++;
478
479 } /*if( (ix != 0) && (ix != j) && (ix <= PCHCON) ) */
480
481 if( iend == 1) break;
482
483 iend=1;
484 jend=1;
485 ix= data.icon2[j];
486 ic=1;
487 segj.jco[0]= jx;
488 xa= data.x2[j];
489 ya= data.y2[j];
490 za= data.z2[j];
491
492 } /* while( TRUE ) */
493
494 } /* for( j = 0; j < data.n; j++ ) */
495
496 mreq = (size_t)segj.maxcon * sizeof(double);
497 mem_realloc( (void **)&segj.ax, mreq, "in geometry.c" );
498 mem_realloc( (void **)&segj.bx, mreq, "in geometry.c" );
499 mem_realloc( (void **)&segj.cx, mreq, "in geometry.c" );
500
501 return( TRUE );
502 }
503
504 /*-----------------------------------------------------------------------*/
505
506 /* subroutine helix generates segment geometry */
507 /* data for a helix of ns segments */
508 void
helix(double s,double hl,double a1,double b1,double a2,double b2,double rad,int ns,int itg)509 helix(
510 double s, double hl,
511 double a1, double b1,
512 double a2, double b2,
513 double rad, int ns, int itg )
514 {
515 int ist, i;
516 size_t mreq;
517 double zinc, copy;
518
519 ist = data.n;
520 data.n += ns;
521 data.np = data.n;
522 data.mp = data.m;
523 data.ipsym = 0;
524
525 if( ns < 1) return;
526
527 zinc= fabs( hl/ ns);
528
529 /* Reallocate tags buffer */
530 mreq = (size_t)data.n * sizeof(int);
531 mem_realloc( (void **)&data.itag, mreq, "in geometry.c" );
532
533 /* Reallocate wire buffers */
534 mreq = (size_t)data.n * sizeof(double);
535 mem_realloc( (void **)&data.x1, mreq, "in geometry.c" );
536 mem_realloc( (void **)&data.y1, mreq, "in geometry.c" );
537 mem_realloc( (void **)&data.z1, mreq, "in geometry.c" );
538 mem_realloc( (void **)&data.x2, mreq, "in geometry.c" );
539 mem_realloc( (void **)&data.y2, mreq, "in geometry.c" );
540 mem_realloc( (void **)&data.z2, mreq, "in geometry.c" );
541 mem_realloc( (void **)&data.bi, mreq, "in geometry.c" );
542
543 data.z1[ist] = 0.0;
544 for( i = ist; i < data.n; i++ )
545 {
546 data.bi[i]= rad;
547 data.itag[i]= itg;
548
549 if( i != ist )
550 data.z1[i]= data.z1[i-1]+ zinc;
551
552 data.z2[i]= data.z1[i]+ zinc;
553
554 if( a2 == a1)
555 {
556 if( b1 == 0.0)
557 b1= a1;
558
559 data.x1[i]= a1* cos(2.0* PI* data.z1[i]/ s);
560 data.y1[i]= b1* sin(2.0* PI* data.z1[i]/ s);
561 data.x2[i]= a1* cos(2.0* PI* data.z2[i]/ s);
562 data.y2[i]= b1* sin(2.0* PI* data.z2[i]/ s);
563 }
564 else
565 {
566 if( b2 == 0.0)
567 b2= a2;
568
569 data.x1[i]=( a1+( a2- a1)* data.z1[i] /
570 fabs( hl))* cos(2.0* PI* data.z1[i]/ s);
571 data.y1[i]=( b1+( b2- b1)* data.z1[i] /
572 fabs( hl))* sin(2.0* PI* data.z1[i]/ s);
573 data.x2[i]=( a1+( a2- a1)* data.z2[i] /
574 fabs( hl))* cos(2.0* PI* data.z2[i]/ s);
575 data.y2[i]=( b1+( b2- b1)* data.z2[i] /
576 fabs( hl))* sin(2.0* PI* data.z2[i]/ s);
577
578 } /* if( a2 == a1) */
579
580 if( hl > 0.0)
581 continue;
582
583 copy= data.x1[i];
584 data.x1[i]= data.y1[i];
585 data.y1[i]= copy;
586 copy= data.x2[i];
587 data.x2[i]= data.y2[i];
588 data.y2[i]= copy;
589
590 } /* for( i = ist; i < data.n; i++ ) */
591
592 return;
593 }
594
595 /*-----------------------------------------------------------------------*/
596
597 /* isegno returns the segment number of the mth segment having the */
598 /* tag number itagi. if itagi=0 segment number m is returned. */
599 int
isegno(int itagi,int mx)600 isegno( int itagi, int mx)
601 {
602 int icnt, iseg;
603
604 if( mx <= 0)
605 {
606 fprintf( stderr,
607 "xnec2c: isegno(): check data, parameter specifying segment\n"
608 "position in a group of equal tags must not be zero\n" );
609 stop( _("isegno(): Parameter specifying segment\n"\
610 "position in a group of equal tags is zero"), ERR_OK );
611 return( -1 );
612 }
613
614 icnt=0;
615 if( itagi == 0)
616 {
617 iseg = mx;
618 return( iseg );
619 }
620
621 if( data.n > 0)
622 {
623 int i;
624
625 for( i = 0; i < data.n; i++ )
626 {
627 if( data.itag[i] != itagi )
628 continue;
629
630 icnt++;
631 if( icnt == mx)
632 {
633 iseg= i+1;
634 return( iseg );
635 }
636
637 } /* for( i = 0; i < data.n; i++ ) */
638
639 } /* if( data.n > 0) */
640
641 fprintf( stderr,
642 "xnec2c: isegno(): no segment has an itag of %d\n", itagi );
643 stop( _("isegno(): Segment tag number error"), ERR_OK );
644
645 return( -1 );
646 }
647
648 /*-----------------------------------------------------------------------*/
649
650 /* subroutine move moves the structure with respect to its */
651 /* coordinate system or reproduces structure in new positions. */
652 /* structure is rotated about x,y,z axes by rox,roy,roz */
653 /* respectively, then shifted by xs,ys,zs */
654 gboolean
move(double rox,double roy,double roz,double xs,double ys,double zs,int its,int nrpt,int itgi)655 move( double rox, double roy,
656 double roz, double xs, double ys,
657 double zs, int its, int nrpt, int itgi )
658 {
659 int nrp, ix, i1, k, i;
660 size_t mreq;
661 double sps, cps, sth, cth, sph, cph, xx, xy;
662 double xz, yx, yy, yz, zx, zy, zz, xi, yi, zi;
663
664 if( fabs( rox)+ fabs( roy) > 1.0e-10)
665 data.ipsym= data.ipsym*3;
666
667 sps= sin( rox);
668 cps= cos( rox);
669 sth= sin( roy);
670 cth= cos( roy);
671 sph= sin( roz);
672 cph= cos( roz);
673 xx= cph* cth;
674 xy= cph* sth* sps- sph* cps;
675 xz= cph* sth* cps+ sph* sps;
676 yx= sph* cth;
677 yy= sph* sth* sps+ cph* cps;
678 yz= sph* sth* cps- cph* sps;
679 zx= -sth;
680 zy= cth* sps;
681 zz= cth* cps;
682
683 if( nrpt == 0) nrp=1;
684 else nrp= nrpt;
685
686 ix=1;
687 if( data.n > 0)
688 {
689 int ir;
690
691 i1= isegno( its, 1);
692 if( i1 < 0 ) return( FALSE ); /* my addition, error */
693 if( i1 < 1) i1= 1;
694
695 ix= i1;
696 if( nrpt == 0) k= i1-1;
697 else
698 {
699 k= data.n;
700 /* Reallocate tags buffer */
701 mreq = (size_t)(data.n + data.m + (data.n + 1 - i1) * nrpt) * sizeof(int);
702 mem_realloc( (void **)&data.itag, mreq, "in geometry.c" );
703
704 /* Reallocate wire buffers */
705 mreq = (size_t)(data.n + (data.n + 1 - i1) * nrpt) * sizeof(double);
706 mem_realloc( (void **)&data.x1, mreq, "in geometry.c" );
707 mem_realloc( (void **)&data.y1, mreq, "in geometry.c" );
708 mem_realloc( (void **)&data.z1, mreq, "in geometry.c" );
709 mem_realloc( (void **)&data.x2, mreq, "in geometry.c" );
710 mem_realloc( (void **)&data.y2, mreq, "in geometry.c" );
711 mem_realloc( (void **)&data.z2, mreq, "in geometry.c" );
712 mem_realloc( (void **)&data.bi, mreq, "in geometry.c" );
713 }
714
715 for( ir = 0; ir < nrp; ir++ )
716 {
717 for( i = i1-1; i < data.n; i++ )
718 {
719 xi= data.x1[i];
720 yi= data.y1[i];
721 zi= data.z1[i];
722 data.x1[k]= xi* xx+ yi* xy+ zi* xz+ xs;
723 data.y1[k]= xi* yx+ yi* yy+ zi* yz+ ys;
724 data.z1[k]= xi* zx+ yi* zy+ zi* zz+ zs;
725 xi= data.x2[i];
726 yi= data.y2[i];
727 zi= data.z2[i];
728 data.x2[k]= xi* xx+ yi* xy+ zi* xz+ xs;
729 data.y2[k]= xi* yx+ yi* yy+ zi* yz+ ys;
730 data.z2[k]= xi* zx+ yi* zy+ zi* zz+ zs;
731 data.bi[k]= data.bi[i];
732 data.itag[k]= data.itag[i];
733 if( data.itag[i] != 0)
734 data.itag[k]= data.itag[i]+ itgi;
735
736 k++;
737
738 } /* for( i = i1; i < data.n; i++ ) */
739
740 i1= data.n+1;
741 data.n= k;
742
743 } /* for( ir = 0; ir < nrp; ir++ ) */
744
745 } /* if( data.n >= n2) */
746
747 if( data.m > 0)
748 {
749 int ii;
750 i1 = 0;
751 if( nrpt == 0) k= 0;
752 else k = data.m;
753
754 /* Reallocate patch buffers */
755 mreq = (size_t)(data.m * (nrpt + 1)) * sizeof(double);
756 mem_realloc( (void **)&data.px, mreq, "in geometry.c" );
757 mem_realloc( (void **)&data.py, mreq, "in geometry.c" );
758 mem_realloc( (void **)&data.pz, mreq, "in geometry.c" );
759 mem_realloc( (void **)&data.t1x, mreq, "in geometry.c" );
760 mem_realloc( (void **)&data.t1y, mreq, "in geometry.c" );
761 mem_realloc( (void **)&data.t1z, mreq, "in geometry.c" );
762 mem_realloc( (void **)&data.t2x, mreq, "in geometry.c" );
763 mem_realloc( (void **)&data.t2y, mreq, "in geometry.c" );
764 mem_realloc( (void **)&data.t2z, mreq, "in geometry.c" );
765 mem_realloc( (void **)&data.pbi, mreq, "in geometry.c" );
766 mem_realloc( (void **)&data.psalp, mreq, "in geometry.c" );
767
768 for( ii = 0; ii < nrp; ii++ )
769 {
770 for( i = i1; i < data.m; i++ )
771 {
772 xi= data.px[i];
773 yi= data.py[i];
774 zi= data.pz[i];
775 data.px[k]= xi* xx+ yi* xy+ zi* xz+ xs;
776 data.py[k]= xi* yx+ yi* yy+ zi* yz+ ys;
777 data.pz[k]= xi* zx+ yi* zy+ zi* zz+ zs;
778 xi= data.t1x[i];
779 yi= data.t1y[i];
780 zi= data.t1z[i];
781 data.t1x[k]= xi* xx+ yi* xy+ zi* xz;
782 data.t1y[k]= xi* yx+ yi* yy+ zi* yz;
783 data.t1z[k]= xi* zx+ yi* zy+ zi* zz;
784 xi= data.t2x[i];
785 yi= data.t2y[i];
786 zi= data.t2z[i];
787 data.t2x[k]= xi* xx+ yi* xy+ zi* xz;
788 data.t2y[k]= xi* yx+ yi* yy+ zi* yz;
789 data.t2z[k]= xi* zx+ yi* zy+ zi* zz;
790 data.psalp[k]= data.psalp[i];
791 data.pbi[k]= data.pbi[i];
792 k++;
793
794 } /* for( i = i1; i < data.m; i++ ) */
795
796 i1= data.m;
797 data.m = k;
798
799 } /* for( ii = 0; ii < nrp; ii++ ) */
800
801 } /* if( data.m >= m2) */
802
803 if( (nrpt == 0) && (ix == 1) )
804 return( TRUE );
805
806 data.np= data.n;
807 data.mp= data.m;
808 data.ipsym=0;
809
810 return( TRUE );
811 }
812
813 /*-----------------------------------------------------------------------*/
814
815 /* patch generates and modifies patch geometry data */
816 gboolean
patch(int nx,int ny,double ax1,double ay1,double az1,double ax2,double ay2,double az2,double ax3,double ay3,double az3,double ax4,double ay4,double az4)817 patch(
818 int nx, int ny,
819 double ax1, double ay1, double az1,
820 double ax2, double ay2, double az2,
821 double ax3, double ay3, double az3,
822 double ax4, double ay4, double az4 )
823 {
824 int mi, ntp;
825 size_t mreq;
826 double s1x=0.0, s1y=0.0, s1z=0.0;
827 double s2x=0.0, s2y=0.0, s2z=0.0, xst=0.0;
828 double znv, xnv, ynv, xa, xn2, yn2;
829 double zn2;
830
831 /* new patches. for nx=0, ny=1,2,3,4 patch is (respectively) */;
832 /* arbitrary, rectagular, triangular, or quadrilateral. */
833 /* for nx and ny > 0 a rectangular surface is produced with */
834 /* nx by ny rectangular patches. */
835
836 data.m++;
837 mi= data.m-1;
838
839 /* Reallocate patch buffers */
840 mreq = (size_t)data.m * sizeof(double);
841 mem_realloc( (void **)&data.px, mreq, "in geometry.c" );
842 mem_realloc( (void **)&data.py, mreq, "in geometry.c" );
843 mem_realloc( (void **)&data.pz, mreq, "in geometry.c" );
844 mem_realloc( (void **)&data.t1x, mreq, "in geometry.c" );
845 mem_realloc( (void **)&data.t1y, mreq, "in geometry.c" );
846 mem_realloc( (void **)&data.t1z, mreq, "in geometry.c" );
847 mem_realloc( (void **)&data.t2x, mreq, "in geometry.c" );
848 mem_realloc( (void **)&data.t2y, mreq, "in geometry.c");
849 mem_realloc( (void **)&data.t2z, mreq, "in geometry.c" );
850 mem_realloc( (void **)&data.pbi, mreq, "in geometry.c" );
851 mem_realloc( (void **)&data.psalp, mreq, "in geometry.c" );
852
853 if( nx > 0) ntp=2;
854 else ntp= ny;
855
856 if( ntp <= 1)
857 {
858 data.px[mi]= ax1;
859 data.py[mi]= ay1;
860 data.pz[mi]= az1;
861 data.pbi[mi]= az2;
862 znv= cos( ax2);
863 xnv= znv* cos( ay2);
864 ynv= znv* sin( ay2);
865 znv= sin( ax2);
866 xa= sqrt( xnv* xnv+ ynv* ynv);
867
868 if( xa >= 1.0e-6)
869 {
870 data.t1x[mi]= -ynv/ xa;
871 data.t1y[mi]= xnv/ xa;
872 data.t1z[mi]=0.0;
873 }
874 else
875 {
876 data.t1x[mi]=1.0;
877 data.t1y[mi]=0.0;
878 data.t1z[mi]=0.0;
879 }
880
881 } /* if( ntp <= 1) */
882 else
883 {
884 s1x= ax2- ax1;
885 s1y= ay2- ay1;
886 s1z= az2- az1;
887 s2x= ax3- ax2;
888 s2y= ay3- ay2;
889 s2z= az3- az2;
890
891 if( nx != 0)
892 {
893 s1x= s1x/ nx;
894 s1y= s1y/ nx;
895 s1z= s1z/ nx;
896 s2x= s2x/ ny;
897 s2y= s2y/ ny;
898 s2z= s2z/ ny;
899 }
900
901 xnv= s1y* s2z- s1z* s2y;
902 ynv= s1z* s2x- s1x* s2z;
903 znv= s1x* s2y- s1y* s2x;
904 xa= sqrt( xnv* xnv+ ynv* ynv+ znv* znv);
905 xnv= xnv/ xa;
906 ynv= ynv/ xa;
907 znv= znv/ xa;
908 xst= sqrt( s1x* s1x+ s1y* s1y+ s1z* s1z);
909 data.t1x[mi]= s1x/ xst;
910 data.t1y[mi]= s1y/ xst;
911 data.t1z[mi]= s1z/ xst;
912
913 if( ntp <= 2)
914 {
915 data.px[mi]= ax1+.5*( s1x+ s2x);
916 data.py[mi]= ay1+.5*( s1y+ s2y);
917 data.pz[mi]= az1+.5*( s1z+ s2z);
918 data.pbi[mi]= xa;
919 }
920 else
921 {
922 if( ntp != 4)
923 {
924 data.px[mi]=( ax1+ ax2+ ax3)/3.0;
925 data.py[mi]=( ay1+ ay2+ ay3)/3.0;
926 data.pz[mi]=( az1+ az2+ az3)/3.0;
927 data.pbi[mi]=.5* xa;
928 }
929 else
930 {
931 s1x= ax3- ax1;
932 s1y= ay3- ay1;
933 s1z= az3- az1;
934 s2x= ax4- ax1;
935 s2y= ay4- ay1;
936 s2z= az4- az1;
937 xn2= s1y* s2z- s1z* s2y;
938 yn2= s1z* s2x- s1x* s2z;
939 zn2= s1x* s2y- s1y* s2x;
940 xst= sqrt( xn2* xn2+ yn2* yn2+ zn2* zn2);
941 double salpn=1.0/(3.0*( xa+ xst));
942 data.px[mi]=( xa*( ax1+ ax2+ ax3) +
943 xst*( ax1+ ax3+ ax4))* salpn;
944 data.py[mi]=( xa*( ay1+ ay2+ ay3) +
945 xst*( ay1+ ay3+ ay4))* salpn;
946 data.pz[mi]=( xa*( az1+ az2+ az3) +
947 xst*( az1+ az3+ az4))* salpn;
948 data.pbi[mi]=.5*( xa+ xst);
949 s1x=( xnv* xn2+ ynv* yn2+ znv* zn2)/ xst;
950
951 if( s1x <= 0.9998)
952 {
953 fprintf( stderr,
954 "xnec2c: patch(): corners of quadrilateral\n"
955 "patch do not lie in a plane\n" );
956 stop( _("patch(): Corners of quadrilateral\n"\
957 "patch do not lie in a plane"), ERR_OK );
958 return( FALSE );
959 }
960
961 } /* if( ntp != 4) */
962
963 } /* if( ntp <= 2) */
964
965 } /* if( ntp <= 1) */
966
967 data.t2x[mi]= ynv* data.t1z[mi]- znv* data.t1y[mi];
968 data.t2y[mi]= znv* data.t1x[mi]- xnv* data.t1z[mi];
969 data.t2z[mi]= xnv* data.t1y[mi]- ynv* data.t1x[mi];
970 data.psalp[mi]=1.0;
971
972 if( nx != 0)
973 {
974 int iy, ix;
975 double xs, ys, zs, xt, yt, zt;
976
977 data.m += nx*ny-1;
978 /* Reallocate patch buffers */
979 mreq = (size_t)data.m * sizeof(double);
980 mem_realloc( (void **)&data.px, mreq, "in geometry.c" );
981 mem_realloc( (void **)&data.py, mreq, "in geometry.c" );
982 mem_realloc( (void **)&data.pz, mreq, "in geometry.c" );
983 mem_realloc( (void **)&data.t1x, mreq, "in geometry.c" );
984 mem_realloc( (void **)&data.t1y, mreq, "in geometry.c" );
985 mem_realloc( (void **)&data.t1z, mreq, "in geometry.c" );
986 mem_realloc( (void **)&data.t2x, mreq, "in geometry.c" );
987 mem_realloc( (void **)&data.t2y, mreq, "in geometry.c" );
988 mem_realloc( (void **)&data.t2z, mreq, "in geometry.c" );
989 mem_realloc( (void **)&data.pbi, mreq, "in geometry.c" );
990 mem_realloc( (void **)&data.psalp, mreq, "in geometry.c" );
991
992 xn2= data.px[mi]- s1x- s2x;
993 yn2= data.py[mi]- s1y- s2y;
994 zn2= data.pz[mi]- s1z- s2z;
995 xs= data.t1x[mi];
996 ys= data.t1y[mi];
997 zs= data.t1z[mi];
998 xt= data.t2x[mi];
999 yt= data.t2y[mi];
1000 zt= data.t2z[mi];
1001
1002 for( iy = 0; iy < ny; iy++ )
1003 {
1004 xn2 += s2x;
1005 yn2 += s2y;
1006 zn2 += s2z;
1007
1008 for( ix = 1; ix <= nx; ix++ )
1009 {
1010 xst= (double)ix;
1011 data.px[mi]= xn2+ xst* s1x;
1012 data.py[mi]= yn2+ xst* s1y;
1013 data.pz[mi]= zn2+ xst* s1z;
1014 data.pbi[mi]= xa;
1015 data.psalp[mi]=1.0;
1016 data.t1x[mi]= xs;
1017 data.t1y[mi]= ys;
1018 data.t1z[mi]= zs;
1019 data.t2x[mi]= xt;
1020 data.t2y[mi]= yt;
1021 data.t2z[mi]= zt;
1022 mi++;
1023 } /* for( ix = 0; ix < nx; ix++ ) */
1024
1025 } /* for( iy = 0; iy < ny; iy++ ) */
1026
1027 } /* if( nx != 0) */
1028
1029 data.ipsym=0;
1030 data.np= data.n;
1031 data.mp= data.m;
1032
1033 return( TRUE );
1034 }
1035
1036 /*-----------------------------------------------------------------------*/
1037
1038 /*** this function was an 'entry point' (part of) 'patch()' ***/
1039 void
subph(int nx,int ny)1040 subph( int nx, int ny )
1041 {
1042 int mia, ix, iy, mi;
1043 size_t mreq;
1044 double xs, ys, zs, xa, xst, s1x, s1y;
1045 double s1z, s2x, s2y, s2z, saln, xt, yt;
1046
1047 /* Reallocate patch buffers */
1048 if( ny == 0 ) data.m += 3;
1049 else data.m += 4;
1050
1051 mreq = (size_t)data.m * sizeof(double);
1052 mem_realloc( (void **)&data.px, mreq, "in geometry.c" );
1053 mem_realloc( (void **)&data.py, mreq, "in geometry.c" );
1054 mem_realloc( (void **)&data.pz, mreq, "in geometry.c" );
1055 mem_realloc( (void **)&data.t1x, mreq, "in geometry.c" );
1056 mem_realloc( (void **)&data.t1y, mreq, "in geometry.c" );
1057 mem_realloc( (void **)&data.t1z, mreq, "in geometry.c" );
1058 mem_realloc( (void **)&data.t2x, mreq, "in geometry.c" );
1059 mem_realloc( (void **)&data.t2y, mreq, "in geometry.c" );
1060 mem_realloc( (void **)&data.t2z, mreq, "in geometry.c" );
1061 mem_realloc( (void **)&data.pbi, mreq, "in geometry.c" );
1062 mem_realloc( (void **)&data.psalp, mreq, "in geometry.c" );
1063 if( !CHILD )
1064 {
1065 mreq = (size_t)(data.n + 2 * data.m) * sizeof(GdkSegment);
1066 mem_realloc( (void **)&structure_segs, mreq, "in geometry.c" );
1067 }
1068
1069 /* Shift patches to make room for new ones */
1070 if( (ny == 0) && (nx != data.m) )
1071 {
1072 for( iy = data.m-1; iy > nx+2; iy-- )
1073 {
1074 ix = iy-3;
1075 data.px[iy]= data.px[ix];
1076 data.py[iy]= data.py[ix];
1077 data.pz[iy]= data.pz[ix];
1078 data.pbi[iy]= data.pbi[ix];
1079 data.psalp[iy]= data.psalp[ix];
1080 data.t1x[iy]= data.t1x[ix];
1081 data.t1y[iy]= data.t1y[ix];
1082 data.t1z[iy]= data.t1z[ix];
1083 data.t2x[iy]= data.t2x[ix];
1084 data.t2y[iy]= data.t2y[ix];
1085 data.t2z[iy]= data.t2z[ix];
1086 }
1087
1088 } /* if( (ny == 0) || (nx != m) ) */
1089
1090 /* divide patch for connection */
1091 mi= nx-1;
1092 xs= data.px[mi];
1093 ys= data.py[mi];
1094 zs= data.pz[mi];
1095 xa= data.pbi[mi]/4.0;
1096 xst= sqrt( xa)/2.0;
1097 s1x= data.t1x[mi];
1098 s1y= data.t1y[mi];
1099 s1z= data.t1z[mi];
1100 s2x= data.t2x[mi];
1101 s2y= data.t2y[mi];
1102 s2z= data.t2z[mi];
1103 saln= data.psalp[mi];
1104 xt= xst;
1105 yt= xst;
1106
1107 if( ny == 0)
1108 mia= mi;
1109 else
1110 {
1111 data.mp++;
1112 mia= data.m-1;
1113 }
1114
1115 for( ix = 1; ix <= 4; ix++ )
1116 {
1117 data.px[mia]= xs+ xt* s1x+ yt* s2x;
1118 data.py[mia]= ys+ xt* s1y+ yt* s2y;
1119 data.pz[mia]= zs+ xt* s1z+ yt* s2z;
1120 data.pbi[mia]= xa;
1121 data.t1x[mia]= s1x;
1122 data.t1y[mia]= s1y;
1123 data.t1z[mia]= s1z;
1124 data.t2x[mia]= s2x;
1125 data.t2y[mia]= s2y;
1126 data.t2z[mia]= s2z;
1127 data.psalp[mia]= saln;
1128
1129 if( ix == 2)
1130 yt= -yt;
1131
1132 if( (ix == 1) || (ix == 3) )
1133 xt= -xt;
1134
1135 mia++;
1136 }
1137
1138 if( nx <= data.mp)
1139 data.mp += 3;
1140
1141 if( ny > 0 )
1142 data.pz[mi]=10000.0;
1143
1144 /* Process new patches created */
1145 if( ! CHILD )
1146 New_Patch_Data();
1147
1148 return;
1149 }
1150
1151 /*-----------------------------------------------------------------------*/
1152
1153 /* reflc reflects partial structure along x,y, or z axes */
1154 /* or rotates structure to complete a symmetric structure. */
1155 gboolean
reflc(int ix,int iy,int iz,int iti,int nop)1156 reflc( int ix, int iy, int iz, int iti, int nop )
1157 {
1158 int i, nx, itagi, k;
1159 size_t mreq;
1160 double e1, e2, fnop, sam, cs, ss, xk, yk;
1161
1162 if( nop == 0) return( TRUE );
1163
1164 data.np= data.n;
1165 data.mp= data.m;
1166 data.ipsym=0;
1167
1168 if( ix >= 0)
1169 {
1170 data.ipsym=1;
1171
1172 /* reflect along z axis */
1173 if( iz != 0)
1174 {
1175 data.ipsym=2;
1176
1177 if( data.n > 0 )
1178 {
1179 /* Reallocate tags buffer */
1180 mreq = (size_t)(2 * data.n + data.m) * sizeof(int);
1181 mem_realloc( (void **)&data.itag, mreq, "in geometry.c" );
1182
1183 /* Reallocate wire buffers */
1184 mreq = (size_t)(2 * data.n) * sizeof(double);
1185 mem_realloc( (void **)&data.x1, mreq, "in geometry.c" );
1186 mem_realloc( (void **)&data.y1, mreq, "in geometry.c" );
1187 mem_realloc( (void **)&data.z1, mreq, "in geometry.c" );
1188 mem_realloc( (void **)&data.x2, mreq, "in geometry.c" );
1189 mem_realloc( (void **)&data.y2, mreq, "in geometry.c" );
1190 mem_realloc( (void **)&data.z2, mreq, "in geometry.c" );
1191 mem_realloc( (void **)&data.bi, mreq, "in geometry.c" );
1192
1193 for( i = 0; i < data.n; i++ )
1194 {
1195 nx= i+ data.n;
1196 e1= data.z1[i];
1197 e2= data.z2[i];
1198
1199 if( (fabs(e1)+fabs(e2) <= 1.0e-5) || (e1*e2 < -1.0e-6) )
1200 {
1201 fprintf( stderr,
1202 "xnec2c: reflc(): geometry data error\n"
1203 "segment %d lies in plane of symmetry\n", i+1 );
1204 stop( _("reflc(): Geometry data error\n"\
1205 "Segment lies in plane of symmetry"), ERR_OK );
1206 return( FALSE );
1207 }
1208
1209 data.x1[nx]= data.x1[i];
1210 data.y1[nx]= data.y1[i];
1211 data.z1[nx]= -e1;
1212 data.x2[nx]= data.x2[i];
1213 data.y2[nx]= data.y2[i];
1214 data.z2[nx]= -e2;
1215 itagi= data.itag[i];
1216
1217 if( itagi == 0)
1218 data.itag[nx]=0;
1219 if( itagi != 0)
1220 data.itag[nx]= itagi+ iti;
1221
1222 data.bi[nx]= data.bi[i];
1223
1224 } /* for( i = 0; i < data.n; i++ ) */
1225
1226 data.n= data.n*2;
1227 iti= iti*2;
1228
1229 } /* if( data.n > 0) */
1230
1231 if( data.m > 0 )
1232 {
1233 /* Reallocate patch buffers */
1234 mreq = (size_t)(2 * data.m) * sizeof(double);
1235 mem_realloc( (void **)&data.px, mreq, "in geometry.c" );
1236 mem_realloc( (void **)&data.py, mreq, "in geometry.c" );
1237 mem_realloc( (void **)&data.pz, mreq, "in geometry.c" );
1238 mem_realloc( (void **)&data.t1x, mreq, "in geometry.c" );
1239 mem_realloc( (void **)&data.t1y, mreq, "in geometry.c" );
1240 mem_realloc( (void **)&data.t1z, mreq, "in geometry.c" );
1241 mem_realloc( (void **)&data.t2x, mreq, "in geometry.c" );
1242 mem_realloc( (void **)&data.t2y, mreq, "in geometry.c" );
1243 mem_realloc( (void **)&data.t2z, mreq, "in geometry.c" );
1244 mem_realloc( (void **)&data.pbi, mreq, "in geometry.c" );
1245 mem_realloc( (void **)&data.psalp, mreq, "in geometry.c" );
1246
1247 for( i = 0; i < data.m; i++ )
1248 {
1249 nx = i+data.m;
1250 if( fabs(data.pz[i]) <= 1.0e-10)
1251 {
1252 fprintf( stderr,
1253 "xnec2c: reflc(): geometry data error\n"
1254 "patch %d lies in plane of symmetry\n", i+1 );
1255 stop( _("reflc(): Geometry data error\n"\
1256 "Patch lies in plane of symmetry"), ERR_OK );
1257 return( FALSE );
1258 }
1259
1260 data.px[nx]= data.px[i];
1261 data.py[nx]= data.py[i];
1262 data.pz[nx]= -data.pz[i];
1263 data.t1x[nx]= data.t1x[i];
1264 data.t1y[nx]= data.t1y[i];
1265 data.t1z[nx]= -data.t1z[i];
1266 data.t2x[nx]= data.t2x[i];
1267 data.t2y[nx]= data.t2y[i];
1268 data.t2z[nx]= -data.t2z[i];
1269 data.psalp[nx]= -data.psalp[i];
1270 data.pbi[nx]= data.pbi[i];
1271 }
1272
1273 data.m= data.m*2;
1274
1275 } /* if( data.m >= m2) */
1276
1277 } /* if( iz != 0) */
1278
1279 /* reflect along y axis */
1280 if( iy != 0)
1281 {
1282 if( data.n > 0)
1283 {
1284 /* Reallocate tags buffer */
1285 mreq = (size_t)(2 * data.n) * sizeof(int);
1286 mem_realloc( (void **)&data.itag, mreq, "in geometry.c" );
1287 /* Reallocate wire buffers */
1288 mreq = (size_t)(2 * data.n) * sizeof(double);
1289 mem_realloc( (void **)&data.x1, mreq, "in geometry.c" );
1290 mem_realloc( (void **)&data.y1, mreq, "in geometry.c" );
1291 mem_realloc( (void **)&data.z1, mreq, "in geometry.c" );
1292 mem_realloc( (void **)&data.x2, mreq, "in geometry.c" );
1293 mem_realloc( (void **)&data.y2, mreq, "in geometry.c" );
1294 mem_realloc( (void **)&data.z2, mreq, "in geometry.c" );
1295 mem_realloc( (void **)&data.bi, mreq, "in geometry.c" );
1296
1297 for( i = 0; i < data.n; i++ )
1298 {
1299 nx= i+ data.n;
1300 e1= data.y1[i];
1301 e2= data.y2[i];
1302
1303 if( (fabs(e1)+fabs(e2) <= 1.0e-5) || (e1*e2 < -1.0e-6) )
1304 {
1305 fprintf( stderr,
1306 "xnec2c: reflc(): geometry data error\n"
1307 "segment %d lies in plane of symmetry\n", i+1 );
1308 stop( _("reflc(): Geometry data error\n"\
1309 "Segment lies in plane of symmetry"), ERR_OK );
1310 return( FALSE );
1311 }
1312
1313 data.x1[nx]= data.x1[i];
1314 data.y1[nx]= -e1;
1315 data.z1[nx]= data.z1[i];
1316 data.x2[nx]= data.x2[i];
1317 data.y2[nx]= -e2;
1318 data.z2[nx]= data.z2[i];
1319 itagi= data.itag[i];
1320
1321 if( itagi == 0)
1322 data.itag[nx]=0;
1323 if( itagi != 0)
1324 data.itag[nx]= itagi+ iti;
1325
1326 data.bi[nx]= data.bi[i];
1327
1328 } /* for( i = n2-1; i < data.n; i++ ) */
1329
1330 data.n= data.n*2;
1331 iti= iti*2;
1332
1333 } /* if( data.n >= n2) */
1334
1335 if( data.m > 0 )
1336 {
1337 /* Reallocate patch buffers */
1338 mreq = (size_t)(2 * data.m) * sizeof(double);
1339 mem_realloc( (void **)&data.px, mreq, "in geometry.c" );
1340 mem_realloc( (void **)&data.py, mreq, "in geometry.c" );
1341 mem_realloc( (void **)&data.pz, mreq, "in geometry.c" );
1342 mem_realloc( (void **)&data.t1x, mreq, "in geometry.c" );
1343 mem_realloc( (void **)&data.t1y, mreq, "in geometry.c" );
1344 mem_realloc( (void **)&data.t1z, mreq, "in geometry.c" );
1345 mem_realloc( (void **)&data.t2x, mreq, "in geometry.c" );
1346 mem_realloc( (void **)&data.t2y, mreq, "in geometry.c" );
1347 mem_realloc( (void **)&data.t2z, mreq, "in geometry.c" );
1348 mem_realloc( (void **)&data.pbi, mreq, "in geometry.c" );
1349 mem_realloc( (void **)&data.psalp, mreq, "in geometry.c" );
1350
1351 for( i = 0; i < data.m; i++ )
1352 {
1353 nx= i+data.m;
1354 if( fabs( data.py[i]) <= 1.0e-10)
1355 {
1356 fprintf( stderr,
1357 "xnec2c: reflc(): geometry data error\n"
1358 "patch %d lies in plane of symmetry\n", i+1 );
1359 stop( _("reflc(): Geometry data error\n"\
1360 "Patch lies in plane of symmetry"), ERR_OK );
1361 return( FALSE );
1362 }
1363
1364 data.px[nx]= data.px[i];
1365 data.py[nx]= -data.py[i];
1366 data.pz[nx]= data.pz[i];
1367 data.t1x[nx]= data.t1x[i];
1368 data.t1y[nx]= -data.t1y[i];
1369 data.t1z[nx]= data.t1z[i];
1370 data.t2x[nx]= data.t2x[i];
1371 data.t2y[nx]= -data.t2y[i];
1372 data.t2z[nx]= data.t2z[i];
1373 data.psalp[nx]= -data.psalp[i];
1374 data.pbi[nx]= data.pbi[i];
1375
1376 } /* for( i = m2; i <= data.m; i++ ) */
1377
1378 data.m= data.m*2;
1379
1380 } /* if( data.m >= m2) */
1381
1382 } /* if( iy != 0) */
1383
1384 /* reflect along x axis */
1385 if( ix == 0 ) return( TRUE );
1386
1387 if( data.n > 0 )
1388 {
1389 /* Reallocate tags buffer */
1390 mreq = (size_t)(2 * data.n) * sizeof(int);
1391 mem_realloc( (void **)&data.itag, mreq, "in geometry.c" );
1392
1393 /* Reallocate wire buffers */
1394 mreq = (size_t)(2 * data.n) * sizeof(double);
1395 mem_realloc( (void **)&data.x1, mreq, "in geometry.c" );
1396 mem_realloc( (void **)&data.y1, mreq, "in geometry.c" );
1397 mem_realloc( (void **)&data.z1, mreq, "in geometry.c" );
1398 mem_realloc( (void **)&data.x2, mreq, "in geometry.c" );
1399 mem_realloc( (void **)&data.y2, mreq, "in geometry.c" );
1400 mem_realloc( (void **)&data.z2, mreq, "in geometry.c" );
1401 mem_realloc( (void **)&data.bi, mreq, "in geometry.c" );
1402
1403 for( i = 0; i < data.n; i++ )
1404 {
1405 nx= i+ data.n;
1406 e1= data.x1[i];
1407 e2= data.x2[i];
1408
1409 if( (fabs(e1)+fabs(e2) <= 1.0e-5) || (e1*e2 < -1.0e-6) )
1410 {
1411 fprintf( stderr,
1412 "xnec2c: reflc(): geometry data error\n"
1413 "segment %d lies in plane of symmetry\n", i+1 );
1414 stop( _("reflc(): Geometry data error\n"\
1415 "Segment lies in plane of symmetry"), ERR_OK );
1416 return( FALSE );
1417 }
1418
1419 data.x1[nx]= -e1;
1420 data.y1[nx]= data.y1[i];
1421 data.z1[nx]= data.z1[i];
1422 data.x2[nx]= -e2;
1423 data.y2[nx]= data.y2[i];
1424 data.z2[nx]= data.z2[i];
1425 itagi= data.itag[i];
1426
1427 if( itagi == 0)
1428 data.itag[nx]=0;
1429 if( itagi != 0)
1430 data.itag[nx]= itagi+ iti;
1431
1432 data.bi[nx]= data.bi[i];
1433 }
1434
1435 data.n= data.n*2;
1436
1437 } /* if( data.n > 0) */
1438
1439 if( data.m == 0 ) return( TRUE );
1440
1441 /* Reallocate patch buffers */
1442 mreq = (size_t)(2 * data.m) * sizeof(double);
1443 mem_realloc( (void **)&data.px, mreq, "in geometry.c" );
1444 mem_realloc( (void **)&data.py, mreq, "in geometry.c" );
1445 mem_realloc( (void **)&data.pz, mreq, "in geometry.c" );
1446 mem_realloc( (void **)&data.t1x, mreq, "in geometry.c" );
1447 mem_realloc( (void **)&data.t1y, mreq, "in geometry.c" );
1448 mem_realloc( (void **)&data.t1z, mreq, "in geometry.c" );
1449 mem_realloc( (void **)&data.t2x, mreq, "in geometry.c" );
1450 mem_realloc( (void **)&data.t2y, mreq, "in geometry.c" );
1451 mem_realloc( (void **)&data.t2z, mreq, "in geometry.c" );
1452 mem_realloc( (void **)&data.pbi, mreq, "in geometry.c" );
1453 mem_realloc( (void **)&data.psalp, mreq, "in geometry.c" );
1454
1455 for( i = 0; i < data.m; i++ )
1456 {
1457 nx= i+data.m;
1458 if( fabs( data.px[i]) <= 1.0e-10)
1459 {
1460 fprintf( stderr,
1461 "xnec2c: reflc(): geometry data error\n"
1462 "patch %d lies in plane of symmetry\n", i+1 );
1463 stop( _("reflc(): Geometry data error\n"\
1464 "Patch lies in plane of symmetry"), ERR_OK );
1465 return( FALSE );
1466 }
1467
1468 data.px[nx]= -data.px[i];
1469 data.py[nx]= data.py[i];
1470 data.pz[nx]= data.pz[i];
1471 data.t1x[nx]= -data.t1x[i];
1472 data.t1y[nx]= data.t1y[i];
1473 data.t1z[nx]= data.t1z[i];
1474 data.t2x[nx]= -data.t2x[i];
1475 data.t2y[nx]= data.t2y[i];
1476 data.t2z[nx]= data.t2z[i];
1477 data.psalp[nx]= -data.psalp[i];
1478 data.pbi[nx]= data.pbi[i];
1479 }
1480
1481 data.m= data.m*2;
1482 return( TRUE );
1483
1484 } /* if( ix >= 0) */
1485
1486 /* reproduce structure with rotation to form cylindrical structure */
1487 fnop= (double)nop;
1488 data.ipsym=-1;
1489 sam=TP/ fnop;
1490 cs= cos( sam);
1491 ss= sin( sam);
1492
1493 if( data.n > 0)
1494 {
1495 data.n *= nop;
1496 nx= data.np;
1497
1498 /* Reallocate tags buffer */
1499 mreq = (size_t)data.n * sizeof(int);
1500 mem_realloc( (void **)&data.itag, mreq, "in geometry.c" );
1501
1502 /* Reallocate wire buffers */
1503 mreq = (size_t)data.n * sizeof(double);
1504 mem_realloc( (void **)&data.x1, mreq, "in geometry.c" );
1505 mem_realloc( (void **)&data.y1, mreq, "in geometry.c" );
1506 mem_realloc( (void **)&data.z1, mreq, "in geometry.c" );
1507 mem_realloc( (void **)&data.x2, mreq, "in geometry.c" );
1508 mem_realloc( (void **)&data.y2, mreq, "in geometry.c" );
1509 mem_realloc( (void **)&data.z2, mreq, "in geometry.c" );
1510 mem_realloc( (void **)&data.bi, mreq, "in geometry.c" );
1511
1512 for( i = nx; i < data.n; i++ )
1513 {
1514 k= i- data.np;
1515 xk= data.x1[k];
1516 yk= data.y1[k];
1517 data.x1[i]= xk* cs- yk* ss;
1518 data.y1[i]= xk* ss+ yk* cs;
1519 data.z1[i]= data.z1[k];
1520 xk= data.x2[k];
1521 yk= data.y2[k];
1522 data.x2[i]= xk* cs- yk* ss;
1523 data.y2[i]= xk* ss+ yk* cs;
1524 data.z2[i]= data.z2[k];
1525 data.bi[i]= data.bi[k];
1526 itagi= data.itag[k];
1527
1528 if( itagi == 0)
1529 data.itag[i]=0;
1530 if( itagi != 0)
1531 data.itag[i]= itagi+ iti;
1532 }
1533
1534 } /* if( data.n >= n2) */
1535
1536 if( data.m == 0 ) return( TRUE );
1537
1538 data.m *= nop;
1539 nx= data.mp;
1540
1541 /* Reallocate patch buffers */
1542 mreq = (size_t)data.m * sizeof(double);
1543 mem_realloc( (void **)&data.px, mreq, "in geometry.c" );
1544 mem_realloc( (void **)&data.py, mreq, "in geometry.c" );
1545 mem_realloc( (void **)&data.pz, mreq, "in geometry.c" );
1546 mem_realloc( (void **)&data.t1x, mreq, "in geometry.c" );
1547 mem_realloc( (void **)&data.t1y, mreq, "in geometry.c" );
1548 mem_realloc( (void **)&data.t1z, mreq, "in geometry.c" );
1549 mem_realloc( (void **)&data.t2x, mreq, "in geometry.c" );
1550 mem_realloc( (void **)&data.t2y, mreq, "in geometry.c" );
1551 mem_realloc( (void **)&data.t2z, mreq, "in geometry.c" );
1552 mem_realloc( (void **)&data.pbi, mreq, "in geometry.c" );
1553 mem_realloc( (void **)&data.psalp, mreq, "in geometry.c" );
1554
1555 for( i = nx; i < data.m; i++ )
1556 {
1557 k = i-data.mp;
1558 xk= data.px[k];
1559 yk= data.py[k];
1560 data.px[i]= xk* cs- yk* ss;
1561 data.py[i]= xk* ss+ yk* cs;
1562 data.pz[i]= data.pz[k];
1563 xk= data.t1x[k];
1564 yk= data.t1y[k];
1565 data.t1x[i]= xk* cs- yk* ss;
1566 data.t1y[i]= xk* ss+ yk* cs;
1567 data.t1z[i]= data.t1z[k];
1568 xk= data.t2x[k];
1569 yk= data.t2y[k];
1570 data.t2x[i]= xk* cs- yk* ss;
1571 data.t2y[i]= xk* ss+ yk* cs;
1572 data.t2z[i]= data.t2z[k];
1573 data.psalp[i]= data.psalp[k];
1574 data.pbi[i]= data.pbi[k];
1575
1576 } /* for( i = nx; i < data.m; i++ ) */
1577
1578 return( TRUE );
1579 }
1580
1581 /*-----------------------------------------------------------------------*/
1582
1583 /* subroutine wire generates segment geometry */
1584 /* data for a straight wire of ns segments. */
1585 void
wire(double xw1,double yw1,double zw1,double xw2,double yw2,double zw2,double rad,double rdel,double rrad,int ns,int itg)1586 wire( double xw1, double yw1, double zw1,
1587 double xw2, double yw2, double zw2, double rad,
1588 double rdel, double rrad, int ns, int itg )
1589 {
1590 int ist, i;
1591 size_t mreq;
1592 double xd, yd, zd, delz, rd, fns, radz;
1593 double xs1, ys1, zs1, xs2, ys2, zs2;
1594
1595 if( ns < 1) return;
1596
1597 ist= data.n;
1598 data.n= data.n+ ns;
1599 data.np= data.n;
1600 data.mp= data.m;
1601 data.ipsym=0;
1602
1603 /* Reallocate tags buffer */
1604 mreq = (size_t)data.n * sizeof(int);
1605 mem_realloc( (void **)&data.itag, mreq, "in geometry.c" );
1606
1607 /* Reallocate wire buffers */
1608 mreq = (size_t)data.n * sizeof(double);
1609 mem_realloc( (void **)&data.x1, mreq, "in geometry.c" );
1610 mem_realloc( (void **)&data.y1, mreq, "in geometry.c" );
1611 mem_realloc( (void **)&data.z1, mreq, "in geometry.c" );
1612 mem_realloc( (void **)&data.x2, mreq, "in geometry.c" );
1613 mem_realloc( (void **)&data.y2, mreq, "in geometry.c" );
1614 mem_realloc( (void **)&data.z2, mreq, "in geometry.c" );
1615 mem_realloc( (void **)&data.bi, mreq, "in geometry.c" );
1616
1617 xd= xw2- xw1;
1618 yd= yw2- yw1;
1619 zd= zw2- zw1;
1620
1621 if( fabs( rdel-1.0) >= 1.0e-6)
1622 {
1623 delz= sqrt( xd* xd+ yd* yd+ zd* zd);
1624 xd= xd/ delz;
1625 yd= yd/ delz;
1626 zd= zd/ delz;
1627 delz= delz*(1.0- rdel)/(1.0- pow(rdel, ns) );
1628 rd= rdel;
1629 }
1630 else
1631 {
1632 fns= ns;
1633 xd= xd/ fns;
1634 yd= yd/ fns;
1635 zd= zd/ fns;
1636 delz=1.0;
1637 rd=1.0;
1638 }
1639
1640 radz= rad;
1641 xs1= xw1;
1642 ys1= yw1;
1643 zs1= zw1;
1644
1645 for( i = ist; i < data.n; i++ )
1646 {
1647 data.itag[i]= itg;
1648 xs2= xs1+ xd* delz;
1649 ys2= ys1+ yd* delz;
1650 zs2= zs1+ zd* delz;
1651 data.x1[i]= xs1;
1652 data.y1[i]= ys1;
1653 data.z1[i]= zs1;
1654 data.x2[i]= xs2;
1655 data.y2[i]= ys2;
1656 data.z2[i]= zs2;
1657 data.bi[i]= radz;
1658 delz= delz* rd;
1659 radz= radz* rrad;
1660 xs1= xs2;
1661 ys1= ys2;
1662 zs1= zs2;
1663 }
1664
1665 data.x2[data.n-1]= xw2;
1666 data.y2[data.n-1]= yw2;
1667 data.z2[data.n-1]= zw2;
1668
1669 return;
1670 }
1671
1672 /*-----------------------------------------------------------------------*/
1673
1674