WebSVN – Nec2c – Blame – /trunk/src/fields.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 /dataj/ */
		29	dataj_t dataj;
		30
		31	/* common /gnd/ */
		32	extern gnd_t gnd;
		33
		34	/* common /incom/ */
		35	extern incom_t incom;
		36
		37	/* common /tmi/ */
		38	tmi_t tmi;
		39
		40	/common /tmh/ /
		41	static tmh_t tmh;
		42
		43	/* common /gwav/ */
		44	extern gwav_t gwav;
		45
		46	/* common /data/ */
		47	extern data_t data;
		48
		49	/* common /crnt/ */
		50	extern crnt_t crnt;
		51
		52	/* common /fpat/ */
		53	extern fpat_t fpat;
		54
		55	/* common /plot/ */
		56	extern plot_t plot;
		57
		58	/* pointers to input/output files */
		59	extern FILE input_fp, output_fp, *plot_fp;
		60
		61	/-------------------------------------------------------------------/
		62
		63	/* compute near e fields of a segment with sine, cosine, and */
		64	/* constant currents. ground effect included. */
		65	void efld (double xi, double yi, double zi, double ai, int ij)
		66	{
		67	#define txk egnd[0]
		68	#define tyk egnd[1]
		69	#define tzk egnd[2]
		70	#define txs egnd[3]
		71	#define tys egnd[4]
		72	#define tzs egnd[5]
		73	#define txc egnd[6]
		74	#define tyc egnd[7]
		75	#define tzc egnd[8]
		76
		77	int ip;
		78	double xij, yij, ijx, rfl, salpr, zij, zp, rhox;
		79	double rhoy, rhoz, rh, r, rmag, cth, px, py;
		80	double xymag, xspec, yspec, rhospc, dmin, shaf;
		81	complex double epx, epy, refs, refps, zrsin, zratx, zscrn;
		82	complex double tezs, ters, tezc, terc, tezk, terk, egnd[9];
		83
		84	xij = xi - dataj.xj;
		85	yij = yi - dataj.yj;
		86	ijx = ij;
		87	rfl = -1.;
		88
		89	for (ip = 0; ip < gnd.ksymp; ip++)
		90	{
		91	if (ip == 1)
		92	ijx = 1;
		93	rfl = -rfl;
		94	salpr = dataj.salpj * rfl;
		95	zij = zi - rfl * dataj.zj;
		96	zp = xij * dataj.cabj + yij * dataj.sabj + zij * salpr;
		97	rhox = xij - dataj.cabj * zp;
		98	rhoy = yij - dataj.sabj * zp;
		99	rhoz = zij - salpr * zp;
		100
		101	rh = sqrtl (rhox * rhox + rhoy * rhoy + rhoz * rhoz + ai * ai);
		102	if (rh <= 1.e-10)
		103	{
		104	rhox = 0.;
		105	rhoy = 0.;
		106	rhoz = 0.;
		107	}
		108	else
		109	{
		110	rhox = rhox / rh;
		111	rhoy = rhoy / rh;
		112	rhoz = rhoz / rh;
		113	}
		114
		115	/* lumped current element approx. for large separations */
		116	r = sqrtl (zp * zp + rh * rh);
		117	if (r >= dataj.rkh)
		118	{
		119	rmag = TP * r;
		120	cth = zp / r;
		121	px = rh / r;
		122	txk = cmplx (cosl (rmag), -sinl (rmag));
		123	py = TP * r * r;
		124	tyk = ETA * cth * txk * cmplx (1.0, -1.0 / rmag) / py;
		125	tzk = ETA * px * txk * cmplx (1.0, rmag - 1.0 / rmag) / (2. * py);
		126	tezk = tyk * cth - tzk * px;
		127	terk = tyk * px + tzk * cth;
		128	rmag = sinl (PI * dataj.s) / PI;
		129	tezc = tezk * rmag;
		130	terc = terk * rmag;
		131	tezk = tezk * dataj.s;
		132	terk = terk * dataj.s;
		133	txs = CPLX_00;
		134	tys = CPLX_00;
		135	tzs = CPLX_00;
		136
		137	} /* if( r >= dataj.rkh) */
		138
		139	if (r < dataj.rkh)
		140	{
		141	/* eksc for thin wire approx. or ekscx for extended t.w. approx. */
		142	if (dataj.iexk != 1)
		143	eksc (
		144	dataj.s,
		145	zp,
		146	rh,
		147	TP,
		148	ijx,
		149	&tezs,
		150	&ters,
		151	&tezc,
		152	&terc,
		153	&tezk,
		154	&terk);
		155	else
		156	ekscx (
		157	dataj.b,
		158	dataj.s,
		159	zp,
		160	rh,
		161	TP,
		162	ijx,
		163	dataj.ind1,
		164	dataj.ind2,
		165	&tezs,
		166	&ters,
		167	&tezc,
		168	&terc,
		169	&tezk,
		170	&terk);
		171
		172	txs = tezs * dataj.cabj + ters * rhox;
		173	tys = tezs * dataj.sabj + ters * rhoy;
		174	tzs = tezs * salpr + ters * rhoz;
		175
		176	} /* if( r < dataj.rkh) */
		177
		178	txk = tezk * dataj.cabj + terk * rhox;
		179	tyk = tezk * dataj.sabj + terk * rhoy;
		180	tzk = tezk * salpr + terk * rhoz;
		181	txc = tezc * dataj.cabj + terc * rhox;
		182	tyc = tezc * dataj.sabj + terc * rhoy;
		183	tzc = tezc * salpr + terc * rhoz;
		184
		185	if (ip == 1)
		186	{
		187	if (gnd.iperf <= 0)
		188	{
		189	zratx = gnd.zrati;
		190	rmag = r;
		191	xymag = sqrtl (xij * xij + yij * yij);
		192
		193	/* set parameters for radial wire ground screen. */
		194	if (gnd.nradl != 0)
		195	{
		196	xspec =
		197	(xi * dataj.zj + zi * dataj.xj) / (zi + dataj.zj);
		198	yspec =
		199	(yi * dataj.zj + zi * dataj.yj) / (zi + dataj.zj);
		200	rhospc = sqrtl (
		201	xspec * xspec + yspec * yspec + gnd.t2 * gnd.t2);
		202
		203	if (rhospc <= gnd.scrwl)
		204	{
		205	zscrn =
		206	gnd.t1 * rhospc * logl (rhospc / gnd.t2);
		207	zratx = (zscrn * gnd.zrati) /
		208	(ETA * gnd.zrati + zscrn);
		209	}
		210	} /* if( gnd.nradl != 0) */
		211
		212	/* calculation of reflection coefficients when ground is
		213	* specified. */
		214	if (xymag <= 1.0e-6)
		215	{
		216	px = 0.;
		217	py = 0.;
		218	cth = 1.;
		219	zrsin = CPLX_10;
		220	}
		221	else
		222	{
		223	px = -yij / xymag;
		224	py = xij / xymag;
		225	cth = zij / rmag;
		226	zrsin =
		227	csqrtl (1.0 - zratx * zratx * (1.0 - cth * cth));
		228
		229	} /* if( xymag <= 1.0e-6) */
		230
		231	refs = (cth - zratx * zrsin) / (cth + zratx * zrsin);
		232	refps = -(zratx * cth - zrsin) / (zratx * cth + zrsin);
		233	refps = refps - refs;
		234	epy = px * txk + py * tyk;
		235	epx = px * epy;
		236	epy = py * epy;
		237	txk = refs * txk + refps * epx;
		238	tyk = refs * tyk + refps * epy;
		239	tzk = refs * tzk;
		240	epy = px * txs + py * tys;
		241	epx = px * epy;
		242	epy = py * epy;
		243	txs = refs * txs + refps * epx;
		244	tys = refs * tys + refps * epy;
		245	tzs = refs * tzs;
		246	epy = px * txc + py * tyc;
		247	epx = px * epy;
		248	epy = py * epy;
		249	txc = refs * txc + refps * epx;
		250	tyc = refs * tyc + refps * epy;
		251	tzc = refs * tzc;
		252
		253	} /* if( gnd.iperf <= 0) */
		254
		255	dataj.exk = dataj.exk - txk * gnd.frati;
		256	dataj.eyk = dataj.eyk - tyk * gnd.frati;
		257	dataj.ezk = dataj.ezk - tzk * gnd.frati;
		258	dataj.exs = dataj.exs - txs * gnd.frati;
		259	dataj.eys = dataj.eys - tys * gnd.frati;
		260	dataj.ezs = dataj.ezs - tzs * gnd.frati;
		261	dataj.exc = dataj.exc - txc * gnd.frati;
		262	dataj.eyc = dataj.eyc - tyc * gnd.frati;
		263	dataj.ezc = dataj.ezc - tzc * gnd.frati;
		264	continue;
		265
		266	} /* if( ip == 1) */
		267
		268	dataj.exk = txk;
		269	dataj.eyk = tyk;
		270	dataj.ezk = tzk;
		271	dataj.exs = txs;
		272	dataj.eys = tys;
		273	dataj.ezs = tzs;
		274	dataj.exc = txc;
		275	dataj.eyc = tyc;
		276	dataj.ezc = tzc;
		277
		278	} /* for( ip = 0; ip < gnd.ksymp; ip++ ) */
		279
		280	if (gnd.iperf != 2)
		281	return;
		282
		283	/* field due to ground using sommerfeld/norton */
		284	incom.sn = sqrtl (dataj.cabj * dataj.cabj + dataj.sabj * dataj.sabj);
		285	if (incom.sn >= 1.0e-5)
		286	{
		287	incom.xsn = dataj.cabj / incom.sn;
		288	incom.ysn = dataj.sabj / incom.sn;
		289	}
		290	else
		291	{
		292	incom.sn = 0.;
		293	incom.xsn = 1.;
		294	incom.ysn = 0.;
		295	}
		296
		297	/* displace observation point for thin wire approximation */
		298	zij = zi + dataj.zj;
		299	salpr = -dataj.salpj;
		300	rhox = dataj.sabj * zij - salpr * yij;
		301	rhoy = salpr * xij - dataj.cabj * zij;
		302	rhoz = dataj.cabj * yij - dataj.sabj * xij;
		303	rh = rhox * rhox + rhoy * rhoy + rhoz * rhoz;
		304
		305	if (rh <= 1.e-10)
		306	{
		307	incom.xo = xi - ai * incom.ysn;
		308	incom.yo = yi + ai * incom.xsn;
		309	incom.zo = zi;
		310	}
		311	else
		312	{
		313	rh = ai / sqrtl (rh);
		314	if (rhoz < 0.)
		315	rh = -rh;
		316	incom.xo = xi + rh * rhox;
		317	incom.yo = yi + rh * rhoy;
		318	incom.zo = zi + rh * rhoz;
		319
		320	} /* if( rh <= 1.e-10) */
		321
		322	r = xij * xij + yij * yij + zij * zij;
		323	if (r <= .95)
		324	{
		325	/* field from interpolation is integrated over segment */
		326	incom.isnor = 1;
		327	dmin = dataj.exk * conjl (dataj.exk) + dataj.eyk * conjl (dataj.eyk) +
		328	dataj.ezk * conjl (dataj.ezk);
		329	dmin = .01 * sqrtl (dmin);
		330	shaf = .5 * dataj.s;
		331	rom2 (-shaf, shaf, egnd, dmin);
		332	}
		333	else
		334	{
		335	/* norton field equations and lumped current element approximation */
		336	incom.isnor = 2;
		337	sflds (0., egnd);
		338	} /* if( r <= .95) */
		339
		340	if (r > .95)
		341	{
		342	zp = xij * dataj.cabj + yij * dataj.sabj + zij * salpr;
		343	rh = r - zp * zp;
		344	if (rh <= 1.e-10)
		345	dmin = 0.;
		346	else
		347	dmin = sqrtl (rh / (rh + ai * ai));
		348
		349	if (dmin <= .95)
		350	{
		351	px = 1. - dmin;
		352	terk = (txk * dataj.cabj + tyk * dataj.sabj + tzk * salpr) * px;
		353	txk = dmin * txk + terk * dataj.cabj;
		354	tyk = dmin * tyk + terk * dataj.sabj;
		355	tzk = dmin * tzk + terk * salpr;
		356	ters = (txs * dataj.cabj + tys * dataj.sabj + tzs * salpr) * px;
		357	txs = dmin * txs + ters * dataj.cabj;
		358	tys = dmin * tys + ters * dataj.sabj;
		359	tzs = dmin * tzs + ters * salpr;
		360	terc = (txc * dataj.cabj + tyc * dataj.sabj + tzc * salpr) * px;
		361	txc = dmin * txc + terc * dataj.cabj;
		362	tyc = dmin * tyc + terc * dataj.sabj;
		363	tzc = dmin * tzc + terc * salpr;
		364
		365	} /* if( dmin <= .95) */
		366
		367	} /* if( r > .95) */
		368
		369	dataj.exk = dataj.exk + txk;
		370	dataj.eyk = dataj.eyk + tyk;
		371	dataj.ezk = dataj.ezk + tzk;
		372	dataj.exs = dataj.exs + txs;
		373	dataj.eys = dataj.eys + tys;
		374	dataj.ezs = dataj.ezs + tzs;
		375	dataj.exc = dataj.exc + txc;
		376	dataj.eyc = dataj.eyc + tyc;
		377	dataj.ezc = dataj.ezc + tzc;
		378
		379	return;
		380	}
		381
		382	/-----------------------------------------------------------------------/
		383
		384	/* compute e field of sine, cosine, and constant */
		385	/* current filaments by thin wire approximation. */
		386	void eksc (
		387	double s,
		388	double z,
		389	double rh,
		390	double xk,
		391	int ij,
		392	complex double *ezs,
		393	complex double *ers,
		394	complex double *ezc,
		395	complex double *erc,
		396	complex double *ezk,
		397	complex double *erk)
		398	{
		399	double rhk, sh, shk, ss, cs, z1a, z2a, cint, sint;
		400	complex double gz1, gz2, gp1, gp2, gzp1, gzp2;
		401
		402	tmi.ij = ij;
		403	tmi.zpk = xk * z;
		404	rhk = xk * rh;
		405	tmi.rkb2 = rhk * rhk;
		406	sh = .5 * s;
		407	shk = xk * sh;
		408	ss = sinl (shk);
		409	cs = cosl (shk);
		410	z2a = sh - z;
		411	z1a = -(sh + z);
		412	gx (z1a, rh, xk, &gz1, &gp1);
		413	gx (z2a, rh, xk, &gz2, &gp2);
		414	gzp1 = gp1 * z1a;
		415	gzp2 = gp2 * z2a;
		416	ezs = CONST1 ((gz2 - gz1) * cs * xk - (gzp2 + gzp1) * ss);
		417	ezc = -CONST1 ((gz2 + gz1) * ss * xk + (gzp2 - gzp1) * cs);
		418	erk = CONST1 (gp2 - gp1) * rh;
		419	intx (-shk, shk, rhk, ij, &cint, &sint);
		420	ezk = -CONST1 (gzp2 - gzp1 + xk * xk * cmplx (cint, -sint));
		421	gzp1 = gzp1 * z1a;
		422	gzp2 = gzp2 * z2a;
		423
		424	if (rh >= 1.0e-10)
		425	{
		426	ers = -CONST1
		427	((gzp2 + gzp1 + gz2 + gz1) * ss - (z2a * gz2 - z1a * gz1) * cs * xk) /
		428	rh;
		429	erc = -CONST1
		430	((gzp2 - gzp1 + gz2 - gz1) * cs + (z2a * gz2 + z1a * gz1) * ss * xk) /
		431	rh;
		432	return;
		433	}
		434
		435	*ers = CPLX_00;
		436	*erc = CPLX_00;
		437
		438	return;
		439	}
		440
		441	/-----------------------------------------------------------------------/
		442
		443	/* compute e field of sine, cosine, and constant current */
		444	/* filaments by extended thin wire approximation. */
		445	void ekscx (
		446	double bx,
		447	double s,
		448	double z,
		449	double rhx,
		450	double xk,
		451	int ij,
		452	int inx1,
		453	int inx2,
		454	complex double *ezs,
		455	complex double *ers,
		456	complex double *ezc,
		457	complex double *erc,
		458	complex double *ezk,
		459	complex double *erk)
		460	{
		461	int ira;
		462	double b, rh, sh, rhk, shk, ss, cs, z1a;
		463	double z2a, a2, bk, bk2, cint, sint;
		464	complex double gz1, gz2, gzp1, gzp2, gr1, gr2;
		465	complex double grp1, grp2, grk1, grk2, gzz1, gzz2;
		466
		467	if (rhx >= bx)
		468	{
		469	rh = rhx;
		470	b = bx;
		471	ira = 0;
		472	}
		473	else
		474	{
		475	rh = bx;
		476	b = rhx;
		477	ira = 1;
		478	}
		479
		480	sh = .5 * s;
		481	tmi.ij = ij;
		482	tmi.zpk = xk * z;
		483	rhk = xk * rh;
		484	tmi.rkb2 = rhk * rhk;
		485	shk = xk * sh;
		486	ss = sinl (shk);
		487	cs = cosl (shk);
		488	z2a = sh - z;
		489	z1a = -(sh + z);
		490	a2 = b * b;
		491
		492	if (inx1 != 2)
		493	gxx (z1a, rh, b, a2, xk, ira, &gz1, &gzp1, &gr1, &grp1, &grk1, &gzz1);
		494	else
		495	{
		496	gx (z1a, rhx, xk, &gz1, &grk1);
		497	gzp1 = grk1 * z1a;
		498	gr1 = gz1 / rhx;
		499	grp1 = gzp1 / rhx;
		500	grk1 = grk1 * rhx;
		501	gzz1 = CPLX_00;
		502	}
		503
		504	if (inx2 != 2)
		505	gxx (z2a, rh, b, a2, xk, ira, &gz2, &gzp2, &gr2, &grp2, &grk2, &gzz2);
		506	else
		507	{
		508	gx (z2a, rhx, xk, &gz2, &grk2);
		509	gzp2 = grk2 * z2a;
		510	gr2 = gz2 / rhx;
		511	grp2 = gzp2 / rhx;
		512	grk2 = grk2 * rhx;
		513	gzz2 = CPLX_00;
		514	}
		515
		516	ezs = CONST1 ((gz2 - gz1) * cs * xk - (gzp2 + gzp1) * ss);
		517	ezc = -CONST1 ((gz2 + gz1) * ss * xk + (gzp2 - gzp1) * cs);
		518	ers = -CONST1 ((z2a * grp2 + z1a * grp1 + gr2 + gr1) * ss -
		519	(z2a * gr2 - z1a * gr1) * cs * xk);
		520	erc = -CONST1 ((z2a * grp2 - z1a * grp1 + gr2 - gr1) * cs +
		521	(z2a * gr2 + z1a * gr1) * ss * xk);
		522	erk = CONST1 (grk2 - grk1);
		523	intx (-shk, shk, rhk, ij, &cint, &sint);
		524	bk = b * xk;
		525	bk2 = bk * bk * .25;
		526	ezk = -CONST1 (gzp2 - gzp1 + xk * xk * (1. - bk2) * cmplx (cint, -sint) -
		527	bk2 * (gzz2 - gzz1));
		528
		529	return;
		530	}
		531
		532	/-----------------------------------------------------------------------/
		533
		534	/* integrand for h field of a wire */
		535	void gh (double zk, double hr, double hi)
		536	{
		537	double rs, r, ckr, skr, rr2, rr3;
		538
		539	rs = zk - tmh.zpka;
		540	rs = tmh.rhks + rs * rs;
		541	r = sqrtl (rs);
		542	ckr = cosl (r);
		543	skr = sinl (r);
		544	rr2 = 1. / rs;
		545	rr3 = rr2 / r;
		546	hr = skr rr2 + ckr * rr3;
		547	hi = ckr rr2 - skr * rr3;
		548
		549	return;
		550	}
		551
		552	/-----------------------------------------------------------------------/
		553
		554	/* gwave computes the electric field, including ground wave, of a */
		555	/* current element over a ground plane using formulas of k.a. norton */
		556	/* (proc. ire, sept., 1937, pp.1203,1236.) */
		557
		558	void gwave (
		559	complex double *erv,
		560	complex double *ezv,
		561	complex double *erh,
		562	complex double *ezh,
		563	complex double *eph)
		564	{
		565	double sppp, sppp2, cppp2, cppp, spp, spp2, cpp2, cpp;
		566	complex double rk1, rk2, t1, t2, t3, t4, p1, rv;
		567	complex double omr, w, f, q1, rh, v, g, xr1, xr2;
		568	complex double x1, x2, x3, x4, x5, x6, x7;
		569
		570	sppp = gwav.zmh / gwav.r1;
		571	sppp2 = sppp * sppp;
		572	cppp2 = 1. - sppp2;
		573
		574	if (cppp2 < 1.0e-20)
		575	cppp2 = 1.0e-20;
		576
		577	cppp = sqrtl (cppp2);
		578	spp = gwav.zph / gwav.r2;
		579	spp2 = spp * spp;
		580	cpp2 = 1. - spp2;
		581
		582	if (cpp2 < 1.0e-20)
		583	cpp2 = 1.0e-20;
		584
		585	cpp = sqrtl (cpp2);
		586	rk1 = -TPJ * gwav.r1;
		587	rk2 = -TPJ * gwav.r2;
		588	t1 = 1. - gwav.u2 * cpp2;
		589	t2 = csqrtl (t1);
		590	t3 = (1. - 1. / rk1) / rk1;
		591	t4 = (1. - 1. / rk2) / rk2;
		592	p1 = rk2 * gwav.u2 * t1 / (2. * cpp2);
		593	rv = (spp - gwav.u * t2) / (spp + gwav.u * t2);
		594	omr = 1. - rv;
		595	w = 1. / omr;
		596	w = (4.0 + 0.0fj) * p1 * w * w;
		597	fbar (w, &f);
		598	q1 = rk2 * t1 / (2. * gwav.u2 * cpp2);
		599	rh = (t2 - gwav.u * spp) / (t2 + gwav.u * spp);
		600	v = 1. / (1. + rh);
		601	v = (4.0 + 0.0fj) * q1 * v * v;
		602	fbar (v, &g);
		603	xr1 = gwav.xx1 / gwav.r1;
		604	xr2 = gwav.xx2 / gwav.r2;
		605	x1 = cppp2 * xr1;
		606	x2 = rv * cpp2 * xr2;
		607	x3 = omr * cpp2 * f * xr2;
		608	x4 = gwav.u * t2 * spp * 2. * xr2 / rk2;
		609	x5 = xr1 * t3 * (1. - 3. * sppp2);
		610	x6 = xr2 * t4 * (1. - 3. * spp2);
		611	ezv = (x1 + x2 + x3 - x4 - x5 - x6) (-CONST4);
		612	x1 = sppp * cppp * xr1;
		613	x2 = rv * spp * cpp * xr2;
		614	x3 = cpp * omr * gwav.u * t2 * f * xr2;
		615	x4 = spp * cpp * omr * xr2 / rk2;
		616	x5 = 3. * sppp * cppp * t3 * xr1;
		617	x6 = cpp * gwav.u * t2 * omr * xr2 / rk2 * .5;
		618	x7 = 3. * spp * cpp * t4 * xr2;
		619	erv = -(x1 + x2 - x3 + x4 - x5 + x6 - x7) (-CONST4);
		620	ezh = -(x1 - x2 + x3 - x4 - x5 - x6 + x7) (-CONST4);
		621	x1 = sppp2 * xr1;
		622	x2 = rv * spp2 * xr2;
		623	x4 = gwav.u2 * t1 * omr * f * xr2;
		624	x5 = t3 * (1. - 3. * cppp2) * xr1;
		625	x6 = t4 * (1. - 3. * cpp2) * (1. - gwav.u2 * (1. + rv) - gwav.u2 * omr * f) * xr2;
		626	x7 = gwav.u2 * cpp2 * omr * (1. - 1. / rk2) *
		627	(f * (gwav.u2 * t1 - spp2 - 1. / rk2) + 1. / rk2) * xr2;
		628	erh = (x1 - x2 - x4 - x5 + x6 + x7) (-CONST4);
		629	x1 = xr1;
		630	x2 = rh * xr2;
		631	x3 = (rh + 1.) * g * xr2;
		632	x4 = t3 * xr1;
		633	x5 = t4 * (1. - gwav.u2 * (1. + rv) - gwav.u2 * omr * f) * xr2;
		634	x6 =
		635	.5 * gwav.u2 * omr * (f * (gwav.u2 * t1 - spp2 - 1. / rk2) + 1. / rk2) * xr2 / rk2;
		636	eph = -(x1 - x2 + x3 - x4 + x5 + x6) (-CONST4);
		637
		638	return;
		639	}
		640
		641	/-----------------------------------------------------------------------/
		642
		643	/* segment end contributions for thin wire approx. */
		644	void gx (double zz, double rh, double xk, complex double gz, complex double gzp)
		645	{
		646	double r, r2, rkz;
		647
		648	r2 = zz * zz + rh * rh;
		649	r = sqrtl (r2);
		650	rkz = xk * r;
		651	*gz = cmplx (cosl (rkz), -sinl (rkz)) / r;
		652	gzp = -cmplx (1.0, rkz) *gz / r2;
		653
		654	return;
		655	}
		656
		657	/-----------------------------------------------------------------------/
		658
		659	/* segment end contributions for ext. thin wire approx. */
		660	void gxx (
		661	double zz,
		662	double rh,
		663	double a,
		664	double a2,
		665	double xk,
		666	int ira,
		667	complex double *g1,
		668	complex double *g1p,
		669	complex double *g2,
		670	complex double *g2p,
		671	complex double *g3,
		672	complex double *gzp)
		673	{
		674	double r, r2, r4, rk, rk2, rh2, t1, t2;
		675	complex double gz, c1, c2, c3;
		676
		677	r2 = zz * zz + rh * rh;
		678	r = sqrtl (r2);
		679	r4 = r2 * r2;
		680	rk = xk * r;
		681	rk2 = rk * rk;
		682	rh2 = rh * rh;
		683	t1 = .25 * a2 * rh2 / r4;
		684	t2 = .5 * a2 / r2;
		685	c1 = cmplx (1.0, rk);
		686	c2 = 3. * c1 - rk2;
		687	c3 = cmplx (6.0, rk) * rk2 - 15. * c1;
		688	gz = cmplx (cosl (rk), -sinl (rk)) / r;
		689	g2 = gz (1. + t1 * c2);
		690	g1 = g2 - t2 * c1 * gz;
		691	gz = gz / r2;
		692	g2p = gz (t1 * c3 - c1);
		693	gzp = t2 c2 * gz;
		694	g3 = g2p + *gzp;
		695	g1p = g3 * zz;
		696
		697	if (ira != 1)
		698	{
		699	g3 = (g3 + gzp) rh;
		700	gzp = -zz c1 * gz;
		701
		702	if (rh <= 1.0e-10)
		703	{
		704	*g2 = 0.;
		705	*g2p = 0.;
		706	return;
		707	}
		708
		709	g2 = g2 / rh;
		710	g2p = g2p * zz / rh;
		711	return;
		712
		713	} /* if( ira != 1) */
		714
		715	t2 = .5 * a;
		716	g2 = -t2 c1 * gz;
		717	g2p = t2 gz * c2 / r2;
		718	g3 = rh2 g2p - a gz * c1;
		719	g2p = g2p * zz;
		720	gzp = -zz c1 * gz;
		721
		722	return;
		723	}
		724
		725	/-----------------------------------------------------------------------/
		726
		727	/* hfk computes the h field of a uniform current */
		728	/* filament by numerical integration */
		729	void hfk (double el1, double el2, double rhk, double zpkx, double sgr, double sgi)
		730	{
		731	int nx = 1, nma = 65536, nts = 4;
		732	int ns, nt;
		733	int flag = TRUE;
		734	double rx = 1.0e-4;
		735	double z, ze, s, ep, zend, dz = 0., zp, dzot = 0., t00r, g1r, g5r = 0, t00i;
		736	double g1i, g5i = 0., t01r, g3r = 0, t01i, g3i = 0, t10r, t10i, te1i, te1r, t02r;
		737	double g2r, g4r, t02i, g2i, g4i, t11r, t11i, t20r, t20i, te2i, te2r;
		738
		739	tmh.zpka = zpkx;
		740	tmh.rhks = rhk * rhk;
		741	z = el1;
		742	ze = el2;
		743	s = ze - z;
		744	ep = s / (10. * nma);
		745	zend = ze - ep;
		746	*sgr = 0.0;
		747	*sgi = 0.0;
		748	ns = nx;
		749	nt = 0;
		750	gh (z, &g1r, &g1i);
		751
		752	while (TRUE)
		753	{
		754	if (flag)
		755	{
		756	dz = s / ns;
		757	zp = z + dz;
		758
		759	if (zp > ze)
		760	{
		761	dz = ze - z;
		762	if (fabsl (dz) <= ep)
		763	{
		764	sgr = sgr * rhk * .5;
		765	sgi = sgi * rhk * .5;
		766	return;
		767	}
		768	}
		769
		770	dzot = dz * .5;
		771	zp = z + dzot;
		772	gh (zp, &g3r, &g3i);
		773	zp = z + dz;
		774	gh (zp, &g5r, &g5i);
		775
		776	} /* if( flag ) */
		777
		778	t00r = (g1r + g5r) * dzot;
		779	t00i = (g1i + g5i) * dzot;
		780	t01r = (t00r + dz * g3r) * 0.5;
		781	t01i = (t00i + dz * g3i) * 0.5;
		782	t10r = (4.0 * t01r - t00r) / 3.0;
		783	t10i = (4.0 * t01i - t00i) / 3.0;
		784
		785	test (t01r, t10r, &te1r, t01i, t10i, &te1i, 0.);
		786	if ((te1i <= rx) && (te1r <= rx))
		787	{
		788	sgr = sgr + t10r;
		789	sgi = sgi + t10i;
		790	nt += 2;
		791
		792	z += dz;
		793	if (z >= zend)
		794	{
		795	sgr = sgr * rhk * .5;
		796	sgi = sgi * rhk * .5;
		797	return;
		798	}
		799
		800	g1r = g5r;
		801	g1i = g5i;
		802	if (nt >= nts)
		803	if (ns > nx)
		804	{
		805	ns = ns / 2;
		806	nt = 1;
		807	}
		808	flag = TRUE;
		809	continue;
		810
		811	} /* if( (te1i <= rx) && (te1r <= rx) ) */
		812
		813	zp = z + dz * 0.25;
		814	gh (zp, &g2r, &g2i);
		815	zp = z + dz * 0.75;
		816	gh (zp, &g4r, &g4i);
		817	t02r = (t01r + dzot * (g2r + g4r)) * 0.5;
		818	t02i = (t01i + dzot * (g2i + g4i)) * 0.5;
		819	t11r = (4.0 * t02r - t01r) / 3.0;
		820	t11i = (4.0 * t02i - t01i) / 3.0;
		821	t20r = (16.0 * t11r - t10r) / 15.0;
		822	t20i = (16.0 * t11i - t10i) / 15.0;
		823
		824	test (t11r, t20r, &te2r, t11i, t20i, &te2i, 0.);
		825	if ((te2i > rx) \|\| (te2r > rx))
		826	{
		827	nt = 0;
		828	if (ns >= nma)
		829	fprintf (output_fp, "\n STEP SIZE LIMITED AT Z= %10.5f", z);
		830	else
		831	{
		832	ns = ns * 2;
		833	dz = s / ns;
		834	dzot = dz * 0.5;
		835	g5r = g3r;
		836	g5i = g3i;
		837	g3r = g2r;
		838	g3i = g2i;
		839
		840	flag = FALSE;
		841	continue;
		842	}
		843
		844	} /* if( (te2i > rx) \|\| (te2r > rx) ) */
		845
		846	sgr = sgr + t20r;
		847	sgi = sgi + t20i;
		848	nt++;
		849
		850	z += dz;
		851	if (z >= zend)
		852	{
		853	sgr = sgr * rhk * .5;
		854	sgi = sgi * rhk * .5;
		855	return;
		856	}
		857
		858	g1r = g5r;
		859	g1i = g5i;
		860	if (nt >= nts)
		861	if (ns > nx)
		862	{
		863	ns = ns / 2;
		864	nt = 1;
		865	}
		866	flag = TRUE;
		867
		868	} /* while( TRUE ) */
		869	}
		870
		871	/-----------------------------------------------------------------------/
		872
		873	/* hintg computes the h field of a patch current */
		874	void hintg (double xi, double yi, double zi)
		875	{
		876	int ip;
		877	double rx, ry, rfl, xymag, pxx, pyy, cth;
		878	double rz, rsq, r, rk, cr, sr, t1zr, t2zr;
		879	complex double gam, f1x, f1y, f1z, f2x, f2y, f2z, rrv, rrh;
		880
		881	rx = xi - dataj.xj;
		882	ry = yi - dataj.yj;
		883	rfl = -1.;
		884	dataj.exk = CPLX_00;
		885	dataj.eyk = CPLX_00;
		886	dataj.ezk = CPLX_00;
		887	dataj.exs = CPLX_00;
		888	dataj.eys = CPLX_00;
		889	dataj.ezs = CPLX_00;
		890
		891	for (ip = 1; ip <= gnd.ksymp; ip++)
		892	{
		893	rfl = -rfl;
		894	rz = zi - dataj.zj * rfl;
		895	rsq = rx * rx + ry * ry + rz * rz;
		896
		897	if (rsq < 1.0e-20)
		898	continue;
		899
		900	r = sqrtl (rsq);
		901	rk = TP * r;
		902	cr = cosl (rk);
		903	sr = sinl (rk);
		904	gam = -(cmplx (cr, -sr) + rk * cmplx (sr, cr)) / (FPI * rsq * r) * dataj.s;
		905	dataj.exc = gam * rx;
		906	dataj.eyc = gam * ry;
		907	dataj.ezc = gam * rz;
		908	t1zr = dataj.t1zj * rfl;
		909	t2zr = dataj.t2zj * rfl;
		910	f1x = dataj.eyc * t1zr - dataj.ezc * dataj.t1yj;
		911	f1y = dataj.ezc * dataj.t1xj - dataj.exc * t1zr;
		912	f1z = dataj.exc * dataj.t1yj - dataj.eyc * dataj.t1xj;
		913	f2x = dataj.eyc * t2zr - dataj.ezc * dataj.t2yj;
		914	f2y = dataj.ezc * dataj.t2xj - dataj.exc * t2zr;
		915	f2z = dataj.exc * dataj.t2yj - dataj.eyc * dataj.t2xj;
		916
		917	if (ip != 1)
		918	{
		919	if (gnd.iperf == 1)
		920	{
		921	f1x = -f1x;
		922	f1y = -f1y;
		923	f1z = -f1z;
		924	f2x = -f2x;
		925	f2y = -f2y;
		926	f2z = -f2z;
		927	}
		928	else
		929	{
		930	xymag = sqrtl (rx * rx + ry * ry);
		931	if (xymag <= 1.0e-6)
		932	{
		933	pxx = 0.;
		934	pyy = 0.;
		935	cth = 1.;
		936	rrv = CPLX_10;
		937	}
		938	else
		939	{
		940	pxx = -ry / xymag;
		941	pyy = rx / xymag;
		942	cth = rz / r;
		943	rrv = csqrtl (
		944	1. - gnd.zrati * gnd.zrati * (1. - cth * cth));
		945
		946	} /* if( xymag <= 1.0e-6) */
		947
		948	rrh = gnd.zrati * cth;
		949	rrh = (rrh - rrv) / (rrh + rrv);
		950	rrv = gnd.zrati * rrv;
		951	rrv = -(cth - rrv) / (cth + rrv);
		952	gam = (f1x * pxx + f1y * pyy) * (rrv - rrh);
		953	f1x = f1x * rrh + gam * pxx;
		954	f1y = f1y * rrh + gam * pyy;
		955	f1z = f1z * rrh;
		956	gam = (f2x * pxx + f2y * pyy) * (rrv - rrh);
		957	f2x = f2x * rrh + gam * pxx;
		958	f2y = f2y * rrh + gam * pyy;
		959	f2z = f2z * rrh;
		960
		961	} /* if( gnd.iperf == 1) */
		962
		963	} /* if( ip != 1) */
		964
		965	dataj.exk += f1x;
		966	dataj.eyk += f1y;
		967	dataj.ezk += f1z;
		968	dataj.exs += f2x;
		969	dataj.eys += f2y;
		970	dataj.ezs += f2z;
		971
		972	} /* for( ip = 1; ip <= gnd.ksymp; ip++ ) */
		973
		974	return;
		975	}
		976
		977	/-----------------------------------------------------------------------/
		978
		979	/* hsfld computes the h field for constant, sine, and */
		980	/* cosine current on a segment including ground effects. */
		981	void hsfld (double xi, double yi, double zi, double ai)
		982	{
		983	int ip;
		984	double xij, yij, rfl, salpr, zij, zp, rhox, rhoy, rhoz, rh, phx;
		985	double phy, phz, rmag, xymag, xspec, yspec, rhospc, px, py, cth;
		986	complex double hpk, hps, hpc, qx, qy, qz, rrv, rrh, zratx;
		987
		988	xij = xi - dataj.xj;
		989	yij = yi - dataj.yj;
		990	rfl = -1.;
		991
		992	for (ip = 0; ip < gnd.ksymp; ip++)
		993	{
		994	rfl = -rfl;
		995	salpr = dataj.salpj * rfl;
		996	zij = zi - rfl * dataj.zj;
		997	zp = xij * dataj.cabj + yij * dataj.sabj + zij * salpr;
		998	rhox = xij - dataj.cabj * zp;
		999	rhoy = yij - dataj.sabj * zp;
		1000	rhoz = zij - salpr * zp;
		1001	rh = sqrtl (rhox * rhox + rhoy * rhoy + rhoz * rhoz + ai * ai);
		1002
		1003	if (rh <= 1.0e-10)
		1004	{
		1005	dataj.exk = 0.;
		1006	dataj.eyk = 0.;
		1007	dataj.ezk = 0.;
		1008	dataj.exs = 0.;
		1009	dataj.eys = 0.;
		1010	dataj.ezs = 0.;
		1011	dataj.exc = 0.;
		1012	dataj.eyc = 0.;
		1013	dataj.ezc = 0.;
		1014	continue;
		1015	}
		1016
		1017	rhox = rhox / rh;
		1018	rhoy = rhoy / rh;
		1019	rhoz = rhoz / rh;
		1020	phx = dataj.sabj * rhoz - salpr * rhoy;
		1021	phy = salpr * rhox - dataj.cabj * rhoz;
		1022	phz = dataj.cabj * rhoy - dataj.sabj * rhox;
		1023
		1024	hsflx (dataj.s, rh, zp, &hpk, &hps, &hpc);
		1025
		1026	if (ip == 1)
		1027	{
		1028	if (gnd.iperf != 1)
		1029	{
		1030	zratx = gnd.zrati;
		1031	rmag = sqrtl (zp * zp + rh * rh);
		1032	xymag = sqrtl (xij * xij + yij * yij);
		1033
		1034	/* set parameters for radial wire ground screen. */
		1035	if (gnd.nradl != 0)
		1036	{
		1037	xspec =
		1038	(xi * dataj.zj + zi * dataj.xj) / (zi + dataj.zj);
		1039	yspec =
		1040	(yi * dataj.zj + zi * dataj.yj) / (zi + dataj.zj);
		1041	rhospc = sqrtl (
		1042	xspec * xspec + yspec * yspec + gnd.t2 * gnd.t2);
		1043
		1044	if (rhospc <= gnd.scrwl)
		1045	{
		1046	rrv = gnd.t1 * rhospc * logl (rhospc / gnd.t2);
		1047	zratx = (rrv * gnd.zrati) /
		1048	(ETA * gnd.zrati + rrv);
		1049	}
		1050	}
		1051
		1052	/* calculation of reflection coefficients when ground is
		1053	* specified. */
		1054	if (xymag <= 1.0e-6)
		1055	{
		1056	px = 0.;
		1057	py = 0.;
		1058	cth = 1.;
		1059	rrv = CPLX_10;
		1060	}
		1061	else
		1062	{
		1063	px = -yij / xymag;
		1064	py = xij / xymag;
		1065	cth = zij / rmag;
		1066	rrv = csqrtl (1. - zratx * zratx * (1. - cth * cth));
		1067	}
		1068
		1069	rrh = zratx * cth;
		1070	rrh = -(rrh - rrv) / (rrh + rrv);
		1071	rrv = zratx * rrv;
		1072	rrv = (cth - rrv) / (cth + rrv);
		1073	qy = (phx * px + phy * py) * (rrv - rrh);
		1074	qx = qy * px + phx * rrh;
		1075	qy = qy * py + phy * rrh;
		1076	qz = phz * rrh;
		1077	dataj.exk = dataj.exk - hpk * qx;
		1078	dataj.eyk = dataj.eyk - hpk * qy;
		1079	dataj.ezk = dataj.ezk - hpk * qz;
		1080	dataj.exs = dataj.exs - hps * qx;
		1081	dataj.eys = dataj.eys - hps * qy;
		1082	dataj.ezs = dataj.ezs - hps * qz;
		1083	dataj.exc = dataj.exc - hpc * qx;
		1084	dataj.eyc = dataj.eyc - hpc * qy;
		1085	dataj.ezc = dataj.ezc - hpc * qz;
		1086	continue;
		1087
		1088	} /* if( gnd.iperf != 1 ) */
		1089
		1090	dataj.exk = dataj.exk - hpk * phx;
		1091	dataj.eyk = dataj.eyk - hpk * phy;
		1092	dataj.ezk = dataj.ezk - hpk * phz;
		1093	dataj.exs = dataj.exs - hps * phx;
		1094	dataj.eys = dataj.eys - hps * phy;
		1095	dataj.ezs = dataj.ezs - hps * phz;
		1096	dataj.exc = dataj.exc - hpc * phx;
		1097	dataj.eyc = dataj.eyc - hpc * phy;
		1098	dataj.ezc = dataj.ezc - hpc * phz;
		1099	continue;
		1100
		1101	} /* if( ip == 1 ) */
		1102
		1103	dataj.exk = hpk * phx;
		1104	dataj.eyk = hpk * phy;
		1105	dataj.ezk = hpk * phz;
		1106	dataj.exs = hps * phx;
		1107	dataj.eys = hps * phy;
		1108	dataj.ezs = hps * phz;
		1109	dataj.exc = hpc * phx;
		1110	dataj.eyc = hpc * phy;
		1111	dataj.ezc = hpc * phz;
		1112
		1113	} /* for( ip = 0; ip < gnd.ksymp; ip++ ) */
		1114
		1115	return;
		1116	}
		1117
		1118	/-----------------------------------------------------------------------/
		1119
		1120	/* calculates h field of sine cosine, and constant current of segment */
		1121	void hsflx (
		1122	double s,
		1123	double rh,
		1124	double zpx,
		1125	complex double *hpk,
		1126	complex double *hps,
		1127	complex double *hpc)
		1128	{
		1129	double r1, r2, zp, z2a, hss, dh, z1;
		1130	double rhz, dk, cdk, sdk, hkr, hki, rh2;
		1131	complex double fjk, ekr1, ekr2, t1, t2, cons;
		1132
		1133	fjk = -TPJ;
		1134	if (rh >= 1.0e-10)
		1135	{
		1136	if (zpx >= 0.)
		1137	{
		1138	zp = zpx;
		1139	hss = 1.;
		1140	}
		1141	else
		1142	{
		1143	zp = -zpx;
		1144	hss = -1.;
		1145	}
		1146
		1147	dh = .5 * s;
		1148	z1 = zp + dh;
		1149	z2a = zp - dh;
		1150	if (z2a >= 1.0e-7)
		1151	rhz = rh / z2a;
		1152	else
		1153	rhz = 1.;
		1154
		1155	dk = TP * dh;
		1156	cdk = cosl (dk);
		1157	sdk = sinl (dk);
		1158	hfk (-dk, dk, rh * TP, zp * TP, &hkr, &hki);
		1159	*hpk = cmplx (hkr, hki);
		1160
		1161	if (rhz >= 1.0e-3)
		1162	{
		1163	rh2 = rh * rh;
		1164	r1 = sqrtl (rh2 + z1 * z1);
		1165	r2 = sqrtl (rh2 + z2a * z2a);
		1166	ekr1 = cexp (fjk * r1);
		1167	ekr2 = cexp (fjk * r2);
		1168	t1 = z1 * ekr1 / r1;
		1169	t2 = z2a * ekr2 / r2;
		1170	hps = (cdk (ekr2 - ekr1) - CPLX_01 * sdk * (t2 + t1)) * hss;
		1171	hpc = -sdk (ekr2 + ekr1) - CPLX_01 * cdk * (t2 - t1);
		1172	cons = -CPLX_01 / (2. * TP * rh);
		1173	hps = cons *hps;
		1174	hpc = cons *hpc;
		1175	return;
		1176
		1177	} /* if( rhz >= 1.0e-3) */
		1178
		1179	ekr1 = cmplx (cdk, sdk) / (z2a * z2a);
		1180	ekr2 = cmplx (cdk, -sdk) / (z1 * z1);
		1181	t1 = TP * (1. / z1 - 1. / z2a);
		1182	t2 = cexp (fjk * zp) * rh / PI8;
		1183	hps = t2 (t1 + (ekr1 + ekr2) * sdk) * hss;
		1184	hpc = t2 (-CPLX_01 * t1 + (ekr1 - ekr2) * cdk);
		1185	return;
		1186
		1187	} /* if( rh >= 1.0e-10) */
		1188
		1189	*hps = CPLX_00;
		1190	*hpc = CPLX_00;
		1191	*hpk = CPLX_00;
		1192
		1193	return;
		1194	}
		1195
		1196	/-----------------------------------------------------------------------/
		1197
		1198	/* nefld computes the near field at specified points in space after */
		1199	/* the structure currents have been computed. */
		1200	void nefld (
		1201	double xob,
		1202	double yob,
		1203	double zob,
		1204	complex double *ex,
		1205	complex double *ey,
		1206	complex double *ez)
		1207	{
		1208	int i, ix, ipr, iprx, jc, ipa;
		1209	double zp, xi, ax;
		1210	complex double acx, bcx, ccx;
		1211
		1212	*ex = CPLX_00;
		1213	*ey = CPLX_00;
		1214	*ez = CPLX_00;
		1215	ax = 0.;
		1216
		1217	if (data.n != 0)
		1218	{
		1219	for (i = 0; i < data.n; i++)
		1220	{
		1221	dataj.xj = xob - data.x[i];
		1222	dataj.yj = yob - data.y[i];
		1223	dataj.zj = zob - data.z[i];
		1224	zp = data.cab[i] * dataj.xj + data.sab[i] * dataj.yj +
		1225	data.salp[i] * dataj.zj;
		1226
		1227	if (fabsl (zp) > 0.5001 * data.si[i])
		1228	continue;
		1229
		1230	zp = dataj.xj * dataj.xj + dataj.yj * dataj.yj + dataj.zj * dataj.zj -
		1231	zp * zp;
		1232	dataj.xj = data.bi[i];
		1233
		1234	if (zp > 0.9 * dataj.xj * dataj.xj)
		1235	continue;
		1236
		1237	ax = dataj.xj;
		1238	break;
		1239
		1240	} /* for( i = 0; i < n; i++ ) */
		1241
		1242	for (i = 0; i < data.n; i++)
		1243	{
		1244	ix = i + 1;
		1245	dataj.s = data.si[i];
		1246	dataj.b = data.bi[i];
		1247	dataj.xj = data.x[i];
		1248	dataj.yj = data.y[i];
		1249	dataj.zj = data.z[i];
		1250	dataj.cabj = data.cab[i];
		1251	dataj.sabj = data.sab[i];
		1252	dataj.salpj = data.salp[i];
		1253
		1254	if (dataj.iexk != 0)
		1255	{
		1256	ipr = data.icon1[i];
		1257
		1258	if (ipr > PCHCON)
		1259	dataj.ind1 = 2;
		1260	else if (ipr < 0)
		1261	{
		1262	ipr = -ipr;
		1263	iprx = ipr - 1;
		1264
		1265	if (-data.icon1[iprx] != ix)
		1266	dataj.ind1 = 2;
		1267	else
		1268	{
		1269	xi = fabsl (
		1270	dataj.cabj * data.cab[iprx] +
		1271	dataj.sabj * data.sab[iprx] +
		1272	dataj.salpj * data.salp[iprx]);
		1273	if ((xi < 0.999999) \|\|
		1274	(fabsl (data.bi[iprx] / dataj.b - 1.) >
		1275	1.0e-6))
		1276	dataj.ind1 = 2;
		1277	else
		1278	dataj.ind1 = 0;
		1279	}
		1280	} /* if( ipr < 0 ) */
		1281	else if (ipr == 0)
		1282	dataj.ind1 = 1;
		1283	else
		1284	{
		1285	iprx = ipr - 1;
		1286
		1287	if (ipr != ix)
		1288	{
		1289	if (data.icon2[iprx] != ix)
		1290	dataj.ind1 = 2;
		1291	else
		1292	{
		1293	xi = fabsl (
		1294	dataj.cabj * data.cab[iprx] +
		1295	dataj.sabj * data.sab[iprx] +
		1296	dataj.salpj * data.salp[iprx]);
		1297	if ((xi < 0.999999) \|\|
		1298	(fabsl (
		1299	data.bi[iprx] / dataj.b -
		1300	1.) > 1.0e-6))
		1301	dataj.ind1 = 2;
		1302	else
		1303	dataj.ind1 = 0;
		1304	}
		1305	} /* if( ipr != ix ) */
		1306	else
		1307	{
		1308	if (dataj.cabj * dataj.cabj +
		1309	dataj.sabj * dataj.sabj >
		1310	1.0e-8)
		1311	dataj.ind1 = 2;
		1312	else
		1313	dataj.ind1 = 0;
		1314	}
		1315	} /* else */
		1316
		1317	ipr = data.icon2[i];
		1318
		1319	if (ipr > PCHCON)
		1320	dataj.ind2 = 2;
		1321	else if (ipr < 0)
		1322	{
		1323	ipr = -ipr;
		1324	iprx = ipr - 1;
		1325
		1326	if (-data.icon2[iprx] != ix)
		1327	dataj.ind1 = 2;
		1328	else
		1329	{
		1330	xi = fabsl (
		1331	dataj.cabj * data.cab[iprx] +
		1332	dataj.sabj * data.sab[iprx] +
		1333	dataj.salpj * data.salp[iprx]);
		1334	if ((xi < 0.999999) \|\|
		1335	(fabsl (data.bi[iprx] / dataj.b - 1.) >
		1336	1.0e-6))
		1337	dataj.ind1 = 2;
		1338	else
		1339	dataj.ind1 = 0;
		1340	}
		1341	} /* if( ipr < 0 ) */
		1342	else if (ipr == 0)
		1343	dataj.ind2 = 1;
		1344	else
		1345	{
		1346	iprx = ipr - 1;
		1347
		1348	if (ipr != ix)
		1349	{
		1350	if (data.icon1[iprx] != ix)
		1351	dataj.ind2 = 2;
		1352	else
		1353	{
		1354	xi = fabsl (
		1355	dataj.cabj * data.cab[iprx] +
		1356	dataj.sabj * data.sab[iprx] +
		1357	dataj.salpj * data.salp[iprx]);
		1358	if ((xi < 0.999999) \|\|
		1359	(fabsl (
		1360	data.bi[iprx] / dataj.b -
		1361	1.) > 1.0e-6))
		1362	dataj.ind2 = 2;
		1363	else
		1364	dataj.ind2 = 0;
		1365	}
		1366	} /* if( ipr != (i+1) ) */
		1367	else
		1368	{
		1369	if (dataj.cabj * dataj.cabj +
		1370	dataj.sabj * dataj.sabj >
		1371	1.0e-8)
		1372	dataj.ind1 = 2;
		1373	else
		1374	dataj.ind1 = 0;
		1375	}
		1376
		1377	} /* else */
		1378
		1379	} /* if( dataj.iexk != 0) */
		1380
		1381	efld (xob, yob, zob, ax, 1);
		1382	acx = cmplx (crnt.air[i], crnt.aii[i]);
		1383	bcx = cmplx (crnt.bir[i], crnt.bii[i]);
		1384	ccx = cmplx (crnt.cir[i], crnt.cii[i]);
		1385	ex += dataj.exk acx + dataj.exs * bcx + dataj.exc * ccx;
		1386	ey += dataj.eyk acx + dataj.eys * bcx + dataj.eyc * ccx;
		1387	ez += dataj.ezk acx + dataj.ezs * bcx + dataj.ezc * ccx;
		1388
		1389	} /* for( i = 0; i < n; i++ ) */
		1390
		1391	if (data.m == 0)
		1392	return;
		1393
		1394	} /* if( n != 0) */
		1395
		1396	jc = data.n - 1;
		1397	for (i = 0; i < data.m; i++)
		1398	{
		1399	dataj.s = data.pbi[i];
		1400	dataj.xj = data.px[i];
		1401	dataj.yj = data.py[i];
		1402	dataj.zj = data.pz[i];
		1403	dataj.t1xj = data.t1x[i];
		1404	dataj.t1yj = data.t1y[i];
		1405	dataj.t1zj = data.t1z[i];
		1406	dataj.t2xj = data.t2x[i];
		1407	dataj.t2yj = data.t2y[i];
		1408	dataj.t2zj = data.t2z[i];
		1409	jc += 3;
		1410	acx = dataj.t1xj * crnt.cur[jc - 2] + dataj.t1yj * crnt.cur[jc - 1] +
		1411	dataj.t1zj * crnt.cur[jc];
		1412	bcx = dataj.t2xj * crnt.cur[jc - 2] + dataj.t2yj * crnt.cur[jc - 1] +
		1413	dataj.t2zj * crnt.cur[jc];
		1414
		1415	for (ipa = 0; ipa < gnd.ksymp; ipa++)
		1416	{
		1417	dataj.ipgnd = ipa + 1;
		1418	unere (xob, yob, zob);
		1419	ex = ex + acx * dataj.exk + bcx * dataj.exs;
		1420	ey = ey + acx * dataj.eyk + bcx * dataj.eys;
		1421	ez = ez + acx * dataj.ezk + bcx * dataj.ezs;
		1422	}
		1423
		1424	} /* for( i = 0; i < m; i++ ) */
		1425
		1426	return;
		1427	}
		1428
		1429	/-----------------------------------------------------------------------/
		1430
		1431	/* compute near e or h fields over a range of points */
		1432	void nfpat (void)
		1433	{
		1434	int i, j, kk;
		1435	double znrt, cth = 0., sth = 0., ynrt, cph = 0., sph = 0., xnrt, xob, yob;
		1436	double zob, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, xxx;
		1437	complex double ex, ey, ez;
		1438
		1439	if (fpat.nfeh != 1)
		1440	{
		1441	fprintf (
		1442	output_fp,
		1443	"\n\n\n"
		1444	" "
		1445	"-------- NEAR ELECTRIC FIELDS --------\n"
		1446	" ------- LOCATION ------- ------- EX ------ ------- EY ------ "
		1447	" ------- EZ ------\n"
		1448	" X Y Z MAGNITUDE PHASE MAGNITUDE PHASE "
		1449	" MAGNITUDE PHASE\n"
		1450	" METERS METERS METERS VOLTS/M DEGREES VOLTS/M DEGREES "
		1451	" VOLTS/M DEGREES");
		1452	}
		1453	else
		1454	{
		1455	fprintf (
		1456	output_fp,
		1457	"\n\n\n"
		1458	" "
		1459	"-------- NEAR MAGNETIC FIELDS ---------\n\n"
		1460	" ------- LOCATION ------- ------- HX ------ ------- HY ------ "
		1461	" ------- HZ ------\n"
		1462	" X Y Z MAGNITUDE PHASE MAGNITUDE PHASE "
		1463	" MAGNITUDE PHASE\n"
		1464	" METERS METERS METERS AMPS/M DEGREES AMPS/M DEGREES "
		1465	" AMPS/M DEGREES");
		1466	}
		1467
		1468	znrt = fpat.znr - fpat.dznr;
		1469	for (i = 0; i < fpat.nrz; i++)
		1470	{
		1471	znrt += fpat.dznr;
		1472	if (fpat.near != 0)
		1473	{
		1474	cth = cosl (TA * znrt);
		1475	sth = sinl (TA * znrt);
		1476	}
		1477
		1478	ynrt = fpat.ynr - fpat.dynr;
		1479	for (j = 0; j < fpat.nry; j++)
		1480	{
		1481	ynrt += fpat.dynr;
		1482	if (fpat.near != 0)
		1483	{
		1484	cph = cosl (TA * ynrt);
		1485	sph = sinl (TA * ynrt);
		1486	}
		1487
		1488	xnrt = fpat.xnr - fpat.dxnr;
		1489	for (kk = 0; kk < fpat.nrx; kk++)
		1490	{
		1491	xnrt += fpat.dxnr;
		1492	if (fpat.near != 0)
		1493	{
		1494	xob = xnrt * sth * cph;
		1495	yob = xnrt * sth * sph;
		1496	zob = xnrt * cth;
		1497	}
		1498	else
		1499	{
		1500	xob = xnrt;
		1501	yob = ynrt;
		1502	zob = znrt;
		1503	}
		1504
		1505	tmp1 = xob / data.wlam;
		1506	tmp2 = yob / data.wlam;
		1507	tmp3 = zob / data.wlam;
		1508
		1509	if (fpat.nfeh != 1)
		1510	nefld (tmp1, tmp2, tmp3, &ex, &ey, &ez);
		1511	else
		1512	nhfld (tmp1, tmp2, tmp3, &ex, &ey, &ez);
		1513
		1514	tmp1 = cabsl (ex);
		1515	tmp2 = cang (ex);
		1516	tmp3 = cabsl (ey);
		1517	tmp4 = cang (ey);
		1518	tmp5 = cabsl (ez);
		1519	tmp6 = cang (ez);
		1520
		1521	fprintf (
		1522	output_fp,
		1523	"\n"
		1524	" %9.4f %9.4f %9.4f %11.4E %7.2f %11.4E %7.2f %11.4E "
		1525	"%7.2f",
		1526	xob,
		1527	yob,
		1528	zob,
		1529	tmp1,
		1530	tmp2,
		1531	tmp3,
		1532	tmp4,
		1533	tmp5,
		1534	tmp6);
		1535
		1536	if (plot.iplp1 != 2)
		1537	continue;
		1538
		1539	if (plot.iplp4 < 0)
		1540	xxx = xob;
		1541	else if (plot.iplp4 == 0)
		1542	xxx = yob;
		1543	else
		1544	xxx = zob;
		1545
		1546	if (plot.iplp2 == 2)
		1547	{
		1548	switch (plot.iplp3)
		1549	{
		1550	case 1:
		1551	fprintf (
		1552	plot_fp,
		1553	"%12.4E %12.4E %12.4E\n",
		1554	xxx,
		1555	tmp1,
		1556	tmp2);
		1557	break;
		1558	case 2:
		1559	fprintf (
		1560	plot_fp,
		1561	"%12.4E %12.4E %12.4E\n",
		1562	xxx,
		1563	tmp3,
		1564	tmp4);
		1565	break;
		1566	case 3:
		1567	fprintf (
		1568	plot_fp,
		1569	"%12.4E %12.4E %12.4E\n",
		1570	xxx,
		1571	tmp5,
		1572	tmp6);
		1573	break;
		1574	case 4:
		1575	fprintf (
		1576	plot_fp,
		1577	"%12.4E %12.4E %12.4E %12.4E %12.4E "
		1578	"%12.4E %12.4E\n",
		1579	xxx,
		1580	tmp1,
		1581	tmp2,
		1582	tmp3,
		1583	tmp4,
		1584	tmp5,
		1585	tmp6);
		1586	}
		1587	continue;
		1588	}
		1589
		1590	if (plot.iplp2 != 1)
		1591	continue;
		1592
		1593	switch (plot.iplp3)
		1594	{
		1595	case 1:
		1596	fprintf (
		1597	plot_fp,
		1598	"%12.4E %12.4E %12.4E\n",
		1599	xxx,
		1600	print_creall (ex),
		1601	print_cimagl (ex));
		1602	break;
		1603	case 2:
		1604	fprintf (
		1605	plot_fp,
		1606	"%12.4E %12.4E %12.4E\n",
		1607	xxx,
		1608	print_creall (ey),
		1609	print_cimagl (ey));
		1610	break;
		1611	case 3:
		1612	fprintf (
		1613	plot_fp,
		1614	"%12.4E %12.4E %12.4E\n",
		1615	xxx,
		1616	print_creall (ez),
		1617	print_cimagl (ez));
		1618	break;
		1619	case 4:
		1620	fprintf (
		1621	plot_fp,
		1622	"%12.4E %12.4E %12.4E %12.4E %12.4E %12.4E "
		1623	"%12.4E\n",
		1624	xxx,
		1625	print_creall (ex),
		1626	print_cimagl (ex),
		1627	print_creall (ey),
		1628	print_cimagl (ey),
		1629	print_creall (ez),
		1630	print_cimagl (ez));
		1631	}
		1632	} /* for( kk = 0; kk < fpat.nrx; kk++ ) */
		1633
		1634	} /* for( j = 0; j < fpat.nry; j++ ) */
		1635
		1636	} /* for( i = 0; i < fpat.nrz; i++ ) */
		1637
		1638	return;
		1639	}
		1640
		1641	/-----------------------------------------------------------------------/
		1642
		1643	/* nhfld computes the near field at specified points in space after */
		1644	/* the structure currents have been computed. */
		1645
		1646	void nhfld (
		1647	double xob,
		1648	double yob,
		1649	double zob,
		1650	complex double *hx,
		1651	complex double *hy,
		1652	complex double *hz)
		1653	{
		1654	int i, jc;
		1655	double ax, zp;
		1656	complex double acx, bcx, ccx;
		1657
		1658	*hx = CPLX_00;
		1659	*hy = CPLX_00;
		1660	*hz = CPLX_00;
		1661	ax = 0.;
		1662
		1663	if (data.n != 0)
		1664	{
		1665	for (i = 0; i < data.n; i++)
		1666	{
		1667	dataj.xj = xob - data.x[i];
		1668	dataj.yj = yob - data.y[i];
		1669	dataj.zj = zob - data.z[i];
		1670	zp = data.cab[i] * dataj.xj + data.sab[i] * dataj.yj +
		1671	data.salp[i] * dataj.zj;
		1672
		1673	if (fabsl (zp) > 0.5001 * data.si[i])
		1674	continue;
		1675
		1676	zp = dataj.xj * dataj.xj + dataj.yj * dataj.yj + dataj.zj * dataj.zj -
		1677	zp * zp;
		1678	dataj.xj = data.bi[i];
		1679
		1680	if (zp > 0.9 * dataj.xj * dataj.xj)
		1681	continue;
		1682
		1683	ax = dataj.xj;
		1684	break;
		1685	}
		1686
		1687	for (i = 0; i < data.n; i++)
		1688	{
		1689	dataj.s = data.si[i];
		1690	dataj.b = data.bi[i];
		1691	dataj.xj = data.x[i];
		1692	dataj.yj = data.y[i];
		1693	dataj.zj = data.z[i];
		1694	dataj.cabj = data.cab[i];
		1695	dataj.sabj = data.sab[i];
		1696	dataj.salpj = data.salp[i];
		1697	hsfld (xob, yob, zob, ax);
		1698	acx = cmplx (crnt.air[i], crnt.aii[i]);
		1699	bcx = cmplx (crnt.bir[i], crnt.bii[i]);
		1700	ccx = cmplx (crnt.cir[i], crnt.cii[i]);
		1701	hx += dataj.exk acx + dataj.exs * bcx + dataj.exc * ccx;
		1702	hy += dataj.eyk acx + dataj.eys * bcx + dataj.eyc * ccx;
		1703	hz += dataj.ezk acx + dataj.ezs * bcx + dataj.ezc * ccx;
		1704	}
		1705
		1706	if (data.m == 0)
		1707	return;
		1708
		1709	} /* if( data.n != 0) */
		1710
		1711	jc = data.n - 1;
		1712	for (i = 0; i < data.m; i++)
		1713	{
		1714	dataj.s = data.pbi[i];
		1715	dataj.xj = data.px[i];
		1716	dataj.yj = data.py[i];
		1717	dataj.zj = data.pz[i];
		1718	dataj.t1xj = data.t1x[i];
		1719	dataj.t1yj = data.t1y[i];
		1720	dataj.t1zj = data.t1z[i];
		1721	dataj.t2xj = data.t2x[i];
		1722	dataj.t2yj = data.t2y[i];
		1723	dataj.t2zj = data.t2z[i];
		1724	hintg (xob, yob, zob);
		1725	jc += 3;
		1726	acx = dataj.t1xj * crnt.cur[jc - 2] + dataj.t1yj * crnt.cur[jc - 1] +
		1727	dataj.t1zj * crnt.cur[jc];
		1728	bcx = dataj.t2xj * crnt.cur[jc - 2] + dataj.t2yj * crnt.cur[jc - 1] +
		1729	dataj.t2zj * crnt.cur[jc];
		1730	hx = hx + acx * dataj.exk + bcx * dataj.exs;
		1731	hy = hy + acx * dataj.eyk + bcx * dataj.eys;
		1732	hz = hz + acx * dataj.ezk + bcx * dataj.ezs;
		1733	}
		1734
		1735	return;
		1736	}
		1737
		1738	/-----------------------------------------------------------------------/
		1739
		1740	/* integrate over patches at wire connection point */
		1741	void pcint (
		1742	double xi, double yi, double zi, double cabi, double sabi, double salpi, complex double *e)
		1743	{
		1744	int nint, i1, i2;
		1745	double d, ds, da, gcon, fcon, xxj, xyj, xzj, xs, s1;
		1746	double xss, yss, zss, s2x, s2, g1, g2, g3, g4, f2, f1;
		1747	complex double e1, e2, e3, e4, e5, e6, e7, e8, e9;
		1748
		1749	nint = 10;
		1750	d = sqrtl (dataj.s) * .5;
		1751	ds = 4. * d / (double) nint;
		1752	da = ds * ds;
		1753	gcon = 1. / dataj.s;
		1754	fcon = 1. / (2. * TP * d);
		1755	xxj = dataj.xj;
		1756	xyj = dataj.yj;
		1757	xzj = dataj.zj;
		1758	xs = dataj.s;
		1759	dataj.s = da;
		1760	s1 = d + ds * .5;
		1761	xss = dataj.xj + s1 * (dataj.t1xj + dataj.t2xj);
		1762	yss = dataj.yj + s1 * (dataj.t1yj + dataj.t2yj);
		1763	zss = dataj.zj + s1 * (dataj.t1zj + dataj.t2zj);
		1764	s1 = s1 + d;
		1765	s2x = s1;
		1766	e1 = CPLX_00;
		1767	e2 = CPLX_00;
		1768	e3 = CPLX_00;
		1769	e4 = CPLX_00;
		1770	e5 = CPLX_00;
		1771	e6 = CPLX_00;
		1772	e7 = CPLX_00;
		1773	e8 = CPLX_00;
		1774	e9 = CPLX_00;
		1775
		1776	for (i1 = 0; i1 < nint; i1++)
		1777	{
		1778	s1 = s1 - ds;
		1779	s2 = s2x;
		1780	xss = xss - ds * dataj.t1xj;
		1781	yss = yss - ds * dataj.t1yj;
		1782	zss = zss - ds * dataj.t1zj;
		1783	dataj.xj = xss;
		1784	dataj.yj = yss;
		1785	dataj.zj = zss;
		1786
		1787	for (i2 = 0; i2 < nint; i2++)
		1788	{
		1789	s2 = s2 - ds;
		1790	dataj.xj = dataj.xj - ds * dataj.t2xj;
		1791	dataj.yj = dataj.yj - ds * dataj.t2yj;
		1792	dataj.zj = dataj.zj - ds * dataj.t2zj;
		1793	unere (xi, yi, zi);
		1794	dataj.exk = dataj.exk * cabi + dataj.eyk * sabi + dataj.ezk * salpi;
		1795	dataj.exs = dataj.exs * cabi + dataj.eys * sabi + dataj.ezs * salpi;
		1796	g1 = (d + s1) * (d + s2) * gcon;
		1797	g2 = (d - s1) * (d + s2) * gcon;
		1798	g3 = (d - s1) * (d - s2) * gcon;
		1799	g4 = (d + s1) * (d - s2) * gcon;
		1800	f2 = (s1 * s1 + s2 * s2) * TP;
		1801	f1 = s1 / f2 - (g1 - g2 - g3 + g4) * fcon;
		1802	f2 = s2 / f2 - (g1 + g2 - g3 - g4) * fcon;
		1803	e1 = e1 + dataj.exk * g1;
		1804	e2 = e2 + dataj.exk * g2;
		1805	e3 = e3 + dataj.exk * g3;
		1806	e4 = e4 + dataj.exk * g4;
		1807	e5 = e5 + dataj.exs * g1;
		1808	e6 = e6 + dataj.exs * g2;
		1809	e7 = e7 + dataj.exs * g3;
		1810	e8 = e8 + dataj.exs * g4;
		1811	e9 = e9 + dataj.exk * f1 + dataj.exs * f2;
		1812
		1813	} /* for( i2 = 0; i2 < nint; i2++ ) */
		1814
		1815	} /* for( i1 = 0; i1 < nint; i1++ ) */
		1816
		1817	e[0] = e1;
		1818	e[1] = e2;
		1819	e[2] = e3;
		1820	e[3] = e4;
		1821	e[4] = e5;
		1822	e[5] = e6;
		1823	e[6] = e7;
		1824	e[7] = e8;
		1825	e[8] = e9;
		1826	dataj.xj = xxj;
		1827	dataj.yj = xyj;
		1828	dataj.zj = xzj;
		1829	dataj.s = xs;
		1830
		1831	return;
		1832	}
		1833
		1834	/-----------------------------------------------------------------------/
		1835
		1836	/* calculates the electric field due to unit current */
		1837	/* in the t1 and t2 directions on a patch */
		1838	void unere (double xob, double yob, double zob)
		1839	{
		1840	double zr, t1zr, t2zr, rx, ry, rz, r, tt1;
		1841	double tt2, rt, xymag, px, py, cth, r2;
		1842	complex double er, q1, q2, rrv, rrh, edp;
		1843
		1844	zr = dataj.zj;
		1845	t1zr = dataj.t1zj;
		1846	t2zr = dataj.t2zj;
		1847
		1848	if (dataj.ipgnd == 2)
		1849	{
		1850	zr = -zr;
		1851	t1zr = -t1zr;
		1852	t2zr = -t2zr;
		1853	}
		1854
		1855	rx = xob - dataj.xj;
		1856	ry = yob - dataj.yj;
		1857	rz = zob - zr;
		1858	r2 = rx * rx + ry * ry + rz * rz;
		1859
		1860	if (r2 <= 1.0e-20)
		1861	{
		1862	dataj.exk = CPLX_00;
		1863	dataj.eyk = CPLX_00;
		1864	dataj.ezk = CPLX_00;
		1865	dataj.exs = CPLX_00;
		1866	dataj.eys = CPLX_00;
		1867	dataj.ezs = CPLX_00;
		1868	return;
		1869	}
		1870
		1871	r = sqrtl (r2);
		1872	tt1 = -TP * r;
		1873	tt2 = tt1 * tt1;
		1874	rt = r2 * r;
		1875	er = cmplx (sinl (tt1), -cosl (tt1)) * (CONST2 * dataj.s);
		1876	q1 = cmplx (tt2 - 1., tt1) * er / rt;
		1877	q2 = cmplx (3. - tt2, -3. * tt1) * er / (rt * r2);
		1878	er = q2 * (dataj.t1xj * rx + dataj.t1yj * ry + t1zr * rz);
		1879	dataj.exk = q1 * dataj.t1xj + er * rx;
		1880	dataj.eyk = q1 * dataj.t1yj + er * ry;
		1881	dataj.ezk = q1 * t1zr + er * rz;
		1882	er = q2 * (dataj.t2xj * rx + dataj.t2yj * ry + t2zr * rz);
		1883	dataj.exs = q1 * dataj.t2xj + er * rx;
		1884	dataj.eys = q1 * dataj.t2yj + er * ry;
		1885	dataj.ezs = q1 * t2zr + er * rz;
		1886
		1887	if (dataj.ipgnd == 1)
		1888	return;
		1889
		1890	if (gnd.iperf == 1)
		1891	{
		1892	dataj.exk = -dataj.exk;
		1893	dataj.eyk = -dataj.eyk;
		1894	dataj.ezk = -dataj.ezk;
		1895	dataj.exs = -dataj.exs;
		1896	dataj.eys = -dataj.eys;
		1897	dataj.ezs = -dataj.ezs;
		1898	return;
		1899	}
		1900
		1901	xymag = sqrtl (rx * rx + ry * ry);
		1902	if (xymag <= 1.0e-6)
		1903	{
		1904	px = 0.;
		1905	py = 0.;
		1906	cth = 1.;
		1907	rrv = CPLX_10;
		1908	}
		1909	else
		1910	{
		1911	px = -ry / xymag;
		1912	py = rx / xymag;
		1913	cth = rz / sqrtl (xymag * xymag + rz * rz);
		1914	rrv = csqrtl (1. - gnd.zrati * gnd.zrati * (1. - cth * cth));
		1915	}
		1916
		1917	rrh = gnd.zrati * cth;
		1918	rrh = (rrh - rrv) / (rrh + rrv);
		1919	rrv = gnd.zrati * rrv;
		1920	rrv = -(cth - rrv) / (cth + rrv);
		1921	edp = (dataj.exk * px + dataj.eyk * py) * (rrh - rrv);
		1922	dataj.exk = dataj.exk * rrv + edp * px;
		1923	dataj.eyk = dataj.eyk * rrv + edp * py;
		1924	dataj.ezk = dataj.ezk * rrv;
		1925	edp = (dataj.exs * px + dataj.eys * py) * (rrh - rrv);
		1926	dataj.exs = dataj.exs * rrv + edp * px;
		1927	dataj.eys = dataj.eys * rrv + edp * py;
		1928	dataj.ezs = dataj.ezs * rrv;
		1929
		1930	return;
		1931	}
		1932
		1933	/-----------------------------------------------------------------------/

Subversion Repositories Nec2c

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