WebSVN – Nec2c – Blame – /trunk/src/geometry.c

Rev	Author	Line No.	Line
2	mjames	1	/* Translated to the C language by N. Kyriazis 20 Aug 2003 *
		2
		3	Program NEC(input,tape5=input,output,tape11,tape12,tape13,tape14,
		4	tape15,tape16,tape20,tape21)
		5
		6	Numerical Electromagnetics Code (NEC2) developed at Lawrence
		7	Livermore lab., Livermore, CA. (contact G. Burke at 415-422-8414
		8	for problems with the NEC code. For problems with the vax implem-
		9	entation, contact J. Breakall at 415-422-8196 or E. Domning at 415
		10	422-5936)
		11	file created 4/11/80.
		12
		13	*********Notice********
		14	This computer code material was prepared as an account of work
		15	sponsored by the United States government. Neither the United
		16	States nor the United States Department Of Energy, nor any of
		17	their employees, nor any of their contractors, subcontractors,
		18	or their employees, makes any warranty, express or implied, or
		19	assumes any legal liability or responsibility for the accuracy,
		20	completeness or usefulness of any information, apparatus, product
		21	or process disclosed, or represents that its use would not infringe
		22	privately-owned rights.
		23
		24	*******************************************************************/
		25
		26	#include "nec2c.h"
		27
		28	/* common /data/ */
		29	extern data_t data;
		30
		31	/* common /segj/ */
		32	extern segj_t segj;
		33
		34	/* pointers to input/output files */
		35	extern FILE input_fp, output_fp, *plot_fp;
		36
		37	/* common /plot/ */
		38	extern plot_t plot;
		39
		40	/-------------------------------------------------------------------/
		41
		42	/* arc generates segment geometry data for an arc of ns segments */
		43	void arc (int itg, int ns, double rada, double ang1, double ang2, double rad)
		44	{
		45	int ist, i, mreq;
		46	double ang, dang, xs1, xs2, zs1, zs2;
		47
		48	ist = data.n;
		49	data.n += ns;
		50	data.np = data.n;
		51	data.mp = data.m;
		52	data.ipsym = 0;
		53
		54	if (ns < 1)
		55	return;
		56
		57	if (fabsl (ang2 - ang1) < 360.00001)
		58	{
		59	/* Reallocate tags buffer */
		60	mem_realloc ((void ) &data.itag, (data.n + data.m) sizeof (int));
		61
		62	/* Reallocate wire buffers */
		63	mreq = data.n * sizeof (double);
		64	mem_realloc ((void *) &data.x1, mreq);
		65	mem_realloc ((void *) &data.y1, mreq);
		66	mem_realloc ((void *) &data.z1, mreq);
		67	mem_realloc ((void *) &data.x2, mreq);
		68	mem_realloc ((void *) &data.y2, mreq);
		69	mem_realloc ((void *) &data.z2, mreq);
		70	mem_realloc ((void *) &data.bi, mreq);
		71
		72	ang = ang1 * TA;
		73	dang = (ang2 - ang1) * TA / ns;
		74	xs1 = rada * cosl (ang);
		75	zs1 = rada * sinl (ang);
		76
		77	for (i = ist; i < data.n; i++)
		78	{
		79	ang += dang;
		80	xs2 = rada * cosl (ang);
		81	zs2 = rada * sinl (ang);
		82	data.x1[i] = xs1;
		83
		84	data.y1[i] = 0.;
		85	data.z1[i] = zs1;
		86	data.x2[i] = xs2;
		87	data.y2[i] = 0.;
		88	data.z2[i] = zs2;
		89	xs1 = xs2;
		90	zs1 = zs2;
		91	data.bi[i] = rad;
		92	data.itag[i] = itg;
		93
		94	} /* for( i = ist; i < data.n; i++ ) */
		95
		96	} /* if( fabsl( ang2- ang1) < 360.00001) */
		97	else
		98	{
		99	fprintf (output_fp, "\n ERROR -- ARC ANGLE EXCEEDS 360 DEGREES");
		100	stop (-1);
		101	}
		102
		103	return;
		104	}
		105
		106	/-----------------------------------------------------------------------/
		107
		108	/* connect sets up segment connection data in arrays icon1 and */
		109	/* icon2 by searching for segment ends that are in contact. */
		110	void conect (int ignd)
		111	{
		112	int i, iz, ic, j, jx, ix, ixx, iseg, iend, jend, jump, ipf;
		113	// int nsflg;
		114	double sep = 0., xi1, yi1, zi1, xi2, yi2, zi2;
		115	double slen, xa, ya, za, xs, ys, zs;
		116
		117	segj.maxcon = 1;
		118
		119	if (ignd != 0)
		120	{
		121	fprintf (output_fp, "\n\n GROUND PLANE SPECIFIED.\n");
		122
		123	if (ignd > 0)
		124	fprintf (
		125	output_fp,
		126	"\n WHERE WIRE ENDS TOUCH GROUND, CURRENT WILL"
		127	" BE INTERPOLATED TO IMAGE IN GROUND PLANE.\n");
		128
		129	if (data.ipsym == 2)
		130	{
		131	data.np = 2 * data.np;
		132	data.mp = 2 * data.mp;
		133	}
		134
		135	if (abs (data.ipsym) > 2)
		136	{
		137	data.np = data.n;
		138	data.mp = data.m;
		139	}
		140
		141	/* possibly should be error condition?? /
		142	if (data.np > data.n)
		143	{
		144	fprintf (output_fp, "\n ERROR: NP > N IN CONECT()");
		145	stop (-1);
		146	}
		147
		148	if ((data.np == data.n) && (data.mp == data.m))
		149	data.ipsym = 0;
		150
		151	} /* if( ignd != 0) */
		152
		153	if (data.n != 0)
		154	{
		155	/* Allocate memory to connections */
		156	mem_realloc ((void ) &data.icon1, (data.n + data.m) sizeof (int));
		157	mem_realloc ((void ) &data.icon2, (data.n + data.m) sizeof (int));
		158
		159	for (i = 0; i < data.n; i++)
		160	{
		161	data.icon1[i] = data.icon2[i] = 0;
		162	iz = i + 1;
		163	xi1 = data.x1[i];
		164	yi1 = data.y1[i];
		165	zi1 = data.z1[i];
		166	xi2 = data.x2[i];
		167	yi2 = data.y2[i];
		168	zi2 = data.z2[i];
		169	slen = sqrtl (
		170	(xi2 - xi1) * (xi2 - xi1) + (yi2 - yi1) * (yi2 - yi1) +
		171	(zi2 - zi1) * (zi2 - zi1)) *
		172	SMIN;
		173
		174	/* determine connection data for end 1 of segment. */
		175	jump = FALSE;
		176	if (ignd > 0)
		177	{
		178	if (zi1 <= -slen)
		179	{
		180	fprintf (
		181	output_fp,
		182	"\n GEOMETRY DATA ERROR -- SEGMENT"
		183	" %d EXTENDS BELOW GROUND",
		184	iz);
		185	stop (-1);
		186	}
		187
		188	if (zi1 <= slen)
		189	{
		190	data.icon1[i] = iz;
		191	data.z1[i] = 0.;
		192	jump = TRUE;
		193
		194	} /* if( zi1 <= slen) */
		195
		196	} /* if( ignd > 0) */
		197
		198	if (!jump)
		199	{
		200	ic = i;
		201	for (j = 1; j < data.n; j++)
		202	{
		203	ic++;
		204	if (ic >= data.n)
		205	ic = 0;
		206
		207	sep = fabsl (xi1 - data.x1[ic]) +
		208	fabsl (yi1 - data.y1[ic]) +
		209	fabsl (zi1 - data.z1[ic]);
		210	if (sep <= slen)
		211	{
		212	data.icon1[i] = -(ic + 1);
		213	break;
		214	}
		215
		216	sep = fabsl (xi1 - data.x2[ic]) +
		217	fabsl (yi1 - data.y2[ic]) +
		218	fabsl (zi1 - data.z2[ic]);
		219	if (sep <= slen)
		220	{
		221	data.icon1[i] = (ic + 1);
		222	break;
		223	}
		224
		225	} /* for( j = 1; j < data.n; j++) */
		226
		227	} /* if( ! jump ) */
		228
		229	/* determine connection data for end 2 of segment. */
		230	if ((ignd > 0) \|\| jump)
		231	{
		232	if (zi2 <= -slen)
		233	{
		234	fprintf (
		235	output_fp,
		236	"\n GEOMETRY DATA ERROR -- SEGMENT"
		237	" %d EXTENDS BELOW GROUND",
		238	iz);
		239	stop (-1);
		240	}
		241
		242	if (zi2 <= slen)
		243	{
		244	if (data.icon1[i] == iz)
		245	{
		246	fprintf (
		247	output_fp,
		248	"\n GEOMETRY DATA ERROR -- SEGMENT"
		249	" %d LIES IN GROUND PLANE",
		250	iz);
		251	stop (-1);
		252	}
		253
		254	data.icon2[i] = iz;
		255	data.z2[i] = 0.;
		256	continue;
		257
		258	} /* if( zi2 <= slen) */
		259
		260	} /* if( ignd > 0) */
		261
		262	ic = i;
		263	for (j = 1; j < data.n; j++)
		264	{
		265	ic++;
		266	if (ic >= data.n)
		267	ic = 0;
		268
		269	sep = fabsl (xi2 - data.x1[ic]) + fabsl (yi2 - data.y1[ic]) +
		270	fabsl (zi2 - data.z1[ic]);
		271	if (sep <= slen)
		272	{
		273	data.icon2[i] = (ic + 1);
		274	break;
		275	}
		276
		277	sep = fabsl (xi2 - data.x2[ic]) + fabsl (yi2 - data.y2[ic]) +
		278	fabsl (zi2 - data.z2[ic]);
		279	if (sep <= slen)
		280	{
		281	data.icon2[i] = -(ic + 1);
		282	break;
		283	}
		284
		285	} /* for( j = 1; j < data.n; j++ ) */
		286
		287	} /* for( i = 0; i < data.n; i++ ) */
		288
		289	/* find wire-surface connections for new patches */
		290	if (data.m != 0)
		291	{
		292	ix = -1;
		293	i = 0;
		294	while (++i <= data.m)
		295	{
		296	ix++;
		297	xs = data.px[ix];
		298	ys = data.py[ix];
		299	zs = data.pz[ix];
		300
		301	for (iseg = 0; iseg < data.n; iseg++)
		302	{
		303	xi1 = data.x1[iseg];
		304	yi1 = data.y1[iseg];
		305	zi1 = data.z1[iseg];
		306	xi2 = data.x2[iseg];
		307	yi2 = data.y2[iseg];
		308	zi2 = data.z2[iseg];
		309
		310	/* for first end of segment */
		311	slen = (fabsl (xi2 - xi1) + fabsl (yi2 - yi1) +
		312	fabsl (zi2 - zi1)) *
		313	SMIN;
		314	sep = fabsl (xi1 - xs) + fabsl (yi1 - ys) +
		315	fabsl (zi1 - zs);
		316
		317	/* connection - divide patch into 4 patches at present
		318	* array loc. */
		319	if (sep <= slen)
		320	{
		321	data.icon1[iseg] = PCHCON + i;
		322	ic = 0;
		323	subph (i, ic);
		324	break;
		325	}
		326
		327	sep = fabsl (xi2 - xs) + fabsl (yi2 - ys) +
		328	fabsl (zi2 - zs);
		329	if (sep <= slen)
		330	{
		331	data.icon2[iseg] = PCHCON + i;
		332	ic = 0;
		333	subph (i, ic);
		334	break;
		335	}
		336
		337	} /* for( iseg = 0; iseg < data.n; iseg++ ) */
		338
		339	} /* while( ++i <= data.m ) */
		340
		341	} /* if( data.m != 0) */
		342
		343	} /* if( data.n != 0) */
		344
		345	fprintf (
		346	output_fp,
		347	"\n\n"
		348	" TOTAL SEGMENTS USED=%5d "
		349	"NO. SEG. IN A SYMMETRIC CELL=%5d"
		350	" SYMMETRY FLAG= %d\n",
		351	data.n,
		352	data.np,
		353	data.ipsym);
		354
		355	if (data.m > 0)
		356	fprintf (
		357	output_fp,
		358	"\n"
		359	" TOTAL PATCHES USED: %d PATCHES"
		360	" IN A SYMMETRIC CELL: %d",
		361	data.m,
		362	data.mp);
		363
		364	iseg = (data.n + data.m) / (data.np + data.mp);
		365	if (iseg != 1)
		366	{
		367	/* may be error condition?? */
		368	if (data.ipsym == 0)
		369	{
		370	fprintf (output_fp, "\n ERROR: IPSYM=0 IN CONECT()");
		371	stop (-1);
		372	}
		373
		374	if (data.ipsym < 0)
		375	fprintf (
		376	output_fp,
		377	"\n STRUCTURE HAS %d FOLD ROTATIONAL SYMMETRY\n",
		378	iseg);
		379	else
		380	{
		381	ic = iseg / 2;
		382	if (iseg == 8)
		383	ic = 3;
		384	fprintf (output_fp, "\n STRUCTURE HAS %d PLANES OF SYMMETRY\n", ic);
		385	} /* if( data.ipsym < 0 ) */
		386
		387	} /* if( iseg == 1) */
		388
		389	if (data.n == 0)
		390	return;
		391
		392	/* Allocate to connection buffers */
		393	mem_realloc ((void ) &segj.jco, segj.maxcon sizeof (int));
		394
		395	/* adjust connected seg. ends to exactly coincide. print junctions */
		396	/* of 3 or more seg. also find old seg. connecting to new seg. */
		397	iseg = 0;
		398	ipf = FALSE;
		399	for (j = 0; j < data.n; j++)
		400	{
		401	jx = j + 1;
		402	iend = -1;
		403	jend = -1;
		404	ix = data.icon1[j];
		405	ic = 1;
		406	segj.jco[0] = -jx;
		407	xa = data.x1[j];
		408	ya = data.y1[j];
		409	za = data.z1[j];
		410
		411	while (TRUE)
		412	{
		413	if ((ix != 0) && (ix != (j + 1)) && (ix <= PCHCON))
		414	{
		415	// nsflg=0;
		416
		417	do
		418	{
		419	if (ix == 0)
		420	{
		421	fprintf (
		422	output_fp,
		423	"\n CONNECT - SEGMENT CONNECTION ERROR "
		424	"FOR SEGMENT: %d",
		425	ix);
		426	stop (-1);
		427	}
		428
		429	if (ix < 0)
		430	ix = -ix;
		431	else
		432	jend = -jend;
		433
		434	jump = FALSE;
		435
		436	if (ix == jx)
		437	break;
		438
		439	if (ix < jx)
		440	{
		441	jump = TRUE;
		442	break;
		443	}
		444
		445	/* Record max. no. of connections */
		446	ic++;
		447	if (ic >= segj.maxcon)
		448	{
		449	segj.maxcon = ic + 1;
		450	mem_realloc (
		451	(void *) &segj.jco,
		452	segj.maxcon * sizeof (int));
		453	}
		454	segj.jco[ic - 1] = ix * jend;
		455
		456	/*if( ix > 0){
		457	//nsflg=1;
		458	}*/
		459	ixx = ix - 1;
		460	if (jend != 1)
		461	{
		462	xa = xa + data.x1[ixx];
		463	ya = ya + data.y1[ixx];
		464	za = za + data.z1[ixx];
		465	ix = data.icon1[ixx];
		466	continue;
		467	}
		468
		469	xa = xa + data.x2[ixx];
		470	ya = ya + data.y2[ixx];
		471	za = za + data.z2[ixx];
		472	ix = data.icon2[ixx];
		473
		474	} /* do */
		475	while (ix != 0);
		476
		477	if (jump && (iend == 1))
		478	break;
		479	else if (jump)
		480	{
		481	iend = 1;
		482	jend = 1;
		483	ix = data.icon2[j];
		484	ic = 1;
		485	segj.jco[0] = jx;
		486	xa = data.x2[j];
		487	ya = data.y2[j];
		488	za = data.z2[j];
		489	continue;
		490	}
		491
		492	sep = (double) ic;
		493	xa = xa / sep;
		494	ya = ya / sep;
		495	za = za / sep;
		496
		497	for (i = 0; i < ic; i++)
		498	{
		499	ix = segj.jco[i];
		500	if (ix <= 0)
		501	{
		502	ix = -ix;
		503	ixx = ix - 1;
		504	data.x1[ixx] = xa;
		505	data.y1[ixx] = ya;
		506	data.z1[ixx] = za;
		507	continue;
		508	}
		509
		510	ixx = ix - 1;
		511	data.x2[ixx] = xa;
		512	data.y2[ixx] = ya;
		513	data.z2[ixx] = za;
		514
		515	} /* for( i = 0; i < ic; i++ ) */
		516
		517	if (!ipf)
		518	{
		519	fprintf (
		520	output_fp,
		521	"\n\n"
		522	" - MULTIPLE WIRE JUNCTIONS -\n"
		523	" JUNCTION SEGMENTS (- FOR END 1, + FOR END "
		524	"2)");
		525
		526	ipf = TRUE;
		527	}
		528	if (ic >= 3)
		529	{
		530	iseg++;
		531	fprintf (output_fp, "\n %5d ", iseg);
		532
		533	for (i = 1; i <= ic; i++)
		534	{
		535	fprintf (output_fp, " %4d", segj.jco[i - 1]);
		536	if (!(i % 20))
		537	fprintf (
		538	output_fp, "\n ");
		539	}
		540	} /* if( ic >= 3) */
		541
		542	} /if( (ix != 0) && (ix != j) && (ix <= PCHCON) ) /
		543
		544	if (iend == 1)
		545	break;
		546
		547	iend = 1;
		548	jend = 1;
		549	ix = data.icon2[j];
		550	ic = 1;
		551	segj.jco[0] = jx;
		552	xa = data.x2[j];
		553	ya = data.y2[j];
		554	za = data.z2[j];
		555
		556	} /* while( TRUE ) */
		557
		558	} /* for( j = 0; j < data.n; j++ ) */
		559	/* print NONE if nothing printed */
		560	if (ipf == 1 && iseg == 0)
		561	{
		562	fprintf (output_fp, "\n NONE");
		563	}
		564
		565	mem_realloc ((void ) &segj.ax, segj.maxcon sizeof (double));
		566	mem_realloc ((void ) &segj.bx, segj.maxcon sizeof (double));
		567	mem_realloc ((void ) &segj.cx, segj.maxcon sizeof (double));
		568
		569	return;
		570	}
		571
		572	/-----------------------------------------------------------------------/
		573
		574	/* datagn is the main routine for input of geometry data. */
		575	void datagn (void)
		576	{
		577	char gm[3];
		578	char ifx[2] = {'', 'X'}, ify[2] = {'', 'Y'}, ifz[2] = {'*', 'Z'};
		579	char ipt[4] = {'P', 'R', 'T', 'Q'};
		580
		581	/* input card mnemonic list */
		582	/* "XT" stands for "exit", added for testing */
		583	#define GM_NUM 12
		584	char *atst[GM_NUM] = {
		585	"GW", "GX", "GR", "GS", "GE", "GM", "SP", "SM", "GA", "SC", "GH", "GF"};
		586
		587	int nwire, isct, iphd, i1, i2, itg, iy, iz, mreq;
		588	int ix, i, ns, gm_num; /* geometry card id as a number */
		589	double rad, xs1, xs2, ys1, ys2, zs1, zs2, x4 = 0, y4 = 0, z4 = 0;
		590	double x3 = 0, y3 = 0, z3 = 0, xw1, xw2, yw1, yw2, zw1, zw2;
		591	double dummy;
		592
		593	data.ipsym = 0;
		594	nwire = 0;
		595	data.n = 0;
		596	data.np = 0;
		597	data.m = 0;
		598	data.mp = 0;
		599	isct = 0;
		600	iphd = FALSE;
		601
		602	/* read geometry data card and branch to */
		603	/* section for operation requested */
		604	do
		605	{
		606	readgm (gm, &itg, &ns, &xw1, &yw1, &zw1, &xw2, &yw2, &zw2, &rad);
		607
		608	/* identify card id mnemonic */
		609	for (gm_num = 0; gm_num < GM_NUM; gm_num++)
		610	if (strncmp (gm, atst[gm_num], 2) == 0)
		611	break;
		612
		613	if (iphd == FALSE)
		614	{
		615	fprintf (
		616	output_fp,
		617	"\n\n\n\n"
		618	" "
		619	"- - - STRUCTURE SPECIFICATION - - -\n\n"
		620	" "
		621	"COORDINATES MUST BE INPUT IN\n"
		622	" "
		623	"METERS OR BE SCALED TO METERS\n"
		624	" "
		625	"BEFORE STRUCTURE INPUT IS ENDED\n");
		626
		627	fprintf (
		628	output_fp,
		629	"\n\n"
		630	" WIRE "
		631	" NO. OF FIRST LAST TAG\n"
		632	" NO. X1 Y1 Z1 X2 Y2 "
		633	" Z2 RADIUS SEG. SEG. SEG. NO.");
		634
		635	iphd = 1;
		636	}
		637
		638	if (gm_num != 10)
		639	isct = 0;
		640
		641	switch (gm_num)
		642	{
		643	case 0: /* "gw" card, generate segment data for straight wire. */
		644
		645	nwire++;
		646	i1 = data.n + 1;
		647	i2 = data.n + ns;
		648
		649	fprintf (
		650	output_fp,
		651	"\n"
		652	" %5d %10.5f %10.5f %10.5f %10.5f"
		653	" %10.5f %10.5f %10.5f %5d %5d %5d %5d",
		654	nwire,
		655	xw1,
		656	yw1,
		657	zw1,
		658	xw2,
		659	yw2,
		660	zw2,
		661	rad,
		662	ns,
		663	i1,
		664	i2,
		665	itg);
		666
		667	if (rad != 0)
		668	{
		669	xs1 = 1.;
		670	ys1 = 1.;
		671	}
		672	else
		673	{
		674	readgm (
		675	gm,
		676	&ix,
		677	&iy,
		678	&xs1,
		679	&ys1,
		680	&zs1,
		681	&dummy,
		682	&dummy,
		683	&dummy,
		684	&dummy);
		685
		686	if (strcmp (gm, "GC") != 0)
		687	{
		688	fprintf (output_fp, "\n GEOMETRY DATA CARD ERROR");
		689	stop (-1);
		690	}
		691
		692	fprintf (
		693	output_fp,
		694	"\n ABOVE WIRE IS TAPERED. SEGMENT LENGTH RATIO: %9.5f\n"
		695	" "
		696	"RADIUS FROM: %9.5f TO: %9.5f",
		697	xs1,
		698	ys1,
		699	zs1);
		700
		701	if ((ys1 == 0) \|\| (zs1 == 0))
		702	{
		703	fprintf (output_fp, "\n GEOMETRY DATA CARD ERROR");
		704	stop (-1);
		705	}
		706
		707	rad = ys1;
		708	ys1 = powl ((zs1 / ys1), (1. / (ns - 1.)));
		709	}
		710
		711	wire (xw1, yw1, zw1, xw2, yw2, zw2, rad, xs1, ys1, ns, itg);
		712
		713	continue;
		714
		715	/* reflect structure along x,y, or z */
		716	/* axes or rotate to form cylinder. */
		717	case 1: /* "gx" card */
		718
		719	iy = ns / 10;
		720	iz = ns - iy * 10;
		721	ix = iy / 10;
		722	iy = iy - ix * 10;
		723
		724	if (ix != 0)
		725	ix = 1;
		726	if (iy != 0)
		727	iy = 1;
		728	if (iz != 0)
		729	iz = 1;
		730
		731	fprintf (
		732	output_fp,
		733	"\n STRUCTURE REFLECTED ALONG THE AXES %c %c %c"
		734	" - TAGS INCREMENTED BY %d\n",
		735	ifx[ix],
		736	ify[iy],
		737	ifz[iz],
		738	itg);
		739
		740	reflc (ix, iy, iz, itg, ns);
		741
		742	continue;
		743
		744	case 2: /* "gr" card */
		745
		746	fprintf (
		747	output_fp,
		748	"\n STRUCTURE ROTATED ABOUT Z-AXIS %d TIMES"
		749	" - LABELS INCREMENTED BY %d\n",
		750	ns,
		751	itg);
		752
		753	ix = -1;
		754	iz = 0;
		755	reflc (ix, iy, iz, itg, ns);
		756
		757	continue;
		758
		759	case 3: /* "gs" card, scale structure dimensions by factor xw1. */
		760
		761	if (data.n > 0)
		762	{
		763	for (i = 0; i < data.n; i++)
		764	{
		765	data.x1[i] = data.x1[i] * xw1;
		766	data.y1[i] = data.y1[i] * xw1;
		767	data.z1[i] = data.z1[i] * xw1;
		768	data.x2[i] = data.x2[i] * xw1;
		769	data.y2[i] = data.y2[i] * xw1;
		770	data.z2[i] = data.z2[i] * xw1;
		771	data.bi[i] = data.bi[i] * xw1;
		772	}
		773	} /* if( data.n >= n2) */
		774
		775	if (data.m > 0)
		776	{
		777	yw1 = xw1 * xw1;
		778	for (i = 0; i < data.m; i++)
		779	{
		780	data.px[i] = data.px[i] * xw1;
		781	data.py[i] = data.py[i] * xw1;
		782	data.pz[i] = data.pz[i] * xw1;
		783	data.pbi[i] = data.pbi[i] * yw1;
		784	}
		785	} /* if( data.m >= m2) */
		786
		787	fprintf (output_fp, "\n STRUCTURE SCALED BY FACTOR %10.5f", xw1);
		788
		789	continue;
		790
		791	case 4: /* "ge" card, terminate structure geometry input. */
		792
		793	if (ns != 0)
		794	{
		795	plot.iplp1 = 1;
		796	plot.iplp2 = 1;
		797	}
		798
		799	conect (itg);
		800
		801	if (data.n != 0)
		802	{
		803	/* Allocate wire buffers */
		804	mreq = data.n * sizeof (double);
		805	mem_realloc ((void *) &data.si, mreq);
		806	mem_realloc ((void *) &data.sab, mreq);
		807	mem_realloc ((void *) &data.cab, mreq);
		808	mem_realloc ((void *) &data.salp, mreq);
		809	mem_realloc ((void *) &data.x, mreq);
		810	mem_realloc ((void *) &data.y, mreq);
		811	mem_realloc ((void *) &data.z, mreq);
		812
		813	fprintf (
		814	output_fp,
		815	"\n\n\n\n\n"
		816	" "
		817	" - - - - SEGMENTATION DATA - - - -\n\n"
		818	" "
		819	" COORDINATES IN METERS\n\n"
		820	" "
		821	" I+ AND I- INDICATE THE SEGMENTS BEFORE AND AFTER I\n\n");
		822
		823	fprintf (
		824	output_fp,
		825	"\n"
		826	" SEG. COORDINATES OF SEG. CENTER SEG. "
		827	"ORIENTATION ANGLES WIRE CONNECTION DATA TAG\n"
		828	" NO. X Y Z LENGTH "
		829	" ALPHA BETA RADIUS I- I I+ NO.");
		830
		831	for (i = 0; i < data.n; i++)
		832	{
		833	xw1 = data.x2[i] - data.x1[i];
		834	yw1 = data.y2[i] - data.y1[i];
		835	zw1 = data.z2[i] - data.z1[i];
		836	data.x[i] = (data.x1[i] + data.x2[i]) / 2.;
		837	data.y[i] = (data.y1[i] + data.y2[i]) / 2.;
		838	data.z[i] = (data.z1[i] + data.z2[i]) / 2.;
		839	xw2 = xw1 * xw1 + yw1 * yw1 + zw1 * zw1;
		840	yw2 = sqrtl (xw2);
		841	yw2 = (xw2 / yw2 + yw2) * .5;
		842	data.si[i] = yw2;
		843	data.cab[i] = xw1 / yw2;
		844	data.sab[i] = yw1 / yw2;
		845	xw2 = zw1 / yw2;
		846
		847	if (xw2 > 1.)
		848	xw2 = 1.;
		849	if (xw2 < -1.)
		850	xw2 = -1.;
		851
		852	data.salp[i] = xw2;
		853	xw2 = asinl (xw2) * TD;
		854	yw2 = atan2l (yw1, xw1) * TD;
		855
		856	fprintf (
		857	output_fp,
		858	"\n"
		859	" %5d%10.5f%10.5f%10.5f%10.5f"
		860	" %10.5f%10.5f%10.5f%7d%7d%7d%7d",
		861	i + 1,
		862	data.x[i],
		863	data.y[i],
		864	data.z[i],
		865	data.si[i],
		866	xw2,
		867	yw2,
		868	data.bi[i],
		869	data.icon1[i],
		870	i + 1,
		871	data.icon2[i],
		872	data.itag[i]);
		873
		874	if (plot.iplp1 == 1)
		875	fprintf (
		876	plot_fp,
		877	"%12.5E %12.5E %12.5E "
		878	"%12.5E %12.5E %12.5E %12.5E%7d%7d%7d\n",
		879	data.x[i],
		880	data.y[i],
		881	data.z[i],
		882	data.si[i],
		883	xw2,
		884	yw2,
		885	data.bi[i],
		886	data.icon1[i],
		887	i + 1,
		888	data.icon2[i]);
		889
		890	if ((data.si[i] <= 1.e-20) \|\| (data.bi[i] <= 0.))
		891	{
		892	fprintf (output_fp, "\n SEGMENT DATA ERROR");
		893	stop (-1);
		894	}
		895
		896	} /* for( i = 0; i < data.n; i++ ) */
		897
		898	} /* if( data.n != 0) */
		899
		900	fprintf (output_fp, "\n\n\n");
		901	if (data.m != 0)
		902	{
		903	fprintf (
		904	output_fp,
		905	"\n\n\n"
		906	" "
		907	" - - - - - SURFACE PATCH DATA - - - -\n"
		908	" "
		909	" COORDINATES IN METERS\n\n"
		910	" PATCH COORD. OF PATCH CENTER UNIT NORMAL "
		911	"VECTOR "
		912	" PATCH COMPONENTS OF UNIT TANGENT VECTORS\n"
		913	" No: X Y Z X Y "
		914	" Z "
		915	" AREA X1 Y1 Z1 X2 Y2 "
		916	" Z2");
		917
		918	for (i = 0; i < data.m; i++)
		919	{
		920	xw1 = (data.t1y[i] * data.t2z[i] -
		921	data.t1z[i] * data.t2y[i]) *
		922	data.psalp[i];
		923	yw1 = (data.t1z[i] * data.t2x[i] -
		924	data.t1x[i] * data.t2z[i]) *
		925	data.psalp[i];
		926	zw1 = (data.t1x[i] * data.t2y[i] -
		927	data.t1y[i] * data.t2x[i]) *
		928	data.psalp[i];
		929
		930	fprintf (
		931	output_fp,
		932	"\n"
		933	" %4d %10.5f %10.5f %10.5f %8.5f %8.5f %8.5f"
		934	" %10.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f",
		935	i + 1,
		936	data.px[i],
		937	data.py[i],
		938	data.pz[i],
		939	xw1,
		940	yw1,
		941	zw1,
		942	data.pbi[i],
		943	data.t1x[i],
		944	data.t1y[i],
		945	data.t1z[i],
		946	data.t2x[i],
		947	data.t2y[i],
		948	data.t2z[i]);
		949
		950	} /* for( i = 0; i < data.m; i++ ) */
		951
		952	} /* if( data.m == 0) */
		953
		954	data.npm = data.n + data.m;
		955	data.np2m = data.n + 2 * data.m;
		956	data.np3m = data.n + 3 * data.m;
		957
		958	return;
		959
		960	/* "gm" card, move structure or reproduce */
		961	/* original structure in new positions. */
		962	case 5:
		963
		964	fprintf (
		965	output_fp,
		966	"\n THE STRUCTURE HAS BEEN MOVED, MOVE DATA CARD IS:\n"
		967	" %3d %5d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
		968	itg,
		969	ns,
		970	xw1,
		971	yw1,
		972	zw1,
		973	xw2,
		974	yw2,
		975	zw2,
		976	rad);
		977
		978	xw1 = xw1 * TA;
		979	yw1 = yw1 * TA;
		980	zw1 = zw1 * TA;
		981
		982	move (xw1, yw1, zw1, xw2, yw2, zw2, (int) (rad + .5), ns, itg);
		983	continue;
		984
		985	case 6: /* "sp" card, generate single new patch */
		986
		987	i1 = data.m + 1;
		988	ns++;
		989
		990	if (itg != 0)
		991	{
		992	fprintf (output_fp, "\n PATCH DATA ERROR");
		993	stop (-1);
		994	}
		995
		996	fprintf (
		997	output_fp,
		998	"\n"
		999	" %5d%c %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
		1000	i1,
		1001	ipt[ns - 1],
		1002	xw1,
		1003	yw1,
		1004	zw1,
		1005	xw2,
		1006	yw2,
		1007	zw2);
		1008
		1009	if ((ns == 2) \|\| (ns == 4))
		1010	isct = 1;
		1011
		1012	if (ns > 1)
		1013	{
		1014	readgm (gm, &ix, &iy, &x3, &y3, &z3, &x4, &y4, &z4, &dummy);
		1015
		1016	if ((ns == 2) \|\| (itg > 0))
		1017	{
		1018	x4 = xw1 + x3 - xw2;
		1019	y4 = yw1 + y3 - yw2;
		1020	z4 = zw1 + z3 - zw2;
		1021	}
		1022
		1023	fprintf (
		1024	output_fp,
		1025	"\n"
		1026	" %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
		1027	x3,
		1028	y3,
		1029	z3,
		1030	x4,
		1031	y4,
		1032	z4);
		1033
		1034	if (strcmp (gm, "SC") != 0)
		1035	{
		1036	fprintf (output_fp, "\n PATCH DATA ERROR");
		1037	stop (-1);
		1038	}
		1039
		1040	} /* if( ns > 1) */
		1041	else
		1042	{
		1043	xw2 = xw2 * TA;
		1044	yw2 = yw2 * TA;
		1045	}
		1046
		1047	patch (itg, ns, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
		1048
		1049	continue;
		1050
		1051	case 7: /* "sm" card, generate multiple-patch surface */
		1052
		1053	i1 = data.m + 1;
		1054	fprintf (
		1055	output_fp,
		1056	"\n"
		1057	" %5d%c %10.5f %11.5f %11.5f %11.5f %11.5f %11.5f"
		1058	" SURFACE - %d BY %d PATCHES",
		1059	i1,
		1060	ipt[1],
		1061	xw1,
		1062	yw1,
		1063	zw1,
		1064	xw2,
		1065	yw2,
		1066	zw2,
		1067	itg,
		1068	ns);
		1069
		1070	if ((itg < 1) \|\| (ns < 1))
		1071	{
		1072	fprintf (output_fp, "\n PATCH DATA ERROR");
		1073	stop (-1);
		1074	}
		1075
		1076	readgm (gm, &ix, &iy, &x3, &y3, &z3, &x4, &y4, &z4, &dummy);
		1077
		1078	if ((ns == 2) \|\| (itg > 0))
		1079	{
		1080	x4 = xw1 + x3 - xw2;
		1081	y4 = yw1 + y3 - yw2;
		1082	z4 = zw1 + z3 - zw2;
		1083	}
		1084
		1085	fprintf (
		1086	output_fp,
		1087	"\n"
		1088	" %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
		1089	x3,
		1090	y3,
		1091	z3,
		1092	x4,
		1093	y4,
		1094	z4);
		1095
		1096	if (strcmp (gm, "SC") != 0)
		1097	{
		1098	fprintf (output_fp, "\n PATCH DATA ERROR");
		1099	stop (-1);
		1100	}
		1101
		1102	patch (itg, ns, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
		1103
		1104	continue;
		1105
		1106	case 8: /* "ga" card, generate segment data for wire arc */
		1107
		1108	nwire++;
		1109	i1 = data.n + 1;
		1110	i2 = data.n + ns;
		1111
		1112	fprintf (
		1113	output_fp,
		1114	"\n"
		1115	" %5d ARC RADIUS: %9.5f FROM: %8.3f TO: %8.3f DEGREES"
		1116	" %11.5f %5d %5d %5d %4d",
		1117	nwire,
		1118	xw1,
		1119	yw1,
		1120	zw1,
		1121	xw2,
		1122	ns,
		1123	i1,
		1124	i2,
		1125	itg);
		1126
		1127	arc (itg, ns, xw1, yw1, zw1, xw2);
		1128
		1129	continue;
		1130
		1131	case 9: /* "sc" card */
		1132
		1133	if (isct == 0)
		1134	{
		1135	fprintf (output_fp, "\n PATCH DATA ERROR");
		1136	stop (-1);
		1137	}
		1138
		1139	i1 = data.m + 1;
		1140	ns++;
		1141
		1142	if ((itg != 0) \|\| ((ns != 2) && (ns != 4)))
		1143	{
		1144	fprintf (output_fp, "\n PATCH DATA ERROR");
		1145	stop (-1);
		1146	}
		1147
		1148	xs1 = x4;
		1149	ys1 = y4;
		1150	zs1 = z4;
		1151	xs2 = x3;
		1152	ys2 = y3;
		1153	zs2 = z3;
		1154	x3 = xw1;
		1155	y3 = yw1;
		1156	z3 = zw1;
		1157
		1158	if (ns == 4)
		1159	{
		1160	x4 = xw2;
		1161	y4 = yw2;
		1162	z4 = zw2;
		1163	}
		1164
		1165	xw1 = xs1;
		1166	yw1 = ys1;
		1167	zw1 = zs1;
		1168	xw2 = xs2;
		1169	yw2 = ys2;
		1170	zw2 = zs2;
		1171
		1172	if (ns != 4)
		1173	{
		1174	x4 = xw1 + x3 - xw2;
		1175	y4 = yw1 + y3 - yw2;
		1176	z4 = zw1 + z3 - zw2;
		1177	}
		1178
		1179	fprintf (
		1180	output_fp,
		1181	"\n"
		1182	" %5d%c %10.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
		1183	i1,
		1184	ipt[ns - 1],
		1185	xw1,
		1186	yw1,
		1187	zw1,
		1188	xw2,
		1189	yw2,
		1190	zw2);
		1191
		1192	fprintf (
		1193	output_fp,
		1194	"\n"
		1195	" %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
		1196	x3,
		1197	y3,
		1198	z3,
		1199	x4,
		1200	y4,
		1201	z4);
		1202
		1203	patch (itg, ns, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
		1204
		1205	continue;
		1206
		1207	case 10: /* "gh" card, generate helix */
		1208
		1209	nwire++;
		1210	i1 = data.n + 1;
		1211	i2 = data.n + ns;
		1212
		1213	fprintf (
		1214	output_fp,
		1215	"\n"
		1216	" %5d HELIX STRUCTURE - SPACING OF TURNS: %8.3f AXIAL"
		1217	" LENGTH: %8.3f %8.3f %5d %5d %5d %4d\n "
		1218	" RADIUS X1:%8.3f Y1:%8.3f X2:%8.3f Y2:%8.3f ",
		1219	nwire,
		1220	xw1,
		1221	yw1,
		1222	rad,
		1223	ns,
		1224	i1,
		1225	i2,
		1226	itg,
		1227	zw1,
		1228	xw2,
		1229	yw2,
		1230	zw2);
		1231
		1232	helix (xw1, yw1, zw1, xw2, yw2, zw2, rad, ns, itg);
		1233
		1234	continue;
		1235
		1236	case 11: /* "gf" card, not supported */
		1237	abort_on_error (-5);
		1238	continue;
		1239	default: /* error message */
		1240
		1241	fprintf (output_fp, "\n GEOMETRY DATA CARD ERROR");
		1242	fprintf (
		1243	output_fp,
		1244	"\n"
		1245	" %2s %3d %5d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
		1246	gm,
		1247	itg,
		1248	ns,
		1249	xw1,
		1250	yw1,
		1251	zw1,
		1252	xw2,
		1253	yw2,
		1254	zw2,
		1255	rad);
		1256
		1257	stop (-1);
		1258
		1259	} /* switch( gm_num ) */
		1260
		1261	} /* do */
		1262	while (TRUE);
		1263
		1264	return;
		1265	}
		1266
		1267	/-----------------------------------------------------------------------/
		1268
		1269	/* subroutine helix generates segment geometry */
		1270	/* data for a helix of ns segments */
		1271	void helix (
		1272	double s,
		1273	double hl,
		1274	double a1,
		1275	double b1,
		1276	double a2,
		1277	double b2,
		1278	double rad,
		1279	int ns,
		1280	int itg)
		1281	{
		1282	int ist, i, mreq;
		1283	// double turns;
		1284	double zinc, copy, sangle, hdia, turn, pitch, hmaj, hmin;
		1285
		1286	ist = data.n;
		1287	data.n += ns;
		1288	data.np = data.n;
		1289	data.mp = data.m;
		1290	data.ipsym = 0;
		1291
		1292	if (ns < 1)
		1293	return;
		1294
		1295	// turns= fabsl( hl/ s);
		1296	zinc = fabsl (hl / ns);
		1297
		1298	/* Reallocate tags buffer */
		1299	mem_realloc ((void ) &data.itag, (data.n + data.m) sizeof (int)); /????/
		1300
		1301	/* Reallocate wire buffers */
		1302	mreq = data.n * sizeof (double);
		1303	mem_realloc ((void *) &data.x1, mreq);
		1304	mem_realloc ((void *) &data.y1, mreq);
		1305	mem_realloc ((void *) &data.z1, mreq);
		1306	mem_realloc ((void *) &data.x2, mreq);
		1307	mem_realloc ((void *) &data.y2, mreq);
		1308	mem_realloc ((void *) &data.z2, mreq);
		1309	mem_realloc ((void *) &data.bi, mreq);
		1310
		1311	data.z1[ist] = 0.;
		1312	for (i = ist; i < data.n; i++)
		1313	{
		1314	data.bi[i] = rad;
		1315	data.itag[i] = itg;
		1316
		1317	if (i != ist)
		1318	data.z1[i] = data.z1[i - 1] + zinc;
		1319
		1320	data.z2[i] = data.z1[i] + zinc;
		1321
		1322	if (a2 == a1)
		1323	{
		1324	if (b1 == 0.)
		1325	b1 = a1;
		1326
		1327	data.x1[i] = a1 * cosl (2. * PI * data.z1[i] / s);
		1328	data.y1[i] = b1 * sinl (2. * PI * data.z1[i] / s);
		1329	data.x2[i] = a1 * cosl (2. * PI * data.z2[i] / s);
		1330	data.y2[i] = b1 * sinl (2. * PI * data.z2[i] / s);
		1331	}
		1332	else
		1333	{
		1334	if (b2 == 0.)
		1335	b2 = a2;
		1336
		1337	data.x1[i] = (a1 + (a2 - a1) * data.z1[i] / fabsl (hl)) *
		1338	cosl (2. * PI * data.z1[i] / s);
		1339	data.y1[i] = (b1 + (b2 - b1) * data.z1[i] / fabsl (hl)) *
		1340	sinl (2. * PI * data.z1[i] / s);
		1341	data.x2[i] = (a1 + (a2 - a1) * data.z2[i] / fabsl (hl)) *
		1342	cosl (2. * PI * data.z2[i] / s);
		1343	data.y2[i] = (b1 + (b2 - b1) * data.z2[i] / fabsl (hl)) *
		1344	sinl (2. * PI * data.z2[i] / s);
		1345
		1346	} /* if( a2 == a1) */
		1347
		1348	if (hl > 0.)
		1349	continue;
		1350
		1351	copy = data.x1[i];
		1352	data.x1[i] = data.y1[i];
		1353	data.y1[i] = copy;
		1354	copy = data.x2[i];
		1355	data.x2[i] = data.y2[i];
		1356	data.y2[i] = copy;
		1357
		1358	} /* for( i = ist; i < data.n; i++ ) */
		1359
		1360	if (a2 != a1)
		1361	{
		1362	sangle = atanl (a2 / (fabsl (hl) + (fabsl (hl) * a1) / (a2 - a1)));
		1363	fprintf (output_fp, "\n THE CONE ANGLE OF THE SPIRAL IS %10.4f", sangle);
		1364	return;
		1365	}
		1366
		1367	if (a1 == b1)
		1368	{
		1369	hdia = 2. * a1;
		1370	turn = hdia * PI;
		1371	pitch = atanl (s / (PI * hdia));
		1372	turn = turn / cosl (pitch);
		1373	pitch = 180. * pitch / PI;
		1374	}
		1375	else
		1376	{
		1377	if (a1 >= b1)
		1378	{
		1379	hmaj = 2. * a1;
		1380	hmin = 2. * b1;
		1381	}
		1382	else
		1383	{
		1384	hmaj = 2. * b1;
		1385	hmin = 2. * a1;
		1386	}
		1387
		1388	hdia = sqrtl ((hmaj * hmaj + hmin * hmin) / 2 * hmaj);
		1389	turn = 2. * PI * hdia;
		1390	pitch = (180. / PI) * atanl (s / (PI * hdia));
		1391
		1392	} /* if( a1 == b1) */
		1393
		1394	fprintf (
		1395	output_fp,
		1396	"\n"
		1397	" THE PITCH ANGLE IS: %.4f THE LENGTH OF WIRE/TURN IS: %.4f",
		1398	pitch,
		1399	turn);
		1400
		1401	return;
		1402	}
		1403
		1404	/-----------------------------------------------------------------------/
		1405
		1406	/* isegno returns the segment number of the mth segment having the */
		1407	/* tag number itagi. if itagi=0 segment number m is returned. */
		1408	int isegno (int itagi, int mx)
		1409	{
		1410	int icnt, i, iseg;
		1411
		1412	if (mx <= 0)
		1413	{
		1414	fprintf (
		1415	output_fp,
		1416	"\n CHECK DATA, PARAMETER SPECIFYING SEGMENT"
		1417	" POSITION IN A GROUP OF EQUAL TAGS MUST NOT BE ZERO");
		1418	stop (-1);
		1419	}
		1420
		1421	icnt = 0;
		1422	if (itagi == 0)
		1423	{
		1424	iseg = mx;
		1425	return (iseg);
		1426	}
		1427
		1428	if (data.n > 0)
		1429	{
		1430	for (i = 0; i < data.n; i++)
		1431	{
		1432	if (data.itag[i] != itagi)
		1433	continue;
		1434
		1435	icnt++;
		1436	if (icnt == mx)
		1437	{
		1438	iseg = i + 1;
		1439	return (iseg);
		1440	}
		1441
		1442	} /* for( i = 0; i < data.n; i++ ) */
		1443
		1444	} /* if( data.n > 0) */
		1445
		1446	fprintf (
		1447	output_fp,
		1448	"\n\n"
		1449	" NO SEGMENT HAS AN ITAG OF %d",
		1450	itagi);
		1451	stop (-1);
		1452
		1453	return (0);
		1454	}
		1455
		1456	/-----------------------------------------------------------------------/
		1457
		1458	/* subroutine move moves the structure with respect to its */
		1459	/* coordinate system or reproduces structure in new positions. */
		1460	/* structure is rotated about x,y,z axes by rox,roy,roz */
		1461	/* respectively, then shifted by xs,ys,zs */
		1462	void move (
		1463	double rox,
		1464	double roy,
		1465	double roz,
		1466	double xs,
		1467	double ys,
		1468	double zs,
		1469	int its,
		1470	int nrpt,
		1471	int itgi)
		1472	{
		1473	int nrp, ix, i1, k, ir, i, ii, mreq;
		1474	double sps, cps, sth, cth, sph, cph, xx, xy;
		1475	double xz, yx, yy, yz, zx, zy, zz, xi, yi, zi;
		1476
		1477	if (fabsl (rox) + fabsl (roy) > 1.0e-10)
		1478	data.ipsym = data.ipsym * 3;
		1479
		1480	sps = sinl (rox);
		1481	cps = cosl (rox);
		1482	sth = sinl (roy);
		1483	cth = cosl (roy);
		1484	sph = sinl (roz);
		1485	cph = cosl (roz);
		1486	xx = cph * cth;
		1487	xy = cph * sth * sps - sph * cps;
		1488	xz = cph * sth * cps + sph * sps;
		1489	yx = sph * cth;
		1490	yy = sph * sth * sps + cph * cps;
		1491	yz = sph * sth * cps - cph * sps;
		1492	zx = -sth;
		1493	zy = cth * sps;
		1494	zz = cth * cps;
		1495
		1496	if (nrpt == 0)
		1497	nrp = 1;
		1498	else
		1499	nrp = nrpt;
		1500
		1501	ix = 1;
		1502	if (data.n > 0)
		1503	{
		1504	i1 = isegno (its, 1);
		1505	if (i1 < 1)
		1506	i1 = 1;
		1507
		1508	ix = i1;
		1509	if (nrpt == 0)
		1510	k = i1 - 1;
		1511	else
		1512	{
		1513	k = data.n;
		1514	/* Reallocate tags buffer */
		1515	mreq = data.n + data.m + (data.n + 1 - i1) * nrpt;
		1516	mem_realloc ((void ) &data.itag, mreq sizeof (int));
		1517
		1518	/* Reallocate wire buffers */
		1519	mreq = (data.n + (data.n + 1 - i1) * nrpt) * sizeof (double);
		1520	mem_realloc ((void *) &data.x1, mreq);
		1521	mem_realloc ((void *) &data.y1, mreq);
		1522	mem_realloc ((void *) &data.z1, mreq);
		1523	mem_realloc ((void *) &data.x2, mreq);
		1524	mem_realloc ((void *) &data.y2, mreq);
		1525	mem_realloc ((void *) &data.z2, mreq);
		1526	mem_realloc ((void *) &data.bi, mreq);
		1527	}
		1528
		1529	for (ir = 0; ir < nrp; ir++)
		1530	{
		1531	for (i = i1 - 1; i < data.n; i++)
		1532	{
		1533	xi = data.x1[i];
		1534	yi = data.y1[i];
		1535	zi = data.z1[i];
		1536	data.x1[k] = xi * xx + yi * xy + zi * xz + xs;
		1537	data.y1[k] = xi * yx + yi * yy + zi * yz + ys;
		1538	data.z1[k] = xi * zx + yi * zy + zi * zz + zs;
		1539	xi = data.x2[i];
		1540	yi = data.y2[i];
		1541	zi = data.z2[i];
		1542	data.x2[k] = xi * xx + yi * xy + zi * xz + xs;
		1543	data.y2[k] = xi * yx + yi * yy + zi * yz + ys;
		1544	data.z2[k] = xi * zx + yi * zy + zi * zz + zs;
		1545	data.bi[k] = data.bi[i];
		1546	data.itag[k] = data.itag[i];
		1547	if (data.itag[i] != 0)
		1548	data.itag[k] = data.itag[i] + itgi;
		1549
		1550	k++;
		1551
		1552	} /* for( i = i1; i < data.n; i++ ) */
		1553
		1554	i1 = data.n + 1;
		1555	data.n = k;
		1556
		1557	} /* for( ir = 0; ir < nrp; ir++ ) */
		1558
		1559	} /* if( data.n >= n2) */
		1560
		1561	if (data.m > 0)
		1562	{
		1563	i1 = 0;
		1564	if (nrpt == 0)
		1565	k = 0;
		1566	else
		1567	k = data.m;
		1568
		1569	/* Reallocate patch buffers */
		1570	mreq = data.m * (1 + nrpt) * sizeof (double);
		1571	mem_realloc ((void *) &data.px, mreq);
		1572	mem_realloc ((void *) &data.py, mreq);
		1573	mem_realloc ((void *) &data.pz, mreq);
		1574	mem_realloc ((void *) &data.t1x, mreq);
		1575	mem_realloc ((void *) &data.t1y, mreq);
		1576	mem_realloc ((void *) &data.t1z, mreq);
		1577	mem_realloc ((void *) &data.t2x, mreq);
		1578	mem_realloc ((void *) &data.t2y, mreq);
		1579	mem_realloc ((void *) &data.t2z, mreq);
		1580	mem_realloc ((void *) &data.pbi, mreq);
		1581	mem_realloc ((void *) &data.psalp, mreq);
		1582
		1583	for (ii = 0; ii < nrp; ii++)
		1584	{
		1585	for (i = i1; i < data.m; i++)
		1586	{
		1587	xi = data.px[i];
		1588	yi = data.py[i];
		1589	zi = data.pz[i];
		1590	data.px[k] = xi * xx + yi * xy + zi * xz + xs;
		1591	data.py[k] = xi * yx + yi * yy + zi * yz + ys;
		1592	data.pz[k] = xi * zx + yi * zy + zi * zz + zs;
		1593	xi = data.t1x[i];
		1594	yi = data.t1y[i];
		1595	zi = data.t1z[i];
		1596	data.t1x[k] = xi * xx + yi * xy + zi * xz;
		1597	data.t1y[k] = xi * yx + yi * yy + zi * yz;
		1598	data.t1z[k] = xi * zx + yi * zy + zi * zz;
		1599	xi = data.t2x[i];
		1600	yi = data.t2y[i];
		1601	zi = data.t2z[i];
		1602	data.t2x[k] = xi * xx + yi * xy + zi * xz;
		1603	data.t2y[k] = xi * yx + yi * yy + zi * yz;
		1604	data.t2z[k] = xi * zx + yi * zy + zi * zz;
		1605	data.psalp[k] = data.psalp[i];
		1606	data.pbi[k] = data.pbi[i];
		1607	k++;
		1608
		1609	} /* for( i = i1; i < data.m; i++ ) */
		1610
		1611	i1 = data.m;
		1612	data.m = k;
		1613
		1614	} /* for( ii = 0; ii < nrp; ii++ ) */
		1615
		1616	} /* if( data.m >= m2) */
		1617
		1618	if ((nrpt == 0) && (ix == 1))
		1619	return;
		1620
		1621	data.np = data.n;
		1622	data.mp = data.m;
		1623	data.ipsym = 0;
		1624
		1625	return;
		1626	}
		1627
		1628	/-----------------------------------------------------------------------/
		1629
		1630	/* patch generates and modifies patch geometry data */
		1631	void patch (
		1632	int nx,
		1633	int ny,
		1634	double ax1,
		1635	double ay1,
		1636	double az1,
		1637	double ax2,
		1638	double ay2,
		1639	double az2,
		1640	double ax3,
		1641	double ay3,
		1642	double az3,
		1643	double ax4,
		1644	double ay4,
		1645	double az4)
		1646	{
		1647	int mi, ntp, iy, ix, mreq;
		1648	double s1x = 0., s1y = 0., s1z = 0., s2x = 0., s2y = 0., s2z = 0., xst = 0.;
		1649	double znv, xnv, ynv, xa, xn2, yn2, zn2, salpn, xs, ys, zs, xt, yt, zt;
		1650
		1651	/* new patches. for nx=0, ny=1,2,3,4 patch is (respectively) */
		1652	/* arbitrary, rectagular, triangular, or quadrilateral. */
		1653	/* for nx and ny > 0 a rectangular surface is produced with */
		1654	/* nx by ny rectangular patches. */
		1655
		1656	data.m++;
		1657	mi = data.m - 1;
		1658
		1659	/* Reallocate patch buffers */
		1660	mreq = data.m * sizeof (double);
		1661	mem_realloc ((void *) &data.px, mreq);
		1662	mem_realloc ((void *) &data.py, mreq);
		1663	mem_realloc ((void *) &data.pz, mreq);
		1664	mem_realloc ((void *) &data.t1x, mreq);
		1665	mem_realloc ((void *) &data.t1y, mreq);
		1666	mem_realloc ((void *) &data.t1z, mreq);
		1667	mem_realloc ((void *) &data.t2x, mreq);
		1668	mem_realloc ((void *) &data.t2y, mreq);
		1669	mem_realloc ((void *) &data.t2z, mreq);
		1670	mem_realloc ((void *) &data.pbi, mreq);
		1671	mem_realloc ((void *) &data.psalp, mreq);
		1672
		1673	if (nx > 0)
		1674	ntp = 2;
		1675	else
		1676	ntp = ny;
		1677
		1678	if (ntp <= 1)
		1679	{
		1680	data.px[mi] = ax1;
		1681	data.py[mi] = ay1;
		1682	data.pz[mi] = az1;
		1683	data.pbi[mi] = az2;
		1684	znv = cosl (ax2);
		1685	xnv = znv * cosl (ay2);
		1686	ynv = znv * sinl (ay2);
		1687	znv = sinl (ax2);
		1688	xa = sqrtl (xnv * xnv + ynv * ynv);
		1689
		1690	if (xa >= 1.0e-6)
		1691	{
		1692	data.t1x[mi] = -ynv / xa;
		1693	data.t1y[mi] = xnv / xa;
		1694	data.t1z[mi] = 0.;
		1695	}
		1696	else
		1697	{
		1698	data.t1x[mi] = 1.;
		1699	data.t1y[mi] = 0.;
		1700	data.t1z[mi] = 0.;
		1701	}
		1702
		1703	} /* if( ntp <= 1) */
		1704	else
		1705	{
		1706	s1x = ax2 - ax1;
		1707	s1y = ay2 - ay1;
		1708	s1z = az2 - az1;
		1709	s2x = ax3 - ax2;
		1710	s2y = ay3 - ay2;
		1711	s2z = az3 - az2;
		1712
		1713	if (nx != 0)
		1714	{
		1715	s1x = s1x / nx;
		1716	s1y = s1y / nx;
		1717	s1z = s1z / nx;
		1718	s2x = s2x / ny;
		1719	s2y = s2y / ny;
		1720	s2z = s2z / ny;
		1721	}
		1722
		1723	xnv = s1y * s2z - s1z * s2y;
		1724	ynv = s1z * s2x - s1x * s2z;
		1725	znv = s1x * s2y - s1y * s2x;
		1726	xa = sqrtl (xnv * xnv + ynv * ynv + znv * znv);
		1727	xnv = xnv / xa;
		1728	ynv = ynv / xa;
		1729	znv = znv / xa;
		1730	xst = sqrtl (s1x * s1x + s1y * s1y + s1z * s1z);
		1731	data.t1x[mi] = s1x / xst;
		1732	data.t1y[mi] = s1y / xst;
		1733	data.t1z[mi] = s1z / xst;
		1734
		1735	if (ntp <= 2)
		1736	{
		1737	data.px[mi] = ax1 + .5 * (s1x + s2x);
		1738	data.py[mi] = ay1 + .5 * (s1y + s2y);
		1739	data.pz[mi] = az1 + .5 * (s1z + s2z);
		1740	data.pbi[mi] = xa;
		1741	}
		1742	else
		1743	{
		1744	if (ntp != 4)
		1745	{
		1746	data.px[mi] = (ax1 + ax2 + ax3) / 3.;
		1747	data.py[mi] = (ay1 + ay2 + ay3) / 3.;
		1748	data.pz[mi] = (az1 + az2 + az3) / 3.;
		1749	data.pbi[mi] = .5 * xa;
		1750	}
		1751	else
		1752	{
		1753	s1x = ax3 - ax1;
		1754	s1y = ay3 - ay1;
		1755	s1z = az3 - az1;
		1756	s2x = ax4 - ax1;
		1757	s2y = ay4 - ay1;
		1758	s2z = az4 - az1;
		1759	xn2 = s1y * s2z - s1z * s2y;
		1760	yn2 = s1z * s2x - s1x * s2z;
		1761	zn2 = s1x * s2y - s1y * s2x;
		1762	xst = sqrtl (xn2 * xn2 + yn2 * yn2 + zn2 * zn2);
		1763	salpn = 1. / (3. * (xa + xst));
		1764	data.px[mi] =
		1765	(xa * (ax1 + ax2 + ax3) + xst * (ax1 + ax3 + ax4)) * salpn;
		1766	data.py[mi] =
		1767	(xa * (ay1 + ay2 + ay3) + xst * (ay1 + ay3 + ay4)) * salpn;
		1768	data.pz[mi] =
		1769	(xa * (az1 + az2 + az3) + xst * (az1 + az3 + az4)) * salpn;
		1770	data.pbi[mi] = .5 * (xa + xst);
		1771	s1x = (xnv * xn2 + ynv * yn2 + znv * zn2) / xst;
		1772
		1773	if (s1x <= 0.9998)
		1774	{
		1775	fprintf (
		1776	output_fp,
		1777	"\n ERROR -- CORNERS OF QUADRILATERAL"
		1778	" PATCH DO NOT LIE IN A PLANE");
		1779	stop (-1);
		1780	}
		1781
		1782	} /* if( ntp != 4) */
		1783
		1784	} /* if( ntp <= 2) */
		1785
		1786	} /* if( ntp <= 1) */
		1787
		1788	data.t2x[mi] = ynv * data.t1z[mi] - znv * data.t1y[mi];
		1789	data.t2y[mi] = znv * data.t1x[mi] - xnv * data.t1z[mi];
		1790	data.t2z[mi] = xnv * data.t1y[mi] - ynv * data.t1x[mi];
		1791	data.psalp[mi] = 1.;
		1792
		1793	if (nx != 0)
		1794	{
		1795	data.m += nx * ny - 1;
		1796
		1797	/* Reallocate patch buffers */
		1798	mreq = data.m * sizeof (double);
		1799	mem_realloc ((void *) &data.px, mreq);
		1800	mem_realloc ((void *) &data.py, mreq);
		1801	mem_realloc ((void *) &data.pz, mreq);
		1802	mem_realloc ((void *) &data.t1x, mreq);
		1803	mem_realloc ((void *) &data.t1y, mreq);
		1804	mem_realloc ((void *) &data.t1z, mreq);
		1805	mem_realloc ((void *) &data.t2x, mreq);
		1806	mem_realloc ((void *) &data.t2y, mreq);
		1807	mem_realloc ((void *) &data.t2z, mreq);
		1808	mem_realloc ((void *) &data.pbi, mreq);
		1809	mem_realloc ((void *) &data.psalp, mreq);
		1810
		1811	xn2 = data.px[mi] - s1x - s2x;
		1812	yn2 = data.py[mi] - s1y - s2y;
		1813	zn2 = data.pz[mi] - s1z - s2z;
		1814	xs = data.t1x[mi];
		1815	ys = data.t1y[mi];
		1816	zs = data.t1z[mi];
		1817	xt = data.t2x[mi];
		1818	yt = data.t2y[mi];
		1819	zt = data.t2z[mi];
		1820
		1821	for (iy = 0; iy < ny; iy++)
		1822	{
		1823	xn2 += s2x;
		1824	yn2 += s2y;
		1825	zn2 += s2z;
		1826
		1827	for (ix = 1; ix <= nx; ix++)
		1828	{
		1829	xst = (double) ix;
		1830	data.px[mi] = xn2 + xst * s1x;
		1831	data.py[mi] = yn2 + xst * s1y;
		1832	data.pz[mi] = zn2 + xst * s1z;
		1833	data.pbi[mi] = xa;
		1834	data.psalp[mi] = 1.;
		1835	data.t1x[mi] = xs;
		1836	data.t1y[mi] = ys;
		1837	data.t1z[mi] = zs;
		1838	data.t2x[mi] = xt;
		1839	data.t2y[mi] = yt;
		1840	data.t2z[mi] = zt;
		1841	mi++;
		1842	} /* for( ix = 0; ix < nx; ix++ ) */
		1843
		1844	} /* for( iy = 0; iy < ny; iy++ ) */
		1845
		1846	} /* if( nx != 0) */
		1847
		1848	data.ipsym = 0;
		1849	data.np = data.n;
		1850	data.mp = data.m;
		1851
		1852	return;
		1853	}
		1854
		1855	/-----------------------------------------------------------------------/
		1856
		1857	/* this function was an 'entry point' (part of) 'patch()' */
		1858	void subph (int nx, int ny)
		1859	{
		1860	int mia, ix, iy, mi, mreq;
		1861	double xs, ys, zs, xa, xst, s1x, s1y, s1z, s2x, s2y, s2z, saln, xt, yt;
		1862
		1863	/* Reallocate patch buffers */
		1864	if (ny == 0)
		1865	data.m += 3;
		1866	else
		1867	data.m += 4;
		1868
		1869	mreq = data.m * sizeof (double);
		1870	mem_realloc ((void *) &data.px, mreq);
		1871	mem_realloc ((void *) &data.py, mreq);
		1872	mem_realloc ((void *) &data.pz, mreq);
		1873	mem_realloc ((void *) &data.t1x, mreq);
		1874	mem_realloc ((void *) &data.t1y, mreq);
		1875	mem_realloc ((void *) &data.t1z, mreq);
		1876	mem_realloc ((void *) &data.t2x, mreq);
		1877	mem_realloc ((void *) &data.t2y, mreq);
		1878	mem_realloc ((void *) &data.t2z, mreq);
		1879	mem_realloc ((void *) &data.pbi, mreq);
		1880	mem_realloc ((void *) &data.psalp, mreq);
		1881	mem_realloc ((void ) &data.icon1, (data.n + data.m) sizeof (int));
		1882	mem_realloc ((void ) &data.icon2, (data.n + data.m) sizeof (int));
		1883
		1884	/* Shift patches to make room for new ones */
		1885	if ((ny == 0) && (nx != data.m))
		1886	{
		1887	for (iy = data.m - 1; iy > nx + 2; iy--)
		1888	{
		1889	ix = iy - 3;
		1890	data.px[iy] = data.px[ix];
		1891	data.py[iy] = data.py[ix];
		1892	data.pz[iy] = data.pz[ix];
		1893	data.pbi[iy] = data.pbi[ix];
		1894	data.psalp[iy] = data.psalp[ix];
		1895	data.t1x[iy] = data.t1x[ix];
		1896	data.t1y[iy] = data.t1y[ix];
		1897	data.t1z[iy] = data.t1z[ix];
		1898	data.t2x[iy] = data.t2x[ix];
		1899	data.t2y[iy] = data.t2y[ix];
		1900	data.t2z[iy] = data.t2z[ix];
		1901	}
		1902
		1903	} /* if( (ny == 0) \|\| (nx != m) ) */
		1904
		1905	/* divide patch for connection */
		1906	mi = nx - 1;
		1907	xs = data.px[mi];
		1908	ys = data.py[mi];
		1909	zs = data.pz[mi];
		1910	xa = data.pbi[mi] / 4.;
		1911	xst = sqrtl (xa) / 2.;
		1912	s1x = data.t1x[mi];
		1913	s1y = data.t1y[mi];
		1914	s1z = data.t1z[mi];
		1915	s2x = data.t2x[mi];
		1916	s2y = data.t2y[mi];
		1917	s2z = data.t2z[mi];
		1918	saln = data.psalp[mi];
		1919	xt = xst;
		1920	yt = xst;
		1921
		1922	if (ny == 0)
		1923	mia = mi;
		1924	else
		1925	{
		1926	data.mp++;
		1927	mia = data.m - 1;
		1928	}
		1929
		1930	for (ix = 1; ix <= 4; ix++)
		1931	{
		1932	data.px[mia] = xs + xt * s1x + yt * s2x;
		1933	data.py[mia] = ys + xt * s1y + yt * s2y;
		1934	data.pz[mia] = zs + xt * s1z + yt * s2z;
		1935	data.pbi[mia] = xa;
		1936	data.t1x[mia] = s1x;
		1937	data.t1y[mia] = s1y;
		1938	data.t1z[mia] = s1z;
		1939	data.t2x[mia] = s2x;
		1940	data.t2y[mia] = s2y;
		1941	data.t2z[mia] = s2z;
		1942	data.psalp[mia] = saln;
		1943
		1944	if (ix == 2)
		1945	yt = -yt;
		1946
		1947	if ((ix == 1) \|\| (ix == 3))
		1948	xt = -xt;
		1949
		1950	mia++;
		1951	}
		1952
		1953	if (nx <= data.mp)
		1954	data.mp += 3;
		1955
		1956	if (ny > 0)
		1957	data.pz[mi] = 10000.;
		1958
		1959	return;
		1960	}
		1961
		1962	/-----------------------------------------------------------------------/
		1963
		1964	void readgm (
		1965	char *gm,
		1966	int *i1,
		1967	int *i2,
		1968	double *x1,
		1969	double *y1,
		1970	double *z1,
		1971	double *x2,
		1972	double *y2,
		1973	double *z2,
		1974	double *rad)
		1975	{
		1976	char line_buf[134];
		1977	int nlin, i, line_idx;
		1978	int nint = 2, nflt = 7;
		1979	int iarr[2] = {0, 0};
		1980	double rarr[7] = {0., 0., 0., 0., 0., 0., 0.};
		1981
		1982	/* read a line from input file */
		1983	load_line (line_buf, input_fp);
		1984	{
		1985	gm[0] = 'E';
		1986	gm[1] = 'N';
		1987	gm[2] = '\0';
		1988	nlin = 2;
		1989	}
		1990
		1991	/* get line length */
		1992	nlin = strlen (line_buf);
		1993
		1994	/* abort if card's mnemonic too short or missing */
		1995	if (nlin < 2)
		1996	{
		1997	fprintf (
		1998	output_fp,
		1999	"\n GEOMETRY DATA CARD ERROR:"
		2000	"\n CARD'S MNEMONIC CODE TOO SHORT OR MISSING.");
		2001	stop (-1);
		2002	}
		2003
		2004	/* extract card's mnemonic code */
		2005	/* was strncpy (gm, line_buf, 2) ; */
		2006	gm[0] = line_buf[0];
		2007	gm[1] = line_buf[1];
		2008	gm[2] = '\0';
		2009
		2010	/* Exit if "XT" command read (for testing) */
		2011	if (strcmp (gm, "XT") == 0)
		2012	{
		2013	fprintf (stderr, "\nnec2c: Exiting after an \"XT\" command in readgm()\n");
		2014	fprintf (
		2015	output_fp, "\n\n nec2c: Exiting after an \"XT\" command in readgm()");
		2016	stop (0);
		2017	}
		2018
		2019	/* Return if only mnemonic on card */
		2020	if (nlin == 2)
		2021	{
		2022	i1 = i2 = 0;
		2023	x1 = y1 = z1 = x2 = y2 = z2 = *rad = 0.;
		2024	return;
		2025	}
		2026
		2027	/* read integers from line */
		2028	line_idx = 1;
		2029	for (i = 0; i < nint; i++)
		2030	{
		2031	/* Find first numerical character */
		2032	while (((line_buf[++line_idx] < '0') \|\| (line_buf[line_idx] > '9')) &&
		2033	(line_buf[line_idx] != '+') && (line_buf[line_idx] != '-'))
		2034	if ((line_buf[line_idx] == '\0'))
		2035	{
		2036	*i1 = iarr[0];
		2037	*i2 = iarr[1];
		2038	*x1 = rarr[0];
		2039	*y1 = rarr[1];
		2040	*z1 = rarr[2];
		2041	*x2 = rarr[3];
		2042	*y2 = rarr[4];
		2043	*z2 = rarr[5];
		2044	*rad = rarr[6];
		2045	return;
		2046	}
		2047
		2048	/* read an integer from line */
		2049	iarr[i] = atoi (&line_buf[line_idx]);
		2050
		2051	/* traverse numerical field to next ' ' or ',' or '\0' */
		2052	line_idx--;
		2053	while ((line_buf[++line_idx] != ' ') && (line_buf[line_idx] != ' ') &&
		2054	(line_buf[line_idx] != ',') && (line_buf[line_idx] != '\0'))
		2055	{
		2056	/* test for non-numerical characters */
		2057	if (((line_buf[line_idx] < '0') \|\| (line_buf[line_idx] > '9')) &&
		2058	(line_buf[line_idx] != '+') && (line_buf[line_idx] != '-'))
		2059	{
		2060	fprintf (
		2061	output_fp,
		2062	"\n GEOMETRY DATA CARD \"%s\" ERROR:"
		2063	"\n NON-NUMERICAL CHARACTER '%c' IN INTEGER FIELD AT "
		2064	"CHAR. %d\n",
		2065	gm,
		2066	line_buf[line_idx],
		2067	(line_idx + 1));
		2068	stop (-1);
		2069	}
		2070
		2071	} /* while( (line_buff[++line_idx] ... */
		2072
		2073	/* Return on end of line */
		2074	if (line_buf[line_idx] == '\0')
		2075	{
		2076	*i1 = iarr[0];
		2077	*i2 = iarr[1];
		2078	*x1 = rarr[0];
		2079	*y1 = rarr[1];
		2080	*z1 = rarr[2];
		2081	*x2 = rarr[3];
		2082	*y2 = rarr[4];
		2083	*z2 = rarr[5];
		2084	*rad = rarr[6];
		2085	return;
		2086	}
		2087
		2088	} /* for( i = 0; i < nint; i++ ) */
		2089
		2090	/* read doubles from line */
		2091	for (i = 0; i < nflt; i++)
		2092	{
		2093	/* Find first numerical character */
		2094	while (((line_buf[++line_idx] < '0') \|\| (line_buf[line_idx] > '9')) &&
		2095	(line_buf[line_idx] != '+') && (line_buf[line_idx] != '-') &&
		2096	(line_buf[line_idx] != '.'))
		2097	if ((line_buf[line_idx] == '\0'))
		2098	{
		2099	*i1 = iarr[0];
		2100	*i2 = iarr[1];
		2101	*x1 = rarr[0];
		2102	*y1 = rarr[1];
		2103	*z1 = rarr[2];
		2104	*x2 = rarr[3];
		2105	*y2 = rarr[4];
		2106	*z2 = rarr[5];
		2107	*rad = rarr[6];
		2108	return;
		2109	}
		2110
		2111	/* read a double from line */
		2112	rarr[i] = atof (&line_buf[line_idx]);
		2113
		2114	/* traverse numerical field to next ' ' or ',' or '\0' */
		2115	line_idx--;
		2116	while ((line_buf[++line_idx] != ' ') && (line_buf[line_idx] != ' ') &&
		2117	(line_buf[line_idx] != ',') && (line_buf[line_idx] != '\0'))
		2118	{
		2119	/* test for non-numerical characters */
		2120	if (((line_buf[line_idx] < '0') \|\| (line_buf[line_idx] > '9')) &&
		2121	(line_buf[line_idx] != '.') && (line_buf[line_idx] != '+') &&
		2122	(line_buf[line_idx] != '-') && (line_buf[line_idx] != 'E') &&
		2123	(line_buf[line_idx] != 'e'))
		2124	{
		2125	fprintf (
		2126	output_fp,
		2127	"\n GEOMETRY DATA CARD \"%s\" ERROR:"
		2128	"\n NON-NUMERICAL CHARACTER '%c' IN FLOAT FIELD AT CHAR. "
		2129	"%d.\n",
		2130	gm,
		2131	line_buf[line_idx],
		2132	(line_idx + 1));
		2133	stop (-1);
		2134	}
		2135
		2136	} /* while( (line_buff[++line_idx] ... */
		2137
		2138	/* Return on end of line */
		2139	if (line_buf[line_idx] == '\0')
		2140	{
		2141	*i1 = iarr[0];
		2142	*i2 = iarr[1];
		2143	*x1 = rarr[0];
		2144	*y1 = rarr[1];
		2145	*z1 = rarr[2];
		2146	*x2 = rarr[3];
		2147	*y2 = rarr[4];
		2148	*z2 = rarr[5];
		2149	*rad = rarr[6];
		2150	return;
		2151	}
		2152
		2153	} /* for( i = 0; i < nflt; i++ ) */
		2154
		2155	*i1 = iarr[0];
		2156	*i2 = iarr[1];
		2157	*x1 = rarr[0];
		2158	*y1 = rarr[1];
		2159	*z1 = rarr[2];
		2160	*x2 = rarr[3];
		2161	*y2 = rarr[4];
		2162	*z2 = rarr[5];
		2163	*rad = rarr[6];
		2164
		2165	return;
		2166	}
		2167
		2168	/-----------------------------------------------------------------------/
		2169
		2170	/* reflc reflects partial structure along x,y, or z axes or rotates */
		2171	/* structure to complete a symmetric structure. */
		2172	void reflc (int ix, int iy, int iz, int itx, int nop)
		2173	{
		2174	int iti, i, nx, itagi, k, mreq;
		2175	double e1, e2, fnop, sam, cs, ss, xk, yk;
		2176
		2177	data.np = data.n;
		2178	data.mp = data.m;
		2179	data.ipsym = 0;
		2180	iti = itx;
		2181
		2182	if (ix >= 0)
		2183	{
		2184	if (nop == 0)
		2185	return;
		2186
		2187	data.ipsym = 1;
		2188
		2189	/* reflect along z axis */
		2190	if (iz != 0)
		2191	{
		2192	data.ipsym = 2;
		2193
		2194	if (data.n > 0)
		2195	{
		2196	/* Reallocate tags buffer */
		2197	mem_realloc (
		2198	(void ) &data.itag, (2 data.n + data.m) * sizeof (int));
		2199
		2200	/* Reallocate wire buffers */
		2201	mreq = 2 * data.n * sizeof (double);
		2202	mem_realloc ((void *) &data.x1, mreq);
		2203	mem_realloc ((void *) &data.y1, mreq);
		2204	mem_realloc ((void *) &data.z1, mreq);
		2205	mem_realloc ((void *) &data.x2, mreq);
		2206	mem_realloc ((void *) &data.y2, mreq);
		2207	mem_realloc ((void *) &data.z2, mreq);
		2208	mem_realloc ((void *) &data.bi, mreq);
		2209
		2210	for (i = 0; i < data.n; i++)
		2211	{
		2212	nx = i + data.n;
		2213	e1 = data.z1[i];
		2214	e2 = data.z2[i];
		2215
		2216	if ((fabsl (e1) + fabsl (e2) <= 1.0e-5) \|\|
		2217	(e1 * e2 < -1.0e-6))
		2218	{
		2219	fprintf (
		2220	output_fp,
		2221	"\n GEOMETRY DATA ERROR--SEGMENT %d"
		2222	" LIES IN PLANE OF SYMMETRY",
		2223	i + 1);
		2224	stop (-1);
		2225	}
		2226
		2227	data.x1[nx] = data.x1[i];
		2228	data.y1[nx] = data.y1[i];
		2229	data.z1[nx] = -e1;
		2230	data.x2[nx] = data.x2[i];
		2231	data.y2[nx] = data.y2[i];
		2232	data.z2[nx] = -e2;
		2233	itagi = data.itag[i];
		2234
		2235	if (itagi == 0)
		2236	data.itag[nx] = 0;
		2237	if (itagi != 0)
		2238	data.itag[nx] = itagi + iti;
		2239
		2240	data.bi[nx] = data.bi[i];
		2241
		2242	} /* for( i = 0; i < data.n; i++ ) */
		2243
		2244	data.n = data.n * 2;
		2245	iti = iti * 2;
		2246
		2247	} /* if( data.n > 0) */
		2248
		2249	if (data.m > 0)
		2250	{
		2251	/* Reallocate patch buffers */
		2252	mreq = 2 * data.m * sizeof (double);
		2253	mem_realloc ((void *) &data.px, mreq);
		2254	mem_realloc ((void *) &data.py, mreq);
		2255	mem_realloc ((void *) &data.pz, mreq);
		2256	mem_realloc ((void *) &data.t1x, mreq);
		2257	mem_realloc ((void *) &data.t1y, mreq);
		2258	mem_realloc ((void *) &data.t1z, mreq);
		2259	mem_realloc ((void *) &data.t2x, mreq);
		2260	mem_realloc ((void *) &data.t2y, mreq);
		2261	mem_realloc ((void *) &data.t2z, mreq);
		2262	mem_realloc ((void *) &data.pbi, mreq);
		2263	mem_realloc ((void *) &data.psalp, mreq);
		2264
		2265	for (i = 0; i < data.m; i++)
		2266	{
		2267	nx = i + data.m;
		2268	if (fabsl (data.pz[i]) <= 1.0e-10)
		2269	{
		2270	fprintf (
		2271	output_fp,
		2272	"\n GEOMETRY DATA ERROR--PATCH %d"
		2273	" LIES IN PLANE OF SYMMETRY",
		2274	i + 1);
		2275	stop (-1);
		2276	}
		2277
		2278	data.px[nx] = data.px[i];
		2279	data.py[nx] = data.py[i];
		2280	data.pz[nx] = -data.pz[i];
		2281	data.t1x[nx] = data.t1x[i];
		2282	data.t1y[nx] = data.t1y[i];
		2283	data.t1z[nx] = -data.t1z[i];
		2284	data.t2x[nx] = data.t2x[i];
		2285	data.t2y[nx] = data.t2y[i];
		2286	data.t2z[nx] = -data.t2z[i];
		2287	data.psalp[nx] = -data.psalp[i];
		2288	data.pbi[nx] = data.pbi[i];
		2289	}
		2290
		2291	data.m = data.m * 2;
		2292
		2293	} /* if( data.m >= m2) */
		2294
		2295	} /* if( iz != 0) */
		2296
		2297	/* reflect along y axis */
		2298	if (iy != 0)
		2299	{
		2300	if (data.n > 0)
		2301	{
		2302	/* Reallocate tags buffer */
		2303	mem_realloc (
		2304	(void *) &data.itag,
		2305	(2 * data.n + data.m) * sizeof (int)); /????/
		2306
		2307	/* Reallocate wire buffers */
		2308	mreq = 2 * data.n * sizeof (double);
		2309	mem_realloc ((void *) &data.x1, mreq);
		2310	mem_realloc ((void *) &data.y1, mreq);
		2311	mem_realloc ((void *) &data.z1, mreq);
		2312	mem_realloc ((void *) &data.x2, mreq);
		2313	mem_realloc ((void *) &data.y2, mreq);
		2314	mem_realloc ((void *) &data.z2, mreq);
		2315	mem_realloc ((void *) &data.bi, mreq);
		2316
		2317	for (i = 0; i < data.n; i++)
		2318	{
		2319	nx = i + data.n;
		2320	e1 = data.y1[i];
		2321	e2 = data.y2[i];
		2322
		2323	if ((fabsl (e1) + fabsl (e2) <= 1.0e-5) \|\|
		2324	(e1 * e2 < -1.0e-6))
		2325	{
		2326	fprintf (
		2327	output_fp,
		2328	"\n GEOMETRY DATA ERROR--SEGMENT %d"
		2329	" LIES IN PLANE OF SYMMETRY",
		2330	i + 1);
		2331	stop (-1);
		2332	}
		2333
		2334	data.x1[nx] = data.x1[i];
		2335	data.y1[nx] = -e1;
		2336	data.z1[nx] = data.z1[i];
		2337	data.x2[nx] = data.x2[i];
		2338	data.y2[nx] = -e2;
		2339	data.z2[nx] = data.z2[i];
		2340	itagi = data.itag[i];
		2341
		2342	if (itagi == 0)
		2343	data.itag[nx] = 0;
		2344	if (itagi != 0)
		2345	data.itag[nx] = itagi + iti;
		2346
		2347	data.bi[nx] = data.bi[i];
		2348
		2349	} /* for( i = n2-1; i < data.n; i++ ) */
		2350
		2351	data.n = data.n * 2;
		2352	iti = iti * 2;
		2353
		2354	} /* if( data.n >= n2) */
		2355
		2356	if (data.m > 0)
		2357	{
		2358	/* Reallocate patch buffers */
		2359	mreq = 2 * data.m * sizeof (double);
		2360	mem_realloc ((void *) &data.px, mreq);
		2361	mem_realloc ((void *) &data.py, mreq);
		2362	mem_realloc ((void *) &data.pz, mreq);
		2363	mem_realloc ((void *) &data.t1x, mreq);
		2364	mem_realloc ((void *) &data.t1y, mreq);
		2365	mem_realloc ((void *) &data.t1z, mreq);
		2366	mem_realloc ((void *) &data.t2x, mreq);
		2367	mem_realloc ((void *) &data.t2y, mreq);
		2368	mem_realloc ((void *) &data.t2z, mreq);
		2369	mem_realloc ((void *) &data.pbi, mreq);
		2370	mem_realloc ((void *) &data.psalp, mreq);
		2371
		2372	for (i = 0; i < data.m; i++)
		2373	{
		2374	nx = i + data.m;
		2375	if (fabsl (data.py[i]) <= 1.0e-10)
		2376	{
		2377	fprintf (
		2378	output_fp,
		2379	"\n GEOMETRY DATA ERROR--PATCH %d"
		2380	" LIES IN PLANE OF SYMMETRY",
		2381	i + 1);
		2382	stop (-1);
		2383	}
		2384
		2385	data.px[nx] = data.px[i];
		2386	data.py[nx] = -data.py[i];
		2387	data.pz[nx] = data.pz[i];
		2388	data.t1x[nx] = data.t1x[i];
		2389	data.t1y[nx] = -data.t1y[i];
		2390	data.t1z[nx] = data.t1z[i];
		2391	data.t2x[nx] = data.t2x[i];
		2392	data.t2y[nx] = -data.t2y[i];
		2393	data.t2z[nx] = data.t2z[i];
		2394	data.psalp[nx] = -data.psalp[i];
		2395	data.pbi[nx] = data.pbi[i];
		2396
		2397	} /* for( i = m2; i <= data.m; i++ ) */
		2398
		2399	data.m = data.m * 2;
		2400
		2401	} /* if( data.m >= m2) */
		2402
		2403	} /* if( iy != 0) */
		2404
		2405	/* reflect along x axis */
		2406	if (ix == 0)
		2407	return;
		2408
		2409	if (data.n > 0)
		2410	{
		2411	/* Reallocate tags buffer */
		2412	mem_realloc (
		2413	(void *) &data.itag,
		2414	(2 * data.n + data.m) * sizeof (int)); /????/
		2415
		2416	/* Reallocate wire buffers */
		2417	mreq = 2 * data.n * sizeof (double);
		2418	mem_realloc ((void *) &data.x1, mreq);
		2419	mem_realloc ((void *) &data.y1, mreq);
		2420	mem_realloc ((void *) &data.z1, mreq);
		2421	mem_realloc ((void *) &data.x2, mreq);
		2422	mem_realloc ((void *) &data.y2, mreq);
		2423	mem_realloc ((void *) &data.z2, mreq);
		2424	mem_realloc ((void *) &data.bi, mreq);
		2425
		2426	for (i = 0; i < data.n; i++)
		2427	{
		2428	nx = i + data.n;
		2429	e1 = data.x1[i];
		2430	e2 = data.x2[i];
		2431
		2432	if ((fabsl (e1) + fabsl (e2) <= 1.0e-5) \|\| (e1 * e2 < -1.0e-6))
		2433	{
		2434	fprintf (
		2435	output_fp,
		2436	"\n GEOMETRY DATA ERROR--SEGMENT %d"
		2437	" LIES IN PLANE OF SYMMETRY",
		2438	i + 1);
		2439	stop (-1);
		2440	}
		2441
		2442	data.x1[nx] = -e1;
		2443	data.y1[nx] = data.y1[i];
		2444	data.z1[nx] = data.z1[i];
		2445	data.x2[nx] = -e2;
		2446	data.y2[nx] = data.y2[i];
		2447	data.z2[nx] = data.z2[i];
		2448	itagi = data.itag[i];
		2449
		2450	if (itagi == 0)
		2451	data.itag[nx] = 0;
		2452	if (itagi != 0)
		2453	data.itag[nx] = itagi + iti;
		2454
		2455	data.bi[nx] = data.bi[i];
		2456	}
		2457
		2458	data.n = data.n * 2;
		2459
		2460	} /* if( data.n > 0) */
		2461
		2462	if (data.m == 0)
		2463	return;
		2464
		2465	/* Reallocate patch buffers */
		2466	mreq = 2 * data.m * sizeof (double);
		2467	mem_realloc ((void *) &data.px, mreq);
		2468	mem_realloc ((void *) &data.py, mreq);
		2469	mem_realloc ((void *) &data.pz, mreq);
		2470	mem_realloc ((void *) &data.t1x, mreq);
		2471	mem_realloc ((void *) &data.t1y, mreq);
		2472	mem_realloc ((void *) &data.t1z, mreq);
		2473	mem_realloc ((void *) &data.t2x, mreq);
		2474	mem_realloc ((void *) &data.t2y, mreq);
		2475	mem_realloc ((void *) &data.t2z, mreq);
		2476	mem_realloc ((void *) &data.pbi, mreq);
		2477	mem_realloc ((void *) &data.psalp, mreq);
		2478
		2479	for (i = 0; i < data.m; i++)
		2480	{
		2481	nx = i + data.m;
		2482	if (fabsl (data.px[i]) <= 1.0e-10)
		2483	{
		2484	fprintf (
		2485	output_fp,
		2486	"\n GEOMETRY DATA ERROR--PATCH %d"
		2487	" LIES IN PLANE OF SYMMETRY",
		2488	i + 1);
		2489	stop (-1);
		2490	}
		2491
		2492	data.px[nx] = -data.px[i];
		2493	data.py[nx] = data.py[i];
		2494	data.pz[nx] = data.pz[i];
		2495	data.t1x[nx] = -data.t1x[i];
		2496	data.t1y[nx] = data.t1y[i];
		2497	data.t1z[nx] = data.t1z[i];
		2498	data.t2x[nx] = -data.t2x[i];
		2499	data.t2y[nx] = data.t2y[i];
		2500	data.t2z[nx] = data.t2z[i];
		2501	data.psalp[nx] = -data.psalp[i];
		2502	data.pbi[nx] = data.pbi[i];
		2503	}
		2504
		2505	data.m = data.m * 2;
		2506	return;
		2507
		2508	} /* if( ix >= 0) */
		2509
		2510	/* reproduce structure with rotation to form cylindrical structure */
		2511	fnop = (double) nop;
		2512	data.ipsym = -1;
		2513	sam = TP / fnop;
		2514	cs = cosl (sam);
		2515	ss = sinl (sam);
		2516
		2517	if (data.n > 0)
		2518	{
		2519	data.n *= nop;
		2520	nx = data.np;
		2521
		2522	/* Reallocate tags buffer */
		2523	mem_realloc ((void ) &data.itag, (data.n + data.m) sizeof (int)); /????/
		2524
		2525	/* Reallocate wire buffers */
		2526	mreq = data.n * sizeof (double);
		2527	mem_realloc ((void *) &data.x1, mreq);
		2528	mem_realloc ((void *) &data.y1, mreq);
		2529	mem_realloc ((void *) &data.z1, mreq);
		2530	mem_realloc ((void *) &data.x2, mreq);
		2531	mem_realloc ((void *) &data.y2, mreq);
		2532	mem_realloc ((void *) &data.z2, mreq);
		2533	mem_realloc ((void *) &data.bi, mreq);
		2534
		2535	for (i = nx; i < data.n; i++)
		2536	{
		2537	k = i - data.np;
		2538	xk = data.x1[k];
		2539	yk = data.y1[k];
		2540	data.x1[i] = xk * cs - yk * ss;
		2541	data.y1[i] = xk * ss + yk * cs;
		2542	data.z1[i] = data.z1[k];
		2543	xk = data.x2[k];
		2544	yk = data.y2[k];
		2545	data.x2[i] = xk * cs - yk * ss;
		2546	data.y2[i] = xk * ss + yk * cs;
		2547	data.z2[i] = data.z2[k];
		2548	data.bi[i] = data.bi[k];
		2549	itagi = data.itag[k];
		2550
		2551	if (itagi == 0)
		2552	data.itag[i] = 0;
		2553	if (itagi != 0)
		2554	data.itag[i] = itagi + iti;
		2555	}
		2556
		2557	} /* if( data.n >= n2) */
		2558
		2559	if (data.m == 0)
		2560	return;
		2561
		2562	data.m *= nop;
		2563	nx = data.mp;
		2564
		2565	/* Reallocate patch buffers */
		2566	mreq = data.m * sizeof (double);
		2567	mem_realloc ((void *) &data.px, mreq);
		2568	mem_realloc ((void *) &data.py, mreq);
		2569	mem_realloc ((void *) &data.pz, mreq);
		2570	mem_realloc ((void *) &data.t1x, mreq);
		2571	mem_realloc ((void *) &data.t1y, mreq);
		2572	mem_realloc ((void *) &data.t1z, mreq);
		2573	mem_realloc ((void *) &data.t2x, mreq);
		2574	mem_realloc ((void *) &data.t2y, mreq);
		2575	mem_realloc ((void *) &data.t2z, mreq);
		2576	mem_realloc ((void *) &data.pbi, mreq);
		2577	mem_realloc ((void *) &data.psalp, mreq);
		2578
		2579	for (i = nx; i < data.m; i++)
		2580	{
		2581	k = i - data.mp;
		2582	xk = data.px[k];
		2583	yk = data.py[k];
		2584	data.px[i] = xk * cs - yk * ss;
		2585	data.py[i] = xk * ss + yk * cs;
		2586	data.pz[i] = data.pz[k];
		2587	xk = data.t1x[k];
		2588	yk = data.t1y[k];
		2589	data.t1x[i] = xk * cs - yk * ss;
		2590	data.t1y[i] = xk * ss + yk * cs;
		2591	data.t1z[i] = data.t1z[k];
		2592	xk = data.t2x[k];
		2593	yk = data.t2y[k];
		2594	data.t2x[i] = xk * cs - yk * ss;
		2595	data.t2y[i] = xk * ss + yk * cs;
		2596	data.t2z[i] = data.t2z[k];
		2597	data.psalp[i] = data.psalp[k];
		2598	data.pbi[i] = data.pbi[k];
		2599
		2600	} /* for( i = nx; i < data.m; i++ ) */
		2601
		2602	return;
		2603	}
		2604
		2605	/-----------------------------------------------------------------------/
		2606
		2607	/* subroutine wire generates segment geometry */
		2608	/* data for a straight wire of ns segments. */
		2609	void wire (
		2610	double xw1,
		2611	double yw1,
		2612	double zw1,
		2613	double xw2,
		2614	double yw2,
		2615	double zw2,
		2616	double rad,
		2617	double rdel,
		2618	double rrad,
		2619	int ns,
		2620	int itg)
		2621	{
		2622	int ist, i, mreq;
		2623	double xd, yd, zd, delz, rd, fns, radz;
		2624	double xs1, ys1, zs1, xs2, ys2, zs2;
		2625
		2626	ist = data.n;
		2627	data.n = data.n + ns;
		2628	data.np = data.n;
		2629	data.mp = data.m;
		2630	data.ipsym = 0;
		2631
		2632	if (ns < 1)
		2633	return;
		2634
		2635	/* Reallocate tags buffer */
		2636	mem_realloc ((void ) &data.itag, (data.n + data.m) sizeof (int)); /????/
		2637
		2638	/* Reallocate wire buffers */
		2639	mreq = data.n * sizeof (double);
		2640	mem_realloc ((void *) &data.x1, mreq);
		2641	mem_realloc ((void *) &data.y1, mreq);
		2642	mem_realloc ((void *) &data.z1, mreq);
		2643	mem_realloc ((void *) &data.x2, mreq);
		2644	mem_realloc ((void *) &data.y2, mreq);
		2645	mem_realloc ((void *) &data.z2, mreq);
		2646	mem_realloc ((void *) &data.bi, mreq);
		2647
		2648	xd = xw2 - xw1;
		2649	yd = yw2 - yw1;
		2650	zd = zw2 - zw1;
		2651
		2652	if (fabsl (rdel - 1.) >= 1.0e-6)
		2653	{
		2654	delz = sqrtl (xd * xd + yd * yd + zd * zd);
		2655	xd = xd / delz;
		2656	yd = yd / delz;
		2657	zd = zd / delz;
		2658	delz = delz * (1. - rdel) / (1. - powl (rdel, ns));
		2659	rd = rdel;
		2660	}
		2661	else
		2662	{
		2663	fns = ns;
		2664	xd = xd / fns;
		2665	yd = yd / fns;
		2666	zd = zd / fns;
		2667	delz = 1.;
		2668	rd = 1.;
		2669	}
		2670
		2671	radz = rad;
		2672	xs1 = xw1;
		2673	ys1 = yw1;
		2674	zs1 = zw1;
		2675
		2676	for (i = ist; i < data.n; i++)
		2677	{
		2678	data.itag[i] = itg;
		2679	xs2 = xs1 + xd * delz;
		2680	ys2 = ys1 + yd * delz;
		2681	zs2 = zs1 + zd * delz;
		2682	data.x1[i] = xs1;
		2683	data.y1[i] = ys1;
		2684	data.z1[i] = zs1;
		2685	data.x2[i] = xs2;
		2686	data.y2[i] = ys2;
		2687	data.z2[i] = zs2;
		2688	data.bi[i] = radz;
		2689	delz = delz * rd;
		2690	radz = radz * rrad;
		2691	xs1 = xs2;
		2692	ys1 = ys2;
		2693	zs1 = zs2;
		2694	}
		2695
		2696	data.x2[data.n - 1] = xw2;
		2697	data.y2[data.n - 1] = yw2;
		2698	data.z2[data.n - 1] = zw2;
		2699
		2700	return;
		2701	}
		2702
		2703	/-----------------------------------------------------------------------/

Subversion Repositories Nec2c

(root)/trunk/src/geometry.c – Rev 2