WebSVN – Nec2c – Blame – /trunk/src/calculations.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 /tmi/ */
		29	extern tmi_t tmi;
		30
		31	/common /ggrid/ /
		32	extern ggrid_t ggrid;
		33
		34	/* common /data/ */
		35	extern data_t data;
		36
		37	/* common /crnt/ */
		38	extern crnt_t crnt;
		39
		40	/* common /vsorc/ */
		41	extern vsorc_t vsorc;
		42
		43	/* common /segj/ */
		44	extern segj_t segj;
		45
		46	/* common /yparm/ */
		47	extern yparm_t yparm;
		48
		49	/* common /zload/ */
		50	extern zload_t zload;
		51
		52	/* common /smat/ */
		53	extern smat_t smat;
		54
		55	/* pointers to input/output files */
		56	extern FILE input_fp, output_fp, *plot_fp;
		57
		58	/-----------------------------------------------------------------------/
		59
		60	/* cabc computes coefficients of the constant (a), sine (b), and */
		61	/* cosine (c) terms in the current interpolation functions for the */
		62	/* current vector cur. */
		63	void cabc (complex double *curx)
		64	{
		65	int i, is, j, jx, jco1, jco2;
		66	double ar, ai, sh;
		67	complex double curd, cs1, cs2;
		68
		69	if (data.n != 0)
		70	{
		71	for (i = 0; i < data.n; i++)
		72	{
		73	crnt.air[i] = 0.;
		74	crnt.aii[i] = 0.;
		75	crnt.bir[i] = 0.;
		76	crnt.bii[i] = 0.;
		77	crnt.cir[i] = 0.;
		78	crnt.cii[i] = 0.;
		79	}
		80
		81	for (i = 0; i < data.n; i++)
		82	{
		83	ar = creall (curx[i]);
		84	ai = cimagl (curx[i]);
		85	tbf (i + 1, 1);
		86
		87	for (jx = 0; jx < segj.jsno; jx++)
		88	{
		89	j = segj.jco[jx] - 1;
		90	crnt.air[j] += segj.ax[jx] * ar;
		91	crnt.aii[j] += segj.ax[jx] * ai;
		92	crnt.bir[j] += segj.bx[jx] * ar;
		93	crnt.bii[j] += segj.bx[jx] * ai;
		94	crnt.cir[j] += segj.cx[jx] * ar;
		95	crnt.cii[j] += segj.cx[jx] * ai;
		96	}
		97
		98	} /* for( i = 0; i < n; i++ ) */
		99
		100	if (vsorc.nqds != 0)
		101	{
		102	for (is = 0; is < vsorc.nqds; is++)
		103	{
		104	i = vsorc.iqds[is] - 1;
		105	jx = data.icon1[i];
		106	data.icon1[i] = 0;
		107	tbf (i + 1, 0);
		108	data.icon1[i] = jx;
		109	sh = data.si[i] * .5;
		110	curd = CCJ * vsorc.vqds[is] /
		111	((logl (2. * sh / data.bi[i]) - 1.) *
		112	(segj.bx[segj.jsno - 1] * cosl (TP * sh) +
		113	segj.cx[segj.jsno - 1] * sinl (TP * sh)) *
		114	data.wlam);
		115	ar = print_creall (curd);
		116	ai = print_cimagl (curd);
		117
		118	for (jx = 0; jx < segj.jsno; jx++)
		119	{
		120	j = segj.jco[jx] - 1;
		121	crnt.air[j] = crnt.air[j] + segj.ax[jx] * ar;
		122	crnt.aii[j] = crnt.aii[j] + segj.ax[jx] * ai;
		123	crnt.bir[j] = crnt.bir[j] + segj.bx[jx] * ar;
		124	crnt.bii[j] = crnt.bii[j] + segj.bx[jx] * ai;
		125	crnt.cir[j] = crnt.cir[j] + segj.cx[jx] * ar;
		126	crnt.cii[j] = crnt.cii[j] + segj.cx[jx] * ai;
		127	}
		128
		129	} /* for( is = 0; is < vsorc.nqds; is++ ) */
		130
		131	} /* if( vsorc.nqds != 0) */
		132
		133	for (i = 0; i < data.n; i++)
		134	curx[i] = cmplx (crnt.air[i] + crnt.cir[i], crnt.aii[i] + crnt.cii[i]);
		135
		136	} /* if( n != 0) */
		137
		138	if (data.m == 0)
		139	return;
		140
		141	/* convert surface currents from */
		142	/* t1,t2 components to x,y,z components */
		143	jco1 = data.np2m;
		144	jco2 = jco1 + data.m;
		145	for (i = 1; i <= data.m; i++)
		146	{
		147	jco1 -= 2;
		148	jco2 -= 3;
		149	cs1 = curx[jco1];
		150	cs2 = curx[jco1 + 1];
		151	curx[jco2] = cs1 * data.t1x[data.m - i] + cs2 * data.t2x[data.m - i];
		152	curx[jco2 + 1] = cs1 * data.t1y[data.m - i] + cs2 * data.t2y[data.m - i];
		153	curx[jco2 + 2] = cs1 * data.t1z[data.m - i] + cs2 * data.t2z[data.m - i];
		154	}
		155
		156	return;
		157	}
		158
		159	/-----------------------------------------------------------------------/
		160
		161	/* couple computes the maximum coupling between pairs of segments. */
		162	void couple (complex double *cur, double wlam)
		163	{
		164	int j, j1, j2, l1, i, k, itt1, itt2, its1, its2, isg1, isg2, npm1;
		165	double dbc, c, gmax;
		166	complex double y11, y12, y22, yl, yin, zl, zin, rho;
		167
		168	if ((vsorc.nsant != 1) \|\| (vsorc.nvqd != 0))
		169	return;
		170
		171	j = isegno (yparm.nctag[yparm.icoup], yparm.ncseg[yparm.icoup]);
		172	if (j != vsorc.isant[0])
		173	return;
		174
		175	zin = vsorc.vsant[0];
		176	yparm.icoup++;
		177	mem_realloc ((void ) &yparm.y11a, yparm.icoup sizeof (complex double));
		178	yparm.y11a[yparm.icoup - 1] = cur[j - 1] * wlam / zin;
		179
		180	l1 = (yparm.icoup - 1) * (yparm.ncoup - 1);
		181	for (i = 0; i < yparm.ncoup; i++)
		182	{
		183	if ((i + 1) == yparm.icoup)
		184	continue;
		185
		186	l1++;
		187	mem_realloc ((void ) &yparm.y12a, l1 sizeof (complex double));
		188	k = isegno (yparm.nctag[i], yparm.ncseg[i]);
		189	yparm.y12a[l1 - 1] = cur[k - 1] * wlam / zin;
		190	}
		191
		192	if (yparm.icoup < yparm.ncoup)
		193	return;
		194
		195	fprintf (
		196	output_fp,
		197	"\n\n\n"
		198	" -----------"
		199	" ISOLATION DATA -----------\n\n"
		200	" ------- COUPLING BETWEEN ------ MAXIMUM "
		201	" ---------- FOR MAXIMUM COUPLING ----------\n"
		202	" SEG SEG COUPLING LOAD"
		203	" IMPEDANCE (2ND SEG) INPUT IMPEDANCE \n"
		204	" TAG SEG NO. TAG SEG NO. (DB) "
		205	" REAL IMAGINARY REAL IMAGINARY");
		206
		207	npm1 = yparm.ncoup - 1;
		208
		209	for (i = 0; i < npm1; i++)
		210	{
		211	itt1 = yparm.nctag[i];
		212	its1 = yparm.ncseg[i];
		213	isg1 = isegno (itt1, its1);
		214	l1 = i + 1;
		215
		216	for (j = l1; j < yparm.ncoup; j++)
		217	{
		218	itt2 = yparm.nctag[j];
		219	its2 = yparm.ncseg[j];
		220	isg2 = isegno (itt2, its2);
		221	j1 = j + i * npm1 - 1;
		222	j2 = i + j * npm1;
		223	y11 = yparm.y11a[i];
		224	y22 = yparm.y11a[j];
		225	y12 = .5 * (yparm.y12a[j1] + yparm.y12a[j2]);
		226	yin = y12 * y12;
		227	dbc = cabsl (yin);
		228	c = dbc / (2. * creall (y11) * creall (y22) - creall (yin));
		229
		230	if ((c >= 0.0) && (c <= 1.0))
		231	{
		232	if (c >= .01)
		233	gmax = (1. - sqrtl (1. - c * c)) / c;
		234	else
		235	gmax = .5 * (c + .25 * c * c * c);
		236
		237	rho = gmax * conjl (yin) / dbc;
		238	yl = ((1. - rho) / (1. + rho) + 1.) * creall (y22) - y22;
		239	zl = 1. / yl;
		240	yin = y11 - yin / (y22 + yl);
		241	zin = 1. / yin;
		242	dbc = db10 (gmax);
		243
		244	fprintf (
		245	output_fp,
		246	"\n"
		247	" %4d %4d %5d %4d %4d %5d %9.3f"
		248	" %12.5f %12.5f %12.5f %12.5f",
		249	itt1,
		250	its1,
		251	isg1,
		252	itt2,
		253	its2,
		254	isg2,
		255	dbc,
		256	print_creall (zl),
		257	print_cimagl (zl),
		258	print_creall (zin),
		259	print_cimagl (zin));
		260
		261	continue;
		262
		263	} /* if( (c >= 0.0) && (c <= 1.0) ) */
		264
		265	fprintf (
		266	output_fp,
		267	"\n"
		268	" %4d %4d %5d %4d %4d %5d ERROR "
		269	"COUPLING IS NOT BETWEEN 0 AND 1. (= %12.5f)",
		270	itt1,
		271	its1,
		272	isg1,
		273	itt2,
		274	its2,
		275	isg2,
		276	c);
		277
		278	} /* for( j = l1; j < yparm.ncoup; j++ ) */
		279
		280	} /* for( i = 0; i < npm1; i++ ) */
		281
		282	return;
		283	}
		284
		285	/-----------------------------------------------------------------------/
		286
		287	/* load calculates the impedance of specified */
		288	/* segments for various types of loading */
		289	void load (
		290	int ldtyp, int ldtag, int ldtagf, int ldtagt, double zlr, double zli, double *zlc)
		291	{
		292	int i, iwarn, istep, istepx, l1, l2, ldtags, jump, ichk;
		293	complex double zt = CPLX_00, tpcj;
		294
		295	tpcj = (0.0 + 1.8836989e-9j);
		296	fprintf (
		297	output_fp,
		298	"\n"
		299	" LOCATION RESISTANCE INDUCTANCE CAPACITANCE "
		300	" IMPEDANCE (OHMS) CONDUCTIVITY CIRCUIT\n"
		301	" OHMS HENRYS FARADS "
		302	" REAL IMAGINARY MHOS/METER TYPE");
		303
		304	/* initialize d array, used for temporary */
		305	/* storage of loading information. */
		306	mem_realloc ((void ) &zload.zarray, data.npm sizeof (complex double));
		307	for (i = 0; i < data.n; i++)
		308	zload.zarray[i] = CPLX_00;
		309
		310	iwarn = FALSE;
		311	istep = 0;
		312
		313	/* cycle over loading cards */
		314	while (TRUE)
		315	{
		316	istepx = istep;
		317	istep++;
		318
		319	if (istep > zload.nload)
		320	{
		321	if (iwarn == TRUE)
		322	fprintf (
		323	output_fp,
		324	"\n NOTE, SOME OF THE ABOVE SEGMENTS "
		325	"HAVE BEEN LOADED TWICE - IMPEDANCES ADDED");
		326
		327	smat.nop = data.n / data.np;
		328	if (smat.nop == 1)
		329	return;
		330
		331	for (i = 0; i < data.np; i++)
		332	{
		333	zt = zload.zarray[i];
		334	l1 = i;
		335
		336	for (l2 = 1; l2 < smat.nop; l2++)
		337	{
		338	l1 += data.np;
		339	zload.zarray[l1] = zt;
		340	}
		341	}
		342	return;
		343
		344	} /* if( istep > zload.nload) */
		345
		346	if (ldtyp[istepx] > 5)
		347	{
		348	fprintf (
		349	output_fp,
		350	"\n IMPROPER LOAD TYPE CHOSEN,"
		351	" REQUESTED TYPE IS %d",
		352	ldtyp[istepx]);
		353	stop (-1);
		354	}
		355
		356	/* search segments for proper itags */
		357	ldtags = ldtag[istepx];
		358	jump = ldtyp[istepx] + 1;
		359	ichk = 0;
		360	l1 = 1;
		361	l2 = data.n;
		362
		363	if (ldtags == 0)
		364	{
		365	if ((ldtagf[istepx] != 0) \|\| (ldtagt[istepx] != 0))
		366	{
		367	l1 = ldtagf[istepx];
		368	l2 = ldtagt[istepx];
		369
		370	} /* if( (ldtagf[istepx] != 0) \|\| (ldtagt[istepx] != 0) ) */
		371
		372	} /* if( ldtags == 0) */
		373
		374	for (i = l1 - 1; i < l2; i++)
		375	{
		376	if (ldtags != 0)
		377	{
		378	if (ldtags != data.itag[i])
		379	continue;
		380
		381	if (ldtagf[istepx] != 0)
		382	{
		383	ichk++;
		384	if ((ichk < ldtagf[istepx]) \|\| (ichk > ldtagt[istepx]))
		385	continue;
		386	}
		387	else
		388	ichk = 1;
		389
		390	} /* if( ldtags != 0) */
		391	else
		392	ichk = 1;
		393
		394	/* calculation of lamdaimped. per unit length, /
		395	/* jump to appropriate section for loading type */
		396	switch (jump)
		397	{
		398	case 1:
		399	zt = zlr[istepx] / data.si[i] +
		400	tpcj * zli[istepx] / (data.si[i] * data.wlam);
		401	if (fabsl (zlc[istepx]) > 1.0e-20)
		402	zt += data.wlam / (tpcj * data.si[i] * zlc[istepx]);
		403	break;
		404
		405	case 2:
		406	zt = tpcj * data.si[i] * zlc[istepx] / data.wlam;
		407	if (fabsl (zli[istepx]) > 1.0e-20)
		408	zt += data.si[i] * data.wlam / (tpcj * zli[istepx]);
		409	if (fabsl (zlr[istepx]) > 1.0e-20)
		410	zt += data.si[i] / zlr[istepx];
		411	zt = 1. / zt;
		412	break;
		413
		414	case 3:
		415	zt = zlr[istepx] * data.wlam + tpcj * zli[istepx];
		416	if (fabsl (zlc[istepx]) > 1.0e-20)
		417	zt += 1. /
		418	(tpcj * data.si[i] * data.si[i] * zlc[istepx]);
		419	break;
		420
		421	case 4:
		422	zt = tpcj * data.si[i] * data.si[i] * zlc[istepx];
		423	if (fabsl (zli[istepx]) > 1.0e-20)
		424	zt += 1. / (tpcj * zli[istepx]);
		425	if (fabsl (zlr[istepx]) > 1.0e-20)
		426	zt += 1. / (zlr[istepx] * data.wlam);
		427	zt = 1. / zt;
		428	break;
		429
		430	case 5:
		431	zt = cmplx (zlr[istepx], zli[istepx]) / data.si[i];
		432	break;
		433
		434	case 6:
		435	zint (zlr[istepx] * data.wlam, data.bi[i], &zt);
		436
		437	} /* switch( jump ) */
		438
		439	if ((fabsl (creall (zload.zarray[i])) +
		440	fabsl (cimagl (zload.zarray[i]))) > 1.0e-20)
		441	iwarn = TRUE;
		442	zload.zarray[i] += zt;
		443
		444	} /* for( i = l1-1; i < l2; i++ ) */
		445
		446	if (ichk == 0)
		447	{
		448	fprintf (
		449	output_fp,
		450	"\n LOADING DATA CARD ERROR,"
		451	" NO SEGMENT HAS AN ITAG = %d",
		452	ldtags);
		453	stop (-1);
		454	}
		455
		456	/* printing the segment loading data, jump to proper print */
		457	switch (jump)
		458	{
		459	case 1:
		460	prnt (
		461	ldtags,
		462	ldtagf[istepx],
		463	ldtagt[istepx],
		464	zlr[istepx],
		465	zli[istepx],
		466	zlc[istepx],
		467	0.,
		468	0.,
		469	0.,
		470	" SERIES ",
		471	2);
		472	break;
		473
		474	case 2:
		475	prnt (
		476	ldtags,
		477	ldtagf[istepx],
		478	ldtagt[istepx],
		479	zlr[istepx],
		480	zli[istepx],
		481	zlc[istepx],
		482	0.,
		483	0.,
		484	0.,
		485	"PARALLEL",
		486	2);
		487	break;
		488
		489	case 3:
		490	prnt (
		491	ldtags,
		492	ldtagf[istepx],
		493	ldtagt[istepx],
		494	zlr[istepx],
		495	zli[istepx],
		496	zlc[istepx],
		497	0.,
		498	0.,
		499	0.,
		500	"SERIES (PER METER)",
		501	5);
		502	break;
		503
		504	case 4:
		505	prnt (
		506	ldtags,
		507	ldtagf[istepx],
		508	ldtagt[istepx],
		509	zlr[istepx],
		510	zli[istepx],
		511	zlc[istepx],
		512	0.,
		513	0.,
		514	0.,
		515	"PARALLEL (PER METER)",
		516	5);
		517	break;
		518
		519	case 5:
		520	prnt (
		521	ldtags,
		522	ldtagf[istepx],
		523	ldtagt[istepx],
		524	0.,
		525	0.,
		526	0.,
		527	zlr[istepx],
		528	zli[istepx],
		529	0.,
		530	"FIXED IMPEDANCE ",
		531	4);
		532	break;
		533
		534	case 6:
		535	prnt (
		536	ldtags,
		537	ldtagf[istepx],
		538	ldtagt[istepx],
		539	0.,
		540	0.,
		541	0.,
		542	0.,
		543	0.,
		544	zlr[istepx],
		545	" WIRE ",
		546	2);
		547
		548	} /* switch( jump ) */
		549
		550	} /* while( TRUE ) */
		551
		552	return;
		553	}
		554
		555	/-----------------------------------------------------------------------/
		556
		557	/* gf computes the integrand exp(jkr)/(kr) for numerical integration. */
		558	void gf (double zk, double co, double si)
		559	{
		560	double zdk, rk, rks;
		561
		562	zdk = zk - tmi.zpk;
		563	rk = sqrtl (tmi.rkb2 + zdk * zdk);
		564	*si = sinl (rk) / rk;
		565
		566	if (tmi.ij != 0)
		567	{
		568	*co = cosl (rk) / rk;
		569	return;
		570	}
		571
		572	if (rk >= .2)
		573	{
		574	*co = (cosl (rk) - 1.) / rk;
		575	return;
		576	}
		577
		578	rks = rk * rk;
		579	co = ((-1.38888889e-3 rks + 4.16666667e-2) * rks - .5) * rk;
		580
		581	return;
		582	}
		583
		584	/-----------------------------------------------------------------------/
		585
		586	/* function db10 returns db for magnitude (field) */
		587	double db10 (double x)
		588	{
		589	if (x < 1.e-20)
		590	return (-999.99);
		591
		592	return (10. * log10l (x));
		593	}
		594
		595	/-----------------------------------------------------------------------/
		596
		597	/* function db20 returns db for mag*2 (power) i /
		598	double db20 (double x)
		599	{
		600	if (x < 1.e-20)
		601	return (-999.99);
		602
		603	return (20. * log10l (x));
		604	}
		605
		606	/-----------------------------------------------------------------------/
		607
		608	/* intrp uses bivariate cubic interpolation to obtain */
		609	/* the values of 4 functions at the point (x,y). */
		610	void intrp (
		611	double x,
		612	double y,
		613	complex double *f1,
		614	complex double *f2,
		615	complex double *f3,
		616	complex double *f4)
		617	{
		618	static int ix, iy, ixs = -10, iys = -10, igrs = -10, ixeg = 0, iyeg = 0;
		619	static int nxm2, nym2, nxms, nyms, nd, ndp;
		620	int nda[3] = {11, 17, 9}, ndpa[3] = {110, 85, 72};
		621	int igr, iadd, iadz, i, k, jump;
		622	static double dx = 1., dy = 1., xs = 0., ys = 0., xz, yz;
		623	double xx, yy;
		624	static complex double a[4][4], b[4][4], c[4][4], d[4][4];
		625	complex double p1 = CPLX_00, p2 = CPLX_00, p3 = CPLX_00, p4 = CPLX_00;
		626	complex double fx1, fx2, fx3, fx4;
		627
		628	jump = FALSE;
		629	if ((x < xs) \|\| (y < ys))
		630	jump = TRUE;
		631	else
		632	{
		633	ix = (int) ((x - xs) / dx) + 1;
		634	iy = (int) ((y - ys) / dy) + 1;
		635	}
		636
		637	/* if point lies in same 4 by 4 point region */
		638	/* as previous point, old values are reused. */
		639	if ((ix < ixeg) \|\| (iy < iyeg) \|\| (abs (ix - ixs) >= 2) \|\| (abs (iy - iys) >= 2) \|\|
		640	jump)
		641	{
		642	/* determine correct grid and grid region */
		643	if (x <= ggrid.xsa[1])
		644	igr = 0;
		645	else
		646	{
		647	if (y > ggrid.ysa[2])
		648	igr = 2;
		649	else
		650	igr = 1;
		651	}
		652
		653	if (igr != igrs)
		654	{
		655	igrs = igr;
		656	dx = ggrid.dxa[igrs];
		657	dy = ggrid.dya[igrs];
		658	xs = ggrid.xsa[igrs];
		659	ys = ggrid.ysa[igrs];
		660	nxm2 = ggrid.nxa[igrs] - 2;
		661	nym2 = ggrid.nya[igrs] - 2;
		662	nxms = ((nxm2 + 1) / 3) * 3 + 1;
		663	nyms = ((nym2 + 1) / 3) * 3 + 1;
		664	nd = nda[igrs];
		665	ndp = ndpa[igrs];
		666	ix = (int) ((x - xs) / dx) + 1;
		667	iy = (int) ((y - ys) / dy) + 1;
		668
		669	} /* if( igr != igrs) */
		670
		671	ixs = ((ix - 1) / 3) * 3 + 2;
		672	if (ixs < 2)
		673	ixs = 2;
		674	ixeg = -10000;
		675
		676	if (ixs > nxm2)
		677	{
		678	ixs = nxm2;
		679	ixeg = nxms;
		680	}
		681
		682	iys = ((iy - 1) / 3) * 3 + 2;
		683	if (iys < 2)
		684	iys = 2;
		685	iyeg = -10000;
		686
		687	if (iys > nym2)
		688	{
		689	iys = nym2;
		690	iyeg = nyms;
		691	}
		692
		693	/* compute coefficients of 4 cubic polynomials in x for */
		694	/* the 4 grid values of y for each of the 4 functions */
		695	iadz = ixs + (iys - 3) * nd - ndp;
		696	for (k = 0; k < 4; k++)
		697	{
		698	iadz += ndp;
		699	iadd = iadz;
		700
		701	for (i = 0; i < 4; i++)
		702	{
		703	iadd += nd;
		704
		705	switch (igrs)
		706	{
		707	case 0:
		708	p1 = ggrid.ar1[iadd - 2];
		709	p2 = ggrid.ar1[iadd - 1];
		710	p3 = ggrid.ar1[iadd];
		711	p4 = ggrid.ar1[iadd + 1];
		712	break;
		713
		714	case 1:
		715	p1 = ggrid.ar2[iadd - 2];
		716	p2 = ggrid.ar2[iadd - 1];
		717	p3 = ggrid.ar2[iadd];
		718	p4 = ggrid.ar2[iadd + 1];
		719	break;
		720
		721	case 2:
		722	p1 = ggrid.ar3[iadd - 2];
		723	p2 = ggrid.ar3[iadd - 1];
		724	p3 = ggrid.ar3[iadd];
		725	p4 = ggrid.ar3[iadd + 1];
		726
		727	} /* switch( igrs ) */
		728
		729	a[i][k] = (p4 - p1 + 3. * (p2 - p3)) * .1666666667;
		730	b[i][k] = (p1 - 2. * p2 + p3) * .5;
		731	c[i][k] = p3 - (2. * p1 + 3. * p2 + p4) * .1666666667;
		732	d[i][k] = p2;
		733
		734	} /* for( i = 0; i < 4; i++ ) */
		735
		736	} /* for( k = 0; k < 4; k++ ) */
		737
		738	xz = (ixs - 1) * dx + xs;
		739	yz = (iys - 1) * dy + ys;
		740
		741	} /* if( (abs(ix- ixs) >= 2) \|\| */
		742
		743	/* evaluate polymomials in x and use cubic */
		744	/* interpolation in y for each of the 4 functions. */
		745	xx = (x - xz) / dx;
		746	yy = (y - yz) / dy;
		747	fx1 = ((a[0][0] * xx + b[0][0]) * xx + c[0][0]) * xx + d[0][0];
		748	fx2 = ((a[1][0] * xx + b[1][0]) * xx + c[1][0]) * xx + d[1][0];
		749	fx3 = ((a[2][0] * xx + b[2][0]) * xx + c[2][0]) * xx + d[2][0];
		750	fx4 = ((a[3][0] * xx + b[3][0]) * xx + c[3][0]) * xx + d[3][0];
		751	p1 = fx4 - fx1 + 3. * (fx2 - fx3);
		752	p2 = 3. * (fx1 - 2. * fx2 + fx3);
		753	p3 = 6. * fx3 - 2. * fx1 - 3. * fx2 - fx4;
		754	f1 = ((p1 yy + p2) * yy + p3) * yy * .1666666667 + fx2;
		755	fx1 = ((a[0][1] * xx + b[0][1]) * xx + c[0][1]) * xx + d[0][1];
		756	fx2 = ((a[1][1] * xx + b[1][1]) * xx + c[1][1]) * xx + d[1][1];
		757	fx3 = ((a[2][1] * xx + b[2][1]) * xx + c[2][1]) * xx + d[2][1];
		758	fx4 = ((a[3][1] * xx + b[3][1]) * xx + c[3][1]) * xx + d[3][1];
		759	p1 = fx4 - fx1 + 3. * (fx2 - fx3);
		760	p2 = 3. * (fx1 - 2. * fx2 + fx3);
		761	p3 = 6. * fx3 - 2. * fx1 - 3. * fx2 - fx4;
		762	f2 = ((p1 yy + p2) * yy + p3) * yy * .1666666667 + fx2;
		763	fx1 = ((a[0][2] * xx + b[0][2]) * xx + c[0][2]) * xx + d[0][2];
		764	fx2 = ((a[1][2] * xx + b[1][2]) * xx + c[1][2]) * xx + d[1][2];
		765	fx3 = ((a[2][2] * xx + b[2][2]) * xx + c[2][2]) * xx + d[2][2];
		766	fx4 = ((a[3][2] * xx + b[3][2]) * xx + c[3][2]) * xx + d[3][2];
		767	p1 = fx4 - fx1 + 3. * (fx2 - fx3);
		768	p2 = 3. * (fx1 - 2. * fx2 + fx3);
		769	p3 = 6. * fx3 - 2. * fx1 - 3. * fx2 - fx4;
		770	f3 = ((p1 yy + p2) * yy + p3) * yy * .1666666667 + fx2;
		771	fx1 = ((a[0][3] * xx + b[0][3]) * xx + c[0][3]) * xx + d[0][3];
		772	fx2 = ((a[1][3] * xx + b[1][3]) * xx + c[1][3]) * xx + d[1][3];
		773	fx3 = ((a[2][3] * xx + b[2][3]) * xx + c[2][3]) * xx + d[2][3];
		774	fx4 = ((a[3][3] * xx + b[3][3]) * xx + c[3][3]) * xx + d[3][3];
		775	p1 = fx4 - fx1 + 3. * (fx2 - fx3);
		776	p2 = 3. * (fx1 - 2. * fx2 + fx3);
		777	p3 = 6. * fx3 - 2. * fx1 - 3. * fx2 - fx4;
		778	f4 = ((p1 yy + p2) * yy + p3) * yy * .16666666670 + fx2;
		779
		780	return;
		781	}
		782
		783	/-----------------------------------------------------------------------/
		784
		785	/* intx performs numerical integration of exp(jkr)/r by the method of */
		786	/* variable interval width romberg integration. the integrand value */
		787	/* is supplied by subroutine gf. */
		788	void intx (double el1, double el2, double b, int ij, double sgr, double sgi)
		789	{
		790	int ns, nt;
		791	int nx = 1, nma = 65536, nts = 4;
		792	int flag = TRUE;
		793	double z, s, ze, fnm, ep, zend, fns, dz = 0., zp, dzot = 0., t00r, g1r, g5r, t00i;
		794	double g1i, g5i, t01r, g3r, t01i, g3i, t10r, t10i, te1i, te1r, t02r;
		795	double g2r, g4r, t02i, g2i, g4i, t11r, t11i, t20r, t20i, te2i, te2r;
		796	double rx = 1.0e-4;
		797
		798	z = el1;
		799	ze = el2;
		800	if (ij == 0)
		801	ze = 0.;
		802	s = ze - z;
		803	fnm = nma;
		804	ep = s / (10. * fnm);
		805	zend = ze - ep;
		806	*sgr = 0.;
		807	*sgi = 0.;
		808	ns = nx;
		809	nt = 0;
		810	gf (z, &g1r, &g1i);
		811
		812	while (TRUE)
		813	{
		814	if (flag)
		815	{
		816	fns = ns;
		817	dz = s / fns;
		818	zp = z + dz;
		819
		820	if (zp > ze)
		821	{
		822	dz = ze - z;
		823	if (fabsl (dz) <= ep)
		824	{
		825	/* add contribution of near singularity for diagonal
		826	* term */
		827	if (ij == 0)
		828	{
		829	*sgr =
		830	2. *
		831	(*sgr +
		832	logl ((sqrtl (b * b + s * s) + s) / b));
		833	sgi = 2. *sgi;
		834	}
		835	return;
		836	}
		837
		838	} /* if( zp > ze) */
		839
		840	dzot = dz * .5;
		841	zp = z + dzot;
		842	gf (zp, &g3r, &g3i);
		843	zp = z + dz;
		844	gf (zp, &g5r, &g5i);
		845
		846	} /* if( flag ) */
		847
		848	t00r = (g1r + g5r) * dzot;
		849	t00i = (g1i + g5i) * dzot;
		850	t01r = (t00r + dz * g3r) * 0.5;
		851	t01i = (t00i + dz * g3i) * 0.5;
		852	t10r = (4.0 * t01r - t00r) / 3.0;
		853	t10i = (4.0 * t01i - t00i) / 3.0;
		854
		855	/* test convergence of 3 point romberg result. */
		856	test (t01r, t10r, &te1r, t01i, t10i, &te1i, 0.);
		857	if ((te1i <= rx) && (te1r <= rx))
		858	{
		859	sgr = sgr + t10r;
		860	sgi = sgi + t10i;
		861	nt += 2;
		862
		863	z += dz;
		864	if (z >= zend)
		865	{
		866	/* add contribution of near singularity for diagonal term */
		867	if (ij == 0)
		868	{
		869	sgr = 2.
		870	(sgr + logl ((sqrtl (b b + s * s) + s) / b));
		871	sgi = 2. *sgi;
		872	}
		873	return;
		874	}
		875
		876	g1r = g5r;
		877	g1i = g5i;
		878	if (nt >= nts)
		879	if (ns > nx)
		880	{
		881	/* Double step size */
		882	ns = ns / 2;
		883	nt = 1;
		884	}
		885	flag = TRUE;
		886	continue;
		887
		888	} /* if( (te1i <= rx) && (te1r <= rx) ) */
		889
		890	zp = z + dz * 0.25;
		891	gf (zp, &g2r, &g2i);
		892	zp = z + dz * 0.75;
		893	gf (zp, &g4r, &g4i);
		894	t02r = (t01r + dzot * (g2r + g4r)) * 0.5;
		895	t02i = (t01i + dzot * (g2i + g4i)) * 0.5;
		896	t11r = (4.0 * t02r - t01r) / 3.0;
		897	t11i = (4.0 * t02i - t01i) / 3.0;
		898	t20r = (16.0 * t11r - t10r) / 15.0;
		899	t20i = (16.0 * t11i - t10i) / 15.0;
		900
		901	/* test convergence of 5 point romberg result. */
		902	test (t11r, t20r, &te2r, t11i, t20i, &te2i, 0.);
		903	if ((te2i > rx) \|\| (te2r > rx))
		904	{
		905	nt = 0;
		906	if (ns >= nma)
		907	fprintf (output_fp, "\n STEP SIZE LIMITED AT Z= %10.5G", z);
		908	else
		909	{
		910	/* halve step size */
		911	ns = ns * 2;
		912	fns = ns;
		913	dz = s / fns;
		914	dzot = dz * 0.5;
		915	g5r = g3r;
		916	g5i = g3i;
		917	g3r = g2r;
		918	g3i = g2i;
		919
		920	flag = FALSE;
		921	continue;
		922	}
		923
		924	} /* if( (te2i > rx) \|\| (te2r > rx) ) */
		925
		926	sgr = sgr + t20r;
		927	sgi = sgi + t20i;
		928	nt++;
		929
		930	z += dz;
		931	if (z >= zend)
		932	{
		933	/* add contribution of near singularity for diagonal term */
		934	if (ij == 0)
		935	{
		936	sgr = 2. (sgr + logl ((sqrtl (b b + s * s) + s) / b));
		937	sgi = 2. *sgi;
		938	}
		939	return;
		940	}
		941
		942	g1r = g5r;
		943	g1i = g5i;
		944	if (nt >= nts)
		945	if (ns > nx)
		946	{
		947	/* Double step size */
		948	ns = ns / 2;
		949	nt = 1;
		950	}
		951	flag = TRUE;
		952
		953	} /* while( TRUE ) */
		954	}
		955
		956	/-----------------------------------------------------------------------/
		957
		958	/* returns smallest of two arguments */
		959	int min (int a, int b)
		960	{
		961	if (a < b)
		962	return (a);
		963	else
		964	return (b);
		965	}
		966
		967	/-----------------------------------------------------------------------/
		968
		969	/* test for convergence in numerical integration */
		970	void test (double f1r, double f2r, double tr, double f1i, double f2i, double ti, double dmin)
		971	{
		972	double den;
		973
		974	den = fabsl (f2r);
		975	*tr = fabsl (f2i);
		976
		977	if (den < *tr)
		978	den = *tr;
		979	if (den < dmin)
		980	den = dmin;
		981
		982	if (den < 1.0e-37)
		983	{
		984	*tr = 0.;
		985	*ti = 0.;
		986	return;
		987	}
		988
		989	*tr = fabsl ((f1r - f2r) / den);
		990	*ti = fabsl ((f1i - f2i) / den);
		991
		992	return;
		993	}
		994
		995	/-----------------------------------------------------------------------/
		996
		997	/* compute component of basis function i on segment is. */
		998	void sbf (int i, int is, double aa, double bb, double *cc)
		999	{
		1000	int ix, jsno, june, jcox, jcoxx, jend, iend, njun1 = 0, njun2;
		1001	double d, sig, pp, sdh, cdh, sd, omc, aj, pm = 0, cd, ap, qp, qm, xxi;
		1002
		1003	*aa = 0.;
		1004	*bb = 0.;
		1005	*cc = 0.;
		1006	june = 0;
		1007	jsno = 0;
		1008	pp = 0.;
		1009	ix = i - 1;
		1010
		1011	jcox = data.icon1[ix];
		1012	if (jcox > PCHCON)
		1013	jcox = i;
		1014	jcoxx = jcox - 1;
		1015
		1016	jend = -1;
		1017	iend = -1;
		1018	sig = -1.;
		1019
		1020	do
		1021	{
		1022	if (jcox != 0)
		1023	{
		1024	if (jcox < 0)
		1025	jcox = -jcox;
		1026	else
		1027	{
		1028	sig = -sig;
		1029	jend = -jend;
		1030	}
		1031
		1032	jcoxx = jcox - 1;
		1033	jsno++;
		1034	d = PI * data.si[jcoxx];
		1035	sdh = sinl (d);
		1036	cdh = cosl (d);
		1037	sd = 2. * sdh * cdh;
		1038
		1039	if (d <= 0.015)
		1040	{
		1041	omc = 4. * d * d;
		1042	omc =
		1043	((1.3888889e-3 * omc - 4.1666666667e-2) * omc + .5) * omc;
		1044	}
		1045	else
		1046	omc = 1. - cdh * cdh + sdh * sdh;
		1047
		1048	aj = 1. / (logl (1. / (PI * data.bi[jcoxx])) - .577215664);
		1049	pp -= omc / sd * aj;
		1050
		1051	if (jcox == is)
		1052	{
		1053	aa = aj / sd sig;
		1054	bb = aj / (2. cdh);
		1055	cc = -aj / (2. sdh) * sig;
		1056	june = iend;
		1057	}
		1058
		1059	if (jcox != i)
		1060	{
		1061	if (jend != 1)
		1062	jcox = data.icon1[jcoxx];
		1063	else
		1064	jcox = data.icon2[jcoxx];
		1065
		1066	if (abs (jcox) != i)
		1067	{
		1068	if (jcox == 0)
		1069	{
		1070	fprintf (
		1071	output_fp,
		1072	"\n SBF - SEGMENT CONNECTION ERROR FOR "
		1073	"SEGMENT %d",
		1074	i);
		1075	stop (-1);
		1076	}
		1077	else
		1078	continue;
		1079	}
		1080
		1081	} /* if( jcox != i ) */
		1082	else if (jcox == is)
		1083	bb = -bb;
		1084
		1085	if (iend == 1)
		1086	break;
		1087
		1088	} /* if( jcox != 0 ) */
		1089
		1090	pm = -pp;
		1091	pp = 0.;
		1092	njun1 = jsno;
		1093
		1094	jcox = data.icon2[ix];
		1095	if (jcox > PCHCON)
		1096	jcox = i;
		1097
		1098	jend = 1;
		1099	iend = 1;
		1100	sig = -1.;
		1101
		1102	} /* do */
		1103	while (jcox != 0);
		1104
		1105	njun2 = jsno - njun1;
		1106	d = PI * data.si[ix];
		1107	sdh = sinl (d);
		1108	cdh = cosl (d);
		1109	sd = 2. * sdh * cdh;
		1110	cd = cdh * cdh - sdh * sdh;
		1111
		1112	if (d <= 0.015)
		1113	{
		1114	omc = 4. * d * d;
		1115	omc = ((1.3888889e-3 * omc - 4.1666666667e-2) * omc + .5) * omc;
		1116	}
		1117	else
		1118	omc = 1. - cd;
		1119
		1120	ap = 1. / (logl (1. / (PI * data.bi[ix])) - .577215664);
		1121	aj = ap;
		1122
		1123	if (njun1 == 0)
		1124	{
		1125	if (njun2 == 0)
		1126	{
		1127	*aa = -1.;
		1128	qp = PI * data.bi[ix];
		1129	xxi = qp * qp;
		1130	xxi = qp * (1. - .5 * xxi) / (1. - xxi);
		1131	cc = 1. / (cdh - xxi sdh);
		1132	return;
		1133	}
		1134
		1135	qp = PI * data.bi[ix];
		1136	xxi = qp * qp;
		1137	xxi = qp * (1. - .5 * xxi) / (1. - xxi);
		1138	qp = -(omc + xxi * sd) / (sd * (ap + xxi * pp) + cd * (xxi * ap - pp));
		1139
		1140	if (june == 1)
		1141	{
		1142	aa = -aa * qp;
		1143	bb = bb * qp;
		1144	cc = -cc * qp;
		1145	if (i != is)
		1146	return;
		1147	}
		1148
		1149	*aa -= 1.;
		1150	d = cd - xxi * sd;
		1151	bb += (sdh + ap qp * (cdh - xxi * sdh)) / d;
		1152	cc += (cdh + ap qp * (sdh + xxi * cdh)) / d;
		1153	return;
		1154
		1155	} /* if( njun1 == 0) */
		1156
		1157	if (njun2 == 0)
		1158	{
		1159	qm = PI * data.bi[ix];
		1160	xxi = qm * qm;
		1161	xxi = qm * (1. - .5 * xxi) / (1. - xxi);
		1162	qm = (omc + xxi * sd) / (sd * (aj - xxi * pm) + cd * (pm + xxi * aj));
		1163
		1164	if (june == -1)
		1165	{
		1166	aa = aa * qm;
		1167	bb = bb * qm;
		1168	cc = cc * qm;
		1169	if (i != is)
		1170	return;
		1171	}
		1172
		1173	*aa -= 1.;
		1174	d = cd - xxi * sd;
		1175	bb += (aj qm * (cdh - xxi * sdh) - sdh) / d;
		1176	cc += (cdh - aj qm * (sdh + xxi * cdh)) / d;
		1177	return;
		1178
		1179	} /* if( njun2 == 0) */
		1180
		1181	qp = sd * (pm * pp + aj * ap) + cd * (pm * ap - pp * aj);
		1182	qm = (ap * omc - pp * sd) / qp;
		1183	qp = -(aj * omc + pm * sd) / qp;
		1184
		1185	if (june != 0)
		1186	{
		1187	if (june < 0)
		1188	{
		1189	aa = aa * qm;
		1190	bb = bb * qm;
		1191	cc = cc * qm;
		1192	}
		1193	else
		1194	{
		1195	aa = -aa * qp;
		1196	bb = bb * qp;
		1197	cc = -cc * qp;
		1198	}
		1199
		1200	if (i != is)
		1201	return;
		1202
		1203	} /* if( june != 0 ) */
		1204
		1205	*aa -= 1.;
		1206	bb += (aj qm + ap * qp) * sdh / sd;
		1207	cc += (aj qm - ap * qp) * cdh / sd;
		1208
		1209	return;
		1210	}
		1211
		1212	/-----------------------------------------------------------------------/
		1213
		1214	/* compute basis function i */
		1215	void tbf (int i, int icap)
		1216	{
		1217	int ix, jcox, jcoxx, jend, iend, njun1 = 0, njun2, jsnop, jsnox;
		1218	double pp, sdh, cdh, sd, omc, aj, pm = 0, cd, ap, qp, qm, xxi;
		1219	double d, sig; /* also global */
		1220
		1221	segj.jsno = 0;
		1222	pp = 0.;
		1223	ix = i - 1;
		1224	jcox = data.icon1[ix];
		1225
		1226	if (jcox > PCHCON)
		1227	jcox = i;
		1228
		1229	jend = -1;
		1230	iend = -1;
		1231	sig = -1.;
		1232
		1233	do
		1234	{
		1235	if (jcox != 0)
		1236	{
		1237	if (jcox < 0)
		1238	jcox = -jcox;
		1239	else
		1240	{
		1241	sig = -sig;
		1242	jend = -jend;
		1243	}
		1244
		1245	jcoxx = jcox - 1;
		1246	segj.jsno++;
		1247	jsnox = segj.jsno - 1;
		1248	segj.jco[jsnox] = jcox;
		1249	d = PI * data.si[jcoxx];
		1250	sdh = sinl (d);
		1251	cdh = cosl (d);
		1252	sd = 2. * sdh * cdh;
		1253
		1254	if (d <= 0.015)
		1255	{
		1256	omc = 4. * d * d;
		1257	omc =
		1258	((1.3888889e-3 * omc - 4.1666666667e-2) * omc + .5) * omc;
		1259	}
		1260	else
		1261	omc = 1. - cdh * cdh + sdh * sdh;
		1262
		1263	aj = 1. / (logl (1. / (PI * data.bi[jcoxx])) - .577215664);
		1264	pp = pp - omc / sd * aj;
		1265	segj.ax[jsnox] = aj / sd * sig;
		1266	segj.bx[jsnox] = aj / (2. * cdh);
		1267	segj.cx[jsnox] = -aj / (2. * sdh) * sig;
		1268
		1269	if (jcox != i)
		1270	{
		1271	if (jend == 1)
		1272	jcox = data.icon2[jcoxx];
		1273	else
		1274	jcox = data.icon1[jcoxx];
		1275
		1276	if (abs (jcox) != i)
		1277	{
		1278	if (jcox != 0)
		1279	continue;
		1280	else
		1281	{
		1282	fprintf (
		1283	output_fp,
		1284	"\n TBF - SEGMENT CONNECTION ERROR FOR "
		1285	"SEGMENT %5d",
		1286	i);
		1287	stop (-1);
		1288	}
		1289	}
		1290
		1291	} /* if( jcox != i) */
		1292	else
		1293	segj.bx[jsnox] = -segj.bx[jsnox];
		1294
		1295	if (iend == 1)
		1296	break;
		1297
		1298	} /* if( jcox != 0 ) */
		1299
		1300	pm = -pp;
		1301	pp = 0.;
		1302	njun1 = segj.jsno;
		1303
		1304	jcox = data.icon2[ix];
		1305	if (jcox > PCHCON)
		1306	jcox = i;
		1307
		1308	jend = 1;
		1309	iend = 1;
		1310	sig = -1.;
		1311
		1312	} /* do */
		1313	while (jcox != 0);
		1314
		1315	njun2 = segj.jsno - njun1;
		1316	jsnop = segj.jsno;
		1317	segj.jco[jsnop] = i;
		1318	d = PI * data.si[ix];
		1319	sdh = sinl (d);
		1320	cdh = cosl (d);
		1321	sd = 2. * sdh * cdh;
		1322	cd = cdh * cdh - sdh * sdh;
		1323
		1324	if (d <= 0.015)
		1325	{
		1326	omc = 4. * d * d;
		1327	omc = ((1.3888889e-3 * omc - 4.1666666667e-2) * omc + .5) * omc;
		1328	}
		1329	else
		1330	omc = 1. - cd;
		1331
		1332	ap = 1. / (logl (1. / (PI * data.bi[ix])) - .577215664);
		1333	aj = ap;
		1334
		1335	if (njun1 == 0)
		1336	{
		1337	if (njun2 == 0)
		1338	{
		1339	segj.bx[jsnop] = 0.;
		1340
		1341	if (icap == 0)
		1342	xxi = 0.;
		1343	else
		1344	{
		1345	qp = PI * data.bi[ix];
		1346	xxi = qp * qp;
		1347	xxi = qp * (1. - .5 * xxi) / (1. - xxi);
		1348	}
		1349
		1350	segj.cx[jsnop] = 1. / (cdh - xxi * sdh);
		1351	segj.jsno = jsnop + 1;
		1352	segj.ax[jsnop] = -1.;
		1353	return;
		1354
		1355	} /* if( njun2 == 0) */
		1356
		1357	if (icap == 0)
		1358	xxi = 0.;
		1359	else
		1360	{
		1361	qp = PI * data.bi[ix];
		1362	xxi = qp * qp;
		1363	xxi = qp * (1. - .5 * xxi) / (1. - xxi);
		1364	}
		1365
		1366	qp = -(omc + xxi * sd) / (sd * (ap + xxi * pp) + cd * (xxi * ap - pp));
		1367	d = cd - xxi * sd;
		1368	segj.bx[jsnop] = (sdh + ap * qp * (cdh - xxi * sdh)) / d;
		1369	segj.cx[jsnop] = (cdh + ap * qp * (sdh + xxi * cdh)) / d;
		1370
		1371	for (iend = 0; iend < njun2; iend++)
		1372	{
		1373	segj.ax[iend] = -segj.ax[iend] * qp;
		1374	segj.bx[iend] = segj.bx[iend] * qp;
		1375	segj.cx[iend] = -segj.cx[iend] * qp;
		1376	}
		1377
		1378	segj.jsno = jsnop + 1;
		1379	segj.ax[jsnop] = -1.;
		1380	return;
		1381
		1382	} /* if( njun1 == 0) */
		1383
		1384	if (njun2 == 0)
		1385	{
		1386	if (icap == 0)
		1387	xxi = 0.;
		1388	else
		1389	{
		1390	qm = PI * data.bi[ix];
		1391	xxi = qm * qm;
		1392	xxi = qm * (1. - .5 * xxi) / (1. - xxi);
		1393	}
		1394
		1395	qm = (omc + xxi * sd) / (sd * (aj - xxi * pm) + cd * (pm + xxi * aj));
		1396	d = cd - xxi * sd;
		1397	segj.bx[jsnop] = (aj * qm * (cdh - xxi * sdh) - sdh) / d;
		1398	segj.cx[jsnop] = (cdh - aj * qm * (sdh + xxi * cdh)) / d;
		1399
		1400	for (iend = 0; iend < njun1; iend++)
		1401	{
		1402	segj.ax[iend] = segj.ax[iend] * qm;
		1403	segj.bx[iend] = segj.bx[iend] * qm;
		1404	segj.cx[iend] = segj.cx[iend] * qm;
		1405	}
		1406
		1407	segj.jsno = jsnop + 1;
		1408	segj.ax[jsnop] = -1.;
		1409	return;
		1410
		1411	} /* if( njun2 == 0) */
		1412
		1413	qp = sd * (pm * pp + aj * ap) + cd * (pm * ap - pp * aj);
		1414	qm = (ap * omc - pp * sd) / qp;
		1415	qp = -(aj * omc + pm * sd) / qp;
		1416	segj.bx[jsnop] = (aj * qm + ap * qp) * sdh / sd;
		1417	segj.cx[jsnop] = (aj * qm - ap * qp) * cdh / sd;
		1418
		1419	for (iend = 0; iend < njun1; iend++)
		1420	{
		1421	segj.ax[iend] = segj.ax[iend] * qm;
		1422	segj.bx[iend] = segj.bx[iend] * qm;
		1423	segj.cx[iend] = segj.cx[iend] * qm;
		1424	}
		1425
		1426	jend = njun1;
		1427	for (iend = jend; iend < segj.jsno; iend++)
		1428	{
		1429	segj.ax[iend] = -segj.ax[iend] * qp;
		1430	segj.bx[iend] = segj.bx[iend] * qp;
		1431	segj.cx[iend] = -segj.cx[iend] * qp;
		1432	}
		1433
		1434	segj.jsno = jsnop + 1;
		1435	segj.ax[jsnop] = -1.;
		1436	}
		1437
		1438	/-----------------------------------------------------------------------/
		1439
		1440	/* compute the components of all basis functions on segment j */
		1441	void trio (int j)
		1442	{
		1443	int jcox, jcoxx, jsnox, jx, jend = 0, iend = 0;
		1444
		1445	segj.jsno = 0;
		1446	jx = j - 1;
		1447	jcox = data.icon1[jx];
		1448	jcoxx = jcox - 1;
		1449
		1450	if (jcox <= PCHCON)
		1451	{
		1452	jend = -1;
		1453	iend = -1;
		1454	}
		1455
		1456	if ((jcox == 0) \|\| (jcox > PCHCON))
		1457	{
		1458	jcox = data.icon2[jx];
		1459	jcoxx = jcox - 1;
		1460
		1461	if (jcox <= PCHCON)
		1462	{
		1463	jend = 1;
		1464	iend = 1;
		1465	}
		1466
		1467	if (jcox == 0 \|\| (jcox > PCHCON))
		1468	{
		1469	jsnox = segj.jsno;
		1470	segj.jsno++;
		1471
		1472	/* Allocate to connections buffers */
		1473	if (segj.jsno >= segj.maxcon)
		1474	{
		1475	segj.maxcon = segj.jsno + 1;
		1476	mem_realloc ((void ) &segj.jco, segj.maxcon sizeof (int));
		1477	mem_realloc ((void ) &segj.ax, segj.maxcon sizeof (double));
		1478	mem_realloc ((void ) &segj.bx, segj.maxcon sizeof (double));
		1479	mem_realloc ((void ) &segj.cx, segj.maxcon sizeof (double));
		1480	}
		1481
		1482	sbf (j, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
		1483	segj.jco[jsnox] = j;
		1484	return;
		1485	}
		1486
		1487	} /* if( (jcox == 0) \|\| (jcox > PCHCON) ) */
		1488
		1489	do
		1490	{
		1491	if (jcox < 0)
		1492	jcox = -jcox;
		1493	else
		1494	jend = -jend;
		1495	jcoxx = jcox - 1;
		1496
		1497	if (jcox != j)
		1498	{
		1499	jsnox = segj.jsno;
		1500	segj.jsno++;
		1501
		1502	/* Allocate to connections buffers */
		1503	if (segj.jsno >= segj.maxcon)
		1504	{
		1505	segj.maxcon = segj.jsno + 1;
		1506	mem_realloc ((void ) &segj.jco, segj.maxcon sizeof (int));
		1507	mem_realloc ((void ) &segj.ax, segj.maxcon sizeof (double));
		1508	mem_realloc ((void ) &segj.bx, segj.maxcon sizeof (double));
		1509	mem_realloc ((void ) &segj.cx, segj.maxcon sizeof (double));
		1510	}
		1511
		1512	sbf (jcox, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
		1513	segj.jco[jsnox] = jcox;
		1514
		1515	if (jend != 1)
		1516	jcox = data.icon1[jcoxx];
		1517	else
		1518	jcox = data.icon2[jcoxx];
		1519
		1520	if (jcox == 0)
		1521	{
		1522	fprintf (
		1523	output_fp,
		1524	"\n TRIO - SEGMENT CONNENTION ERROR FOR SEGMENT %5d",
		1525	j);
		1526	stop (-1);
		1527	}
		1528	else
		1529	continue;
		1530
		1531	} /* if( jcox != j) */
		1532
		1533	if (iend == 1)
		1534	break;
		1535
		1536	jcox = data.icon2[jx];
		1537
		1538	if (jcox > PCHCON)
		1539	break;
		1540
		1541	jend = 1;
		1542	iend = 1;
		1543
		1544	} /* do */
		1545	while (jcox != 0);
		1546
		1547	jsnox = segj.jsno;
		1548	segj.jsno++;
		1549
		1550	/* Allocate to connections buffers */
		1551	if (segj.jsno >= segj.maxcon)
		1552	{
		1553	segj.maxcon = segj.jsno + 1;
		1554	mem_realloc ((void ) &segj.jco, segj.maxcon sizeof (int));
		1555	mem_realloc ((void ) &segj.ax, segj.maxcon sizeof (double));
		1556	mem_realloc ((void ) &segj.bx, segj.maxcon sizeof (double));
		1557	mem_realloc ((void ) &segj.cx, segj.maxcon sizeof (double));
		1558	}
		1559
		1560	sbf (j, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
		1561	segj.jco[jsnox] = j;
		1562
		1563	return;
		1564	}
		1565
		1566	/-----------------------------------------------------------------------/
		1567
		1568	/* zint computes the internal impedance of a circular wire */
		1569	void zint (double sigl, double rolam, complex double *zint)
		1570	{
		1571	#define cc1 (6.0e-7 + 1.9e-6fj)
		1572	#define cc2 (-3.4e-6 + 5.1e-6fj)
		1573	#define cc3 (-2.52e-5 + 0.fj)
		1574	#define cc4 (-9.06e-5 - 9.01e-5fj)
		1575	#define cc5 (0. - 9.765e-4fj)
		1576	#define cc6 (.0110486 - .0110485fj)
		1577	#define cc7 (0. - .3926991fj)
		1578	#define cc8 (1.6e-6 - 3.2e-6fj)
		1579	#define cc9 (1.17e-5 - 2.4e-6fj)
		1580	#define cc10 (3.46e-5 + 3.38e-5fj)
		1581	#define cc11 (5.0e-7 + 2.452e-4fj)
		1582	#define cc12 (-1.3813e-3 + 1.3811e-3fj)
		1583	#define cc13 (-6.25001e-2 - 1.0e-7fj)
		1584	#define cc14 (.7071068 + .7071068fj)
		1585	#define cn cc14
		1586
		1587	#define th(d) \
		1588	((((((cc1 * (d) + cc2) * (d) + cc3) * (d) + cc4) * (d) + cc5) * (d) + cc6) * (d) + cc7)
		1589	#define ph(d) \
		1590	((((((cc8 * (d) + cc9) * (d) + cc10) * (d) + cc11) * (d) + cc12) * (d) + cc13) * \
		1591	(d) + \
		1592	cc14)
		1593	#define f(d) (csqrtl (POT / (d)) * cexpl (-cn * (d) + th (-8. / x)))
		1594	#define g(d) (cexpl (cn * (d) + th (8. / x)) / csqrtl (TP * (d)))
		1595
		1596	double x, y, s, ber, bei;
		1597	double tpcmu = 2.368705e+3;
		1598	;
		1599	double cmotp = 60.00;
		1600	complex double br1, br2;
		1601
		1602	x = sqrtl (tpcmu * sigl) * rolam;
		1603	if (x <= 110.)
		1604	{
		1605	if (x <= 8.)
		1606	{
		1607	y = x / 8.;
		1608	y = y * y;
		1609	s = y * y;
		1610
		1611	ber = ((((((-9.01e-6 * s + 1.22552e-3) * s - .08349609) * s +
		1612	2.6419140) *
		1613	s -
		1614	32.363456) *
		1615	s +
		1616	113.77778) *
		1617	s -
		1618	64.) *
		1619	s +
		1620	1.;
		1621
		1622	bei = ((((((1.1346e-4 * s - .01103667) * s + .52185615) * s -
		1623	10.567658) *
		1624	s +
		1625	72.817777) *
		1626	s -
		1627	113.77778) *
		1628	s +
		1629	16.) *
		1630	y;
		1631
		1632	br1 = cmplx (ber, bei);
		1633
		1634	ber = (((((((-3.94e-6 * s + 4.5957e-4) * s - .02609253) * s +
		1635	.66047849) *
		1636	s -
		1637	6.0681481) *
		1638	s +
		1639	14.222222) *
		1640	s -
		1641	4.) *
		1642	y) *
		1643	x;
		1644
		1645	bei = ((((((4.609e-5 * s - 3.79386e-3) * s + .14677204) * s -
		1646	2.3116751) *
		1647	s +
		1648	11.377778) *
		1649	s -
		1650	10.666667) *
		1651	s +
		1652	.5) *
		1653	x;
		1654
		1655	br2 = cmplx (ber, bei);
		1656	br1 = br1 / br2;
		1657	zint = CPLX_01 sqrtl (cmotp / sigl) * br1 / rolam;
		1658
		1659	} /* if( x <= 8.) */
		1660
		1661	br2 = CPLX_01 * f (x) / PI;
		1662	br1 = g (x) + br2;
		1663	br2 = g (x) * ph (8. / x) - br2 * ph (-8. / x);
		1664	br1 = br1 / br2;
		1665	zint = CPLX_01 sqrtl (cmotp / sigl) * br1 / rolam;
		1666
		1667	} /* if( x <= 110.) */
		1668
		1669	br1 = cmplx (.70710678, -.70710678);
		1670	zint = CPLX_01 sqrtl (cmotp / sigl) * br1 / rolam;
		1671	}
		1672
		1673	/-----------------------------------------------------------------------/
		1674
		1675	/* cang returns the phase angle of a complex number in degrees. */
		1676	double cang (complex double z)
		1677	{
		1678	return (cargl (z) * TD);
		1679	}
		1680
		1681	/-----------------------------------------------------------------------/

Subversion Repositories Nec2c

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