Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/bessel.cs @ 5138

Visit:

Last change on this file since 5138 was 3839, checked in by mkommend, 15 years ago
implemented first version of LR (ticket #1012)
File size: 50.8 KB

Rev	Line
[3839]	1	/*************************************************************************
	2	Cephes Math Library Release 2.8: June, 2000
	3	Copyright by Stephen L. Moshier
	4
	5	Contributors:
	6	* Sergey Bochkanov (ALGLIB project). Translation from C to
	7	pseudocode.
	8
	9	See subroutines comments for additional copyrights.
	10
	11	>>> SOURCE LICENSE >>>
	12	This program is free software; you can redistribute it and/or modify
	13	it under the terms of the GNU General Public License as published by
	14	the Free Software Foundation (www.fsf.org); either version 2 of the
	15	License, or (at your option) any later version.
	16
	17	This program is distributed in the hope that it will be useful,
	18	but WITHOUT ANY WARRANTY; without even the implied warranty of
	19	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	20	GNU General Public License for more details.
	21
	22	A copy of the GNU General Public License is available at
	23	http://www.fsf.org/licensing/licenses
	24
	25	>>> END OF LICENSE >>>
	26	*************************************************************************/
	27
	28	using System;
	29
	30	namespace alglib
	31	{
	32	public class bessel
	33	{
	34	/*************************************************************************
	35	Bessel function of order zero
	36
	37	Returns Bessel function of order zero of the argument.
	38
	39	The domain is divided into the intervals [0, 5] and
	40	(5, infinity). In the first interval the following rational
	41	approximation is used:
	42
	43
	44	2 2
	45	(w - r ) (w - r ) P (w) / Q (w)
	46	1 2 3 8
	47
	48	2
	49	where w = x and the two r's are zeros of the function.
	50
	51	In the second interval, the Hankel asymptotic expansion
	52	is employed with two rational functions of degree 6/6
	53	and 7/7.
	54
	55	ACCURACY:
	56
	57	Absolute error:
	58	arithmetic domain # trials peak rms
	59	IEEE 0, 30 60000 4.2e-16 1.1e-16
	60
	61	Cephes Math Library Release 2.8: June, 2000
	62	Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
	63	*************************************************************************/
	64	public static double besselj0(double x)
	65	{
	66	double result = 0;
	67	double xsq = 0;
	68	double nn = 0;
	69	double pzero = 0;
	70	double qzero = 0;
	71	double p1 = 0;
	72	double q1 = 0;
	73
	74	if( (double)(x)<(double)(0) )
	75	{
	76	x = -x;
	77	}
	78	if( (double)(x)>(double)(8.0) )
	79	{
	80	besselasympt0(x, ref pzero, ref qzero);
	81	nn = x-Math.PI/4;
	82	result = Math.Sqrt(2/Math.PI/x)(pzeroMath.Cos(nn)-qzero*Math.Sin(nn));
	83	return result;
	84	}
	85	xsq = AP.Math.Sqr(x);
	86	p1 = 26857.86856980014981415848441;
	87	p1 = -40504123.71833132706360663322+xsq*p1;
	88	p1 = 25071582855.36881945555156435+xsq*p1;
	89	p1 = -8085222034853.793871199468171+xsq*p1;
	90	p1 = 1434354939140344.111664316553+xsq*p1;
	91	p1 = -136762035308817138.6865416609+xsq*p1;
	92	p1 = 6382059341072356562.289432465+xsq*p1;
	93	p1 = -117915762910761053603.8440800+xsq*p1;
	94	p1 = 493378725179413356181.6813446+xsq*p1;
	95	q1 = 1.0;
	96	q1 = 1363.063652328970604442810507+xsq*q1;
	97	q1 = 1114636.098462985378182402543+xsq*q1;
	98	q1 = 669998767.2982239671814028660+xsq*q1;
	99	q1 = 312304311494.1213172572469442+xsq*q1;
	100	q1 = 112775673967979.8507056031594+xsq*q1;
	101	q1 = 30246356167094626.98627330784+xsq*q1;
	102	q1 = 5428918384092285160.200195092+xsq*q1;
	103	q1 = 493378725179413356211.3278438+xsq*q1;
	104	result = p1/q1;
	105	return result;
	106	}
	107
	108
	109	/*************************************************************************
	110	Bessel function of order one
	111
	112	Returns Bessel function of order one of the argument.
	113
	114	The domain is divided into the intervals [0, 8] and
	115	(8, infinity). In the first interval a 24 term Chebyshev
	116	expansion is used. In the second, the asymptotic
	117	trigonometric representation is employed using two
	118	rational functions of degree 5/5.
	119
	120	ACCURACY:
	121
	122	Absolute error:
	123	arithmetic domain # trials peak rms
	124	IEEE 0, 30 30000 2.6e-16 1.1e-16
	125
	126	Cephes Math Library Release 2.8: June, 2000
	127	Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
	128	*************************************************************************/
	129	public static double besselj1(double x)
	130	{
	131	double result = 0;
	132	double s = 0;
	133	double xsq = 0;
	134	double nn = 0;
	135	double pzero = 0;
	136	double qzero = 0;
	137	double p1 = 0;
	138	double q1 = 0;
	139
	140	s = Math.Sign(x);
	141	if( (double)(x)<(double)(0) )
	142	{
	143	x = -x;
	144	}
	145	if( (double)(x)>(double)(8.0) )
	146	{
	147	besselasympt1(x, ref pzero, ref qzero);
	148	nn = x-3*Math.PI/4;
	149	result = Math.Sqrt(2/Math.PI/x)(pzeroMath.Cos(nn)-qzero*Math.Sin(nn));
	150	if( (double)(s)<(double)(0) )
	151	{
	152	result = -result;
	153	}
	154	return result;
	155	}
	156	xsq = AP.Math.Sqr(x);
	157	p1 = 2701.122710892323414856790990;
	158	p1 = -4695753.530642995859767162166+xsq*p1;
	159	p1 = 3413234182.301700539091292655+xsq*p1;
	160	p1 = -1322983480332.126453125473247+xsq*p1;
	161	p1 = 290879526383477.5409737601689+xsq*p1;
	162	p1 = -35888175699101060.50743641413+xsq*p1;
	163	p1 = 2316433580634002297.931815435+xsq*p1;
	164	p1 = -66721065689249162980.20941484+xsq*p1;
	165	p1 = 581199354001606143928.050809+xsq*p1;
	166	q1 = 1.0;
	167	q1 = 1606.931573481487801970916749+xsq*q1;
	168	q1 = 1501793.594998585505921097578+xsq*q1;
	169	q1 = 1013863514.358673989967045588+xsq*q1;
	170	q1 = 524371026216.7649715406728642+xsq*q1;
	171	q1 = 208166122130760.7351240184229+xsq*q1;
	172	q1 = 60920613989175217.46105196863+xsq*q1;
	173	q1 = 11857707121903209998.37113348+xsq*q1;
	174	q1 = 1162398708003212287858.529400+xsq*q1;
	175	result = sxp1/q1;
	176	return result;
	177	}
	178
	179
	180	/*************************************************************************
	181	Bessel function of integer order
	182
	183	Returns Bessel function of order n, where n is a
	184	(possibly negative) integer.
	185
	186	The ratio of jn(x) to j0(x) is computed by backward
	187	recurrence. First the ratio jn/jn-1 is found by a
	188	continued fraction expansion. Then the recurrence
	189	relating successive orders is applied until j0 or j1 is
	190	reached.
	191
	192	If n = 0 or 1 the routine for j0 or j1 is called
	193	directly.
	194
	195	ACCURACY:
	196
	197	Absolute error:
	198	arithmetic range # trials peak rms
	199	IEEE 0, 30 5000 4.4e-16 7.9e-17
	200
	201
	202	Not suitable for large n or x. Use jv() (fractional order) instead.
	203
	204	Cephes Math Library Release 2.8: June, 2000
	205	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	206	*************************************************************************/
	207	public static double besseljn(int n,
	208	double x)
	209	{
	210	double result = 0;
	211	double pkm2 = 0;
	212	double pkm1 = 0;
	213	double pk = 0;
	214	double xk = 0;
	215	double r = 0;
	216	double ans = 0;
	217	int k = 0;
	218	int sg = 0;
	219
	220	if( n<0 )
	221	{
	222	n = -n;
	223	if( n%2==0 )
	224	{
	225	sg = 1;
	226	}
	227	else
	228	{
	229	sg = -1;
	230	}
	231	}
	232	else
	233	{
	234	sg = 1;
	235	}
	236	if( (double)(x)<(double)(0) )
	237	{
	238	if( n%2!=0 )
	239	{
	240	sg = -sg;
	241	}
	242	x = -x;
	243	}
	244	if( n==0 )
	245	{
	246	result = sg*besselj0(x);
	247	return result;
	248	}
	249	if( n==1 )
	250	{
	251	result = sg*besselj1(x);
	252	return result;
	253	}
	254	if( n==2 )
	255	{
	256	if( (double)(x)==(double)(0) )
	257	{
	258	result = 0;
	259	}
	260	else
	261	{
	262	result = sg(2.0besselj1(x)/x-besselj0(x));
	263	}
	264	return result;
	265	}
	266	if( (double)(x)<(double)(AP.Math.MachineEpsilon) )
	267	{
	268	result = 0;
	269	return result;
	270	}
	271	k = 53;
	272	pk = 2*(n+k);
	273	ans = pk;
	274	xk = x*x;
	275	do
	276	{
	277	pk = pk-2.0;
	278	ans = pk-xk/ans;
	279	k = k-1;
	280	}
	281	while( k!=0 );
	282	ans = x/ans;
	283	pk = 1.0;
	284	pkm1 = 1.0/ans;
	285	k = n-1;
	286	r = 2*k;
	287	do
	288	{
	289	pkm2 = (pkm1r-pkx)/x;
	290	pk = pkm1;
	291	pkm1 = pkm2;
	292	r = r-2.0;
	293	k = k-1;
	294	}
	295	while( k!=0 );
	296	if( (double)(Math.Abs(pk))>(double)(Math.Abs(pkm1)) )
	297	{
	298	ans = besselj1(x)/pk;
	299	}
	300	else
	301	{
	302	ans = besselj0(x)/pkm1;
	303	}
	304	result = sg*ans;
	305	return result;
	306	}
	307
	308
	309	/*************************************************************************
	310	Bessel function of the second kind, order zero
	311
	312	Returns Bessel function of the second kind, of order
	313	zero, of the argument.
	314
	315	The domain is divided into the intervals [0, 5] and
	316	(5, infinity). In the first interval a rational approximation
	317	R(x) is employed to compute
	318	y0(x) = R(x) + 2 * log(x) * j0(x) / PI.
	319	Thus a call to j0() is required.
	320
	321	In the second interval, the Hankel asymptotic expansion
	322	is employed with two rational functions of degree 6/6
	323	and 7/7.
	324
	325
	326
	327	ACCURACY:
	328
	329	Absolute error, when y0(x) < 1; else relative error:
	330
	331	arithmetic domain # trials peak rms
	332	IEEE 0, 30 30000 1.3e-15 1.6e-16
	333
	334	Cephes Math Library Release 2.8: June, 2000
	335	Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
	336	*************************************************************************/
	337	public static double bessely0(double x)
	338	{
	339	double result = 0;
	340	double nn = 0;
	341	double xsq = 0;
	342	double pzero = 0;
	343	double qzero = 0;
	344	double p4 = 0;
	345	double q4 = 0;
	346
	347	if( (double)(x)>(double)(8.0) )
	348	{
	349	besselasympt0(x, ref pzero, ref qzero);
	350	nn = x-Math.PI/4;
	351	result = Math.Sqrt(2/Math.PI/x)(pzeroMath.Sin(nn)+qzero*Math.Cos(nn));
	352	return result;
	353	}
	354	xsq = AP.Math.Sqr(x);
	355	p4 = -41370.35497933148554125235152;
	356	p4 = 59152134.65686889654273830069+xsq*p4;
	357	p4 = -34363712229.79040378171030138+xsq*p4;
	358	p4 = 10255208596863.94284509167421+xsq*p4;
	359	p4 = -1648605817185729.473122082537+xsq*p4;
	360	p4 = 137562431639934407.8571335453+xsq*p4;
	361	p4 = -5247065581112764941.297350814+xsq*p4;
	362	p4 = 65874732757195549259.99402049+xsq*p4;
	363	p4 = -27502866786291095837.01933175+xsq*p4;
	364	q4 = 1.0;
	365	q4 = 1282.452772478993804176329391+xsq*q4;
	366	q4 = 1001702.641288906265666651753+xsq*q4;
	367	q4 = 579512264.0700729537480087915+xsq*q4;
	368	q4 = 261306575504.1081249568482092+xsq*q4;
	369	q4 = 91620380340751.85262489147968+xsq*q4;
	370	q4 = 23928830434997818.57439356652+xsq*q4;
	371	q4 = 4192417043410839973.904769661+xsq*q4;
	372	q4 = 372645883898616588198.9980+xsq*q4;
	373	result = p4/q4+2/Math.PIbesselj0(x)Math.Log(x);
	374	return result;
	375	}
	376
	377
	378	/*************************************************************************
	379	Bessel function of second kind of order one
	380
	381	Returns Bessel function of the second kind of order one
	382	of the argument.
	383
	384	The domain is divided into the intervals [0, 8] and
	385	(8, infinity). In the first interval a 25 term Chebyshev
	386	expansion is used, and a call to j1() is required.
	387	In the second, the asymptotic trigonometric representation
	388	is employed using two rational functions of degree 5/5.
	389
	390	ACCURACY:
	391
	392	Absolute error:
	393	arithmetic domain # trials peak rms
	394	IEEE 0, 30 30000 1.0e-15 1.3e-16
	395
	396	Cephes Math Library Release 2.8: June, 2000
	397	Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
	398	*************************************************************************/
	399	public static double bessely1(double x)
	400	{
	401	double result = 0;
	402	double nn = 0;
	403	double xsq = 0;
	404	double pzero = 0;
	405	double qzero = 0;
	406	double p4 = 0;
	407	double q4 = 0;
	408
	409	if( (double)(x)>(double)(8.0) )
	410	{
	411	besselasympt1(x, ref pzero, ref qzero);
	412	nn = x-3*Math.PI/4;
	413	result = Math.Sqrt(2/Math.PI/x)(pzeroMath.Sin(nn)+qzero*Math.Cos(nn));
	414	return result;
	415	}
	416	xsq = AP.Math.Sqr(x);
	417	p4 = -2108847.540133123652824139923;
	418	p4 = 3639488548.124002058278999428+xsq*p4;
	419	p4 = -2580681702194.450950541426399+xsq*p4;
	420	p4 = 956993023992168.3481121552788+xsq*p4;
	421	p4 = -196588746272214065.8820322248+xsq*p4;
	422	p4 = 21931073399177975921.11427556+xsq*p4;
	423	p4 = -1212297555414509577913.561535+xsq*p4;
	424	p4 = 26554738314348543268942.48968+xsq*p4;
	425	p4 = -99637534243069222259967.44354+xsq*p4;
	426	q4 = 1.0;
	427	q4 = 1612.361029677000859332072312+xsq*q4;
	428	q4 = 1563282.754899580604737366452+xsq*q4;
	429	q4 = 1128686837.169442121732366891+xsq*q4;
	430	q4 = 646534088126.5275571961681500+xsq*q4;
	431	q4 = 297663212564727.6729292742282+xsq*q4;
	432	q4 = 108225825940881955.2553850180+xsq*q4;
	433	q4 = 29549879358971486742.90758119+xsq*q4;
	434	q4 = 5435310377188854170800.653097+xsq*q4;
	435	q4 = 508206736694124324531442.4152+xsq*q4;
	436	result = xp4/q4+2/Math.PI(besselj1(x)*Math.Log(x)-1/x);
	437	return result;
	438	}
	439
	440
	441	/*************************************************************************
	442	Bessel function of second kind of integer order
	443
	444	Returns Bessel function of order n, where n is a
	445	(possibly negative) integer.
	446
	447	The function is evaluated by forward recurrence on
	448	n, starting with values computed by the routines
	449	y0() and y1().
	450
	451	If n = 0 or 1 the routine for y0 or y1 is called
	452	directly.
	453
	454	ACCURACY:
	455	Absolute error, except relative
	456	when y > 1:
	457	arithmetic domain # trials peak rms
	458	IEEE 0, 30 30000 3.4e-15 4.3e-16
	459
	460	Cephes Math Library Release 2.8: June, 2000
	461	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	462	*************************************************************************/
	463	public static double besselyn(int n,
	464	double x)
	465	{
	466	double result = 0;
	467	int i = 0;
	468	double a = 0;
	469	double b = 0;
	470	double tmp = 0;
	471	double s = 0;
	472
	473	s = 1;
	474	if( n<0 )
	475	{
	476	n = -n;
	477	if( n%2!=0 )
	478	{
	479	s = -1;
	480	}
	481	}
	482	if( n==0 )
	483	{
	484	result = bessely0(x);
	485	return result;
	486	}
	487	if( n==1 )
	488	{
	489	result = s*bessely1(x);
	490	return result;
	491	}
	492	a = bessely0(x);
	493	b = bessely1(x);
	494	for(i=1; i<=n-1; i++)
	495	{
	496	tmp = b;
	497	b = 2i/xb-a;
	498	a = tmp;
	499	}
	500	result = s*b;
	501	return result;
	502	}
	503
	504
	505	/*************************************************************************
	506	Modified Bessel function of order zero
	507
	508	Returns modified Bessel function of order zero of the
	509	argument.
	510
	511	The function is defined as i0(x) = j0( ix ).
	512
	513	The range is partitioned into the two intervals [0,8] and
	514	(8, infinity). Chebyshev polynomial expansions are employed
	515	in each interval.
	516
	517	ACCURACY:
	518
	519	Relative error:
	520	arithmetic domain # trials peak rms
	521	IEEE 0,30 30000 5.8e-16 1.4e-16
	522
	523	Cephes Math Library Release 2.8: June, 2000
	524	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	525	*************************************************************************/
	526	public static double besseli0(double x)
	527	{
	528	double result = 0;
	529	double y = 0;
	530	double v = 0;
	531	double z = 0;
	532	double b0 = 0;
	533	double b1 = 0;
	534	double b2 = 0;
	535
	536	if( (double)(x)<(double)(0) )
	537	{
	538	x = -x;
	539	}
	540	if( (double)(x)<=(double)(8.0) )
	541	{
	542	y = x/2.0-2.0;
	543	besselmfirstcheb(-4.41534164647933937950E-18, ref b0, ref b1, ref b2);
	544	besselmnextcheb(y, 3.33079451882223809783E-17, ref b0, ref b1, ref b2);
	545	besselmnextcheb(y, -2.43127984654795469359E-16, ref b0, ref b1, ref b2);
	546	besselmnextcheb(y, 1.71539128555513303061E-15, ref b0, ref b1, ref b2);
	547	besselmnextcheb(y, -1.16853328779934516808E-14, ref b0, ref b1, ref b2);
	548	besselmnextcheb(y, 7.67618549860493561688E-14, ref b0, ref b1, ref b2);
	549	besselmnextcheb(y, -4.85644678311192946090E-13, ref b0, ref b1, ref b2);
	550	besselmnextcheb(y, 2.95505266312963983461E-12, ref b0, ref b1, ref b2);
	551	besselmnextcheb(y, -1.72682629144155570723E-11, ref b0, ref b1, ref b2);
	552	besselmnextcheb(y, 9.67580903537323691224E-11, ref b0, ref b1, ref b2);
	553	besselmnextcheb(y, -5.18979560163526290666E-10, ref b0, ref b1, ref b2);
	554	besselmnextcheb(y, 2.65982372468238665035E-9, ref b0, ref b1, ref b2);
	555	besselmnextcheb(y, -1.30002500998624804212E-8, ref b0, ref b1, ref b2);
	556	besselmnextcheb(y, 6.04699502254191894932E-8, ref b0, ref b1, ref b2);
	557	besselmnextcheb(y, -2.67079385394061173391E-7, ref b0, ref b1, ref b2);
	558	besselmnextcheb(y, 1.11738753912010371815E-6, ref b0, ref b1, ref b2);
	559	besselmnextcheb(y, -4.41673835845875056359E-6, ref b0, ref b1, ref b2);
	560	besselmnextcheb(y, 1.64484480707288970893E-5, ref b0, ref b1, ref b2);
	561	besselmnextcheb(y, -5.75419501008210370398E-5, ref b0, ref b1, ref b2);
	562	besselmnextcheb(y, 1.88502885095841655729E-4, ref b0, ref b1, ref b2);
	563	besselmnextcheb(y, -5.76375574538582365885E-4, ref b0, ref b1, ref b2);
	564	besselmnextcheb(y, 1.63947561694133579842E-3, ref b0, ref b1, ref b2);
	565	besselmnextcheb(y, -4.32430999505057594430E-3, ref b0, ref b1, ref b2);
	566	besselmnextcheb(y, 1.05464603945949983183E-2, ref b0, ref b1, ref b2);
	567	besselmnextcheb(y, -2.37374148058994688156E-2, ref b0, ref b1, ref b2);
	568	besselmnextcheb(y, 4.93052842396707084878E-2, ref b0, ref b1, ref b2);
	569	besselmnextcheb(y, -9.49010970480476444210E-2, ref b0, ref b1, ref b2);
	570	besselmnextcheb(y, 1.71620901522208775349E-1, ref b0, ref b1, ref b2);
	571	besselmnextcheb(y, -3.04682672343198398683E-1, ref b0, ref b1, ref b2);
	572	besselmnextcheb(y, 6.76795274409476084995E-1, ref b0, ref b1, ref b2);
	573	v = 0.5*(b0-b2);
	574	result = Math.Exp(x)*v;
	575	return result;
	576	}
	577	z = 32.0/x-2.0;
	578	besselmfirstcheb(-7.23318048787475395456E-18, ref b0, ref b1, ref b2);
	579	besselmnextcheb(z, -4.83050448594418207126E-18, ref b0, ref b1, ref b2);
	580	besselmnextcheb(z, 4.46562142029675999901E-17, ref b0, ref b1, ref b2);
	581	besselmnextcheb(z, 3.46122286769746109310E-17, ref b0, ref b1, ref b2);
	582	besselmnextcheb(z, -2.82762398051658348494E-16, ref b0, ref b1, ref b2);
	583	besselmnextcheb(z, -3.42548561967721913462E-16, ref b0, ref b1, ref b2);
	584	besselmnextcheb(z, 1.77256013305652638360E-15, ref b0, ref b1, ref b2);
	585	besselmnextcheb(z, 3.81168066935262242075E-15, ref b0, ref b1, ref b2);
	586	besselmnextcheb(z, -9.55484669882830764870E-15, ref b0, ref b1, ref b2);
	587	besselmnextcheb(z, -4.15056934728722208663E-14, ref b0, ref b1, ref b2);
	588	besselmnextcheb(z, 1.54008621752140982691E-14, ref b0, ref b1, ref b2);
	589	besselmnextcheb(z, 3.85277838274214270114E-13, ref b0, ref b1, ref b2);
	590	besselmnextcheb(z, 7.18012445138366623367E-13, ref b0, ref b1, ref b2);
	591	besselmnextcheb(z, -1.79417853150680611778E-12, ref b0, ref b1, ref b2);
	592	besselmnextcheb(z, -1.32158118404477131188E-11, ref b0, ref b1, ref b2);
	593	besselmnextcheb(z, -3.14991652796324136454E-11, ref b0, ref b1, ref b2);
	594	besselmnextcheb(z, 1.18891471078464383424E-11, ref b0, ref b1, ref b2);
	595	besselmnextcheb(z, 4.94060238822496958910E-10, ref b0, ref b1, ref b2);
	596	besselmnextcheb(z, 3.39623202570838634515E-9, ref b0, ref b1, ref b2);
	597	besselmnextcheb(z, 2.26666899049817806459E-8, ref b0, ref b1, ref b2);
	598	besselmnextcheb(z, 2.04891858946906374183E-7, ref b0, ref b1, ref b2);
	599	besselmnextcheb(z, 2.89137052083475648297E-6, ref b0, ref b1, ref b2);
	600	besselmnextcheb(z, 6.88975834691682398426E-5, ref b0, ref b1, ref b2);
	601	besselmnextcheb(z, 3.36911647825569408990E-3, ref b0, ref b1, ref b2);
	602	besselmnextcheb(z, 8.04490411014108831608E-1, ref b0, ref b1, ref b2);
	603	v = 0.5*(b0-b2);
	604	result = Math.Exp(x)*v/Math.Sqrt(x);
	605	return result;
	606	}
	607
	608
	609	/*************************************************************************
	610	Modified Bessel function of order one
	611
	612	Returns modified Bessel function of order one of the
	613	argument.
	614
	615	The function is defined as i1(x) = -i j1( ix ).
	616
	617	The range is partitioned into the two intervals [0,8] and
	618	(8, infinity). Chebyshev polynomial expansions are employed
	619	in each interval.
	620
	621	ACCURACY:
	622
	623	Relative error:
	624	arithmetic domain # trials peak rms
	625	IEEE 0, 30 30000 1.9e-15 2.1e-16
	626
	627	Cephes Math Library Release 2.8: June, 2000
	628	Copyright 1985, 1987, 2000 by Stephen L. Moshier
	629	*************************************************************************/
	630	public static double besseli1(double x)
	631	{
	632	double result = 0;
	633	double y = 0;
	634	double z = 0;
	635	double v = 0;
	636	double b0 = 0;
	637	double b1 = 0;
	638	double b2 = 0;
	639
	640	z = Math.Abs(x);
	641	if( (double)(z)<=(double)(8.0) )
	642	{
	643	y = z/2.0-2.0;
	644	besselm1firstcheb(2.77791411276104639959E-18, ref b0, ref b1, ref b2);
	645	besselm1nextcheb(y, -2.11142121435816608115E-17, ref b0, ref b1, ref b2);
	646	besselm1nextcheb(y, 1.55363195773620046921E-16, ref b0, ref b1, ref b2);
	647	besselm1nextcheb(y, -1.10559694773538630805E-15, ref b0, ref b1, ref b2);
	648	besselm1nextcheb(y, 7.60068429473540693410E-15, ref b0, ref b1, ref b2);
	649	besselm1nextcheb(y, -5.04218550472791168711E-14, ref b0, ref b1, ref b2);
	650	besselm1nextcheb(y, 3.22379336594557470981E-13, ref b0, ref b1, ref b2);
	651	besselm1nextcheb(y, -1.98397439776494371520E-12, ref b0, ref b1, ref b2);
	652	besselm1nextcheb(y, 1.17361862988909016308E-11, ref b0, ref b1, ref b2);
	653	besselm1nextcheb(y, -6.66348972350202774223E-11, ref b0, ref b1, ref b2);
	654	besselm1nextcheb(y, 3.62559028155211703701E-10, ref b0, ref b1, ref b2);
	655	besselm1nextcheb(y, -1.88724975172282928790E-9, ref b0, ref b1, ref b2);
	656	besselm1nextcheb(y, 9.38153738649577178388E-9, ref b0, ref b1, ref b2);
	657	besselm1nextcheb(y, -4.44505912879632808065E-8, ref b0, ref b1, ref b2);
	658	besselm1nextcheb(y, 2.00329475355213526229E-7, ref b0, ref b1, ref b2);
	659	besselm1nextcheb(y, -8.56872026469545474066E-7, ref b0, ref b1, ref b2);
	660	besselm1nextcheb(y, 3.47025130813767847674E-6, ref b0, ref b1, ref b2);
	661	besselm1nextcheb(y, -1.32731636560394358279E-5, ref b0, ref b1, ref b2);
	662	besselm1nextcheb(y, 4.78156510755005422638E-5, ref b0, ref b1, ref b2);
	663	besselm1nextcheb(y, -1.61760815825896745588E-4, ref b0, ref b1, ref b2);
	664	besselm1nextcheb(y, 5.12285956168575772895E-4, ref b0, ref b1, ref b2);
	665	besselm1nextcheb(y, -1.51357245063125314899E-3, ref b0, ref b1, ref b2);
	666	besselm1nextcheb(y, 4.15642294431288815669E-3, ref b0, ref b1, ref b2);
	667	besselm1nextcheb(y, -1.05640848946261981558E-2, ref b0, ref b1, ref b2);
	668	besselm1nextcheb(y, 2.47264490306265168283E-2, ref b0, ref b1, ref b2);
	669	besselm1nextcheb(y, -5.29459812080949914269E-2, ref b0, ref b1, ref b2);
	670	besselm1nextcheb(y, 1.02643658689847095384E-1, ref b0, ref b1, ref b2);
	671	besselm1nextcheb(y, -1.76416518357834055153E-1, ref b0, ref b1, ref b2);
	672	besselm1nextcheb(y, 2.52587186443633654823E-1, ref b0, ref b1, ref b2);
	673	v = 0.5*(b0-b2);
	674	z = vzMath.Exp(z);
	675	}
	676	else
	677	{
	678	y = 32.0/z-2.0;
	679	besselm1firstcheb(7.51729631084210481353E-18, ref b0, ref b1, ref b2);
	680	besselm1nextcheb(y, 4.41434832307170791151E-18, ref b0, ref b1, ref b2);
	681	besselm1nextcheb(y, -4.65030536848935832153E-17, ref b0, ref b1, ref b2);
	682	besselm1nextcheb(y, -3.20952592199342395980E-17, ref b0, ref b1, ref b2);
	683	besselm1nextcheb(y, 2.96262899764595013876E-16, ref b0, ref b1, ref b2);
	684	besselm1nextcheb(y, 3.30820231092092828324E-16, ref b0, ref b1, ref b2);
	685	besselm1nextcheb(y, -1.88035477551078244854E-15, ref b0, ref b1, ref b2);
	686	besselm1nextcheb(y, -3.81440307243700780478E-15, ref b0, ref b1, ref b2);
	687	besselm1nextcheb(y, 1.04202769841288027642E-14, ref b0, ref b1, ref b2);
	688	besselm1nextcheb(y, 4.27244001671195135429E-14, ref b0, ref b1, ref b2);
	689	besselm1nextcheb(y, -2.10154184277266431302E-14, ref b0, ref b1, ref b2);
	690	besselm1nextcheb(y, -4.08355111109219731823E-13, ref b0, ref b1, ref b2);
	691	besselm1nextcheb(y, -7.19855177624590851209E-13, ref b0, ref b1, ref b2);
	692	besselm1nextcheb(y, 2.03562854414708950722E-12, ref b0, ref b1, ref b2);
	693	besselm1nextcheb(y, 1.41258074366137813316E-11, ref b0, ref b1, ref b2);
	694	besselm1nextcheb(y, 3.25260358301548823856E-11, ref b0, ref b1, ref b2);
	695	besselm1nextcheb(y, -1.89749581235054123450E-11, ref b0, ref b1, ref b2);
	696	besselm1nextcheb(y, -5.58974346219658380687E-10, ref b0, ref b1, ref b2);
	697	besselm1nextcheb(y, -3.83538038596423702205E-9, ref b0, ref b1, ref b2);
	698	besselm1nextcheb(y, -2.63146884688951950684E-8, ref b0, ref b1, ref b2);
	699	besselm1nextcheb(y, -2.51223623787020892529E-7, ref b0, ref b1, ref b2);
	700	besselm1nextcheb(y, -3.88256480887769039346E-6, ref b0, ref b1, ref b2);
	701	besselm1nextcheb(y, -1.10588938762623716291E-4, ref b0, ref b1, ref b2);
	702	besselm1nextcheb(y, -9.76109749136146840777E-3, ref b0, ref b1, ref b2);
	703	besselm1nextcheb(y, 7.78576235018280120474E-1, ref b0, ref b1, ref b2);
	704	v = 0.5*(b0-b2);
	705	z = v*Math.Exp(z)/Math.Sqrt(z);
	706	}
	707	if( (double)(x)<(double)(0) )
	708	{
	709	z = -z;
	710	}
	711	result = z;
	712	return result;
	713	}
	714
	715
	716	/*************************************************************************
	717	Modified Bessel function, second kind, order zero
	718
	719	Returns modified Bessel function of the second kind
	720	of order zero of the argument.
	721
	722	The range is partitioned into the two intervals [0,8] and
	723	(8, infinity). Chebyshev polynomial expansions are employed
	724	in each interval.
	725
	726	ACCURACY:
	727
	728	Tested at 2000 random points between 0 and 8. Peak absolute
	729	error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
	730	Relative error:
	731	arithmetic domain # trials peak rms
	732	IEEE 0, 30 30000 1.2e-15 1.6e-16
	733
	734	Cephes Math Library Release 2.8: June, 2000
	735	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	736	*************************************************************************/
	737	public static double besselk0(double x)
	738	{
	739	double result = 0;
	740	double y = 0;
	741	double z = 0;
	742	double v = 0;
	743	double b0 = 0;
	744	double b1 = 0;
	745	double b2 = 0;
	746
	747	System.Diagnostics.Debug.Assert((double)(x)>(double)(0), "Domain error in BesselK0: x<=0");
	748	if( (double)(x)<=(double)(2) )
	749	{
	750	y = x*x-2.0;
	751	besselmfirstcheb(1.37446543561352307156E-16, ref b0, ref b1, ref b2);
	752	besselmnextcheb(y, 4.25981614279661018399E-14, ref b0, ref b1, ref b2);
	753	besselmnextcheb(y, 1.03496952576338420167E-11, ref b0, ref b1, ref b2);
	754	besselmnextcheb(y, 1.90451637722020886025E-9, ref b0, ref b1, ref b2);
	755	besselmnextcheb(y, 2.53479107902614945675E-7, ref b0, ref b1, ref b2);
	756	besselmnextcheb(y, 2.28621210311945178607E-5, ref b0, ref b1, ref b2);
	757	besselmnextcheb(y, 1.26461541144692592338E-3, ref b0, ref b1, ref b2);
	758	besselmnextcheb(y, 3.59799365153615016266E-2, ref b0, ref b1, ref b2);
	759	besselmnextcheb(y, 3.44289899924628486886E-1, ref b0, ref b1, ref b2);
	760	besselmnextcheb(y, -5.35327393233902768720E-1, ref b0, ref b1, ref b2);
	761	v = 0.5*(b0-b2);
	762	v = v-Math.Log(0.5x)besseli0(x);
	763	}
	764	else
	765	{
	766	z = 8.0/x-2.0;
	767	besselmfirstcheb(5.30043377268626276149E-18, ref b0, ref b1, ref b2);
	768	besselmnextcheb(z, -1.64758043015242134646E-17, ref b0, ref b1, ref b2);
	769	besselmnextcheb(z, 5.21039150503902756861E-17, ref b0, ref b1, ref b2);
	770	besselmnextcheb(z, -1.67823109680541210385E-16, ref b0, ref b1, ref b2);
	771	besselmnextcheb(z, 5.51205597852431940784E-16, ref b0, ref b1, ref b2);
	772	besselmnextcheb(z, -1.84859337734377901440E-15, ref b0, ref b1, ref b2);
	773	besselmnextcheb(z, 6.34007647740507060557E-15, ref b0, ref b1, ref b2);
	774	besselmnextcheb(z, -2.22751332699166985548E-14, ref b0, ref b1, ref b2);
	775	besselmnextcheb(z, 8.03289077536357521100E-14, ref b0, ref b1, ref b2);
	776	besselmnextcheb(z, -2.98009692317273043925E-13, ref b0, ref b1, ref b2);
	777	besselmnextcheb(z, 1.14034058820847496303E-12, ref b0, ref b1, ref b2);
	778	besselmnextcheb(z, -4.51459788337394416547E-12, ref b0, ref b1, ref b2);
	779	besselmnextcheb(z, 1.85594911495471785253E-11, ref b0, ref b1, ref b2);
	780	besselmnextcheb(z, -7.95748924447710747776E-11, ref b0, ref b1, ref b2);
	781	besselmnextcheb(z, 3.57739728140030116597E-10, ref b0, ref b1, ref b2);
	782	besselmnextcheb(z, -1.69753450938905987466E-9, ref b0, ref b1, ref b2);
	783	besselmnextcheb(z, 8.57403401741422608519E-9, ref b0, ref b1, ref b2);
	784	besselmnextcheb(z, -4.66048989768794782956E-8, ref b0, ref b1, ref b2);
	785	besselmnextcheb(z, 2.76681363944501510342E-7, ref b0, ref b1, ref b2);
	786	besselmnextcheb(z, -1.83175552271911948767E-6, ref b0, ref b1, ref b2);
	787	besselmnextcheb(z, 1.39498137188764993662E-5, ref b0, ref b1, ref b2);
	788	besselmnextcheb(z, -1.28495495816278026384E-4, ref b0, ref b1, ref b2);
	789	besselmnextcheb(z, 1.56988388573005337491E-3, ref b0, ref b1, ref b2);
	790	besselmnextcheb(z, -3.14481013119645005427E-2, ref b0, ref b1, ref b2);
	791	besselmnextcheb(z, 2.44030308206595545468E0, ref b0, ref b1, ref b2);
	792	v = 0.5*(b0-b2);
	793	v = v*Math.Exp(-x)/Math.Sqrt(x);
	794	}
	795	result = v;
	796	return result;
	797	}
	798
	799
	800	/*************************************************************************
	801	Modified Bessel function, second kind, order one
	802
	803	Computes the modified Bessel function of the second kind
	804	of order one of the argument.
	805
	806	The range is partitioned into the two intervals [0,2] and
	807	(2, infinity). Chebyshev polynomial expansions are employed
	808	in each interval.
	809
	810	ACCURACY:
	811
	812	Relative error:
	813	arithmetic domain # trials peak rms
	814	IEEE 0, 30 30000 1.2e-15 1.6e-16
	815
	816	Cephes Math Library Release 2.8: June, 2000
	817	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	818	*************************************************************************/
	819	public static double besselk1(double x)
	820	{
	821	double result = 0;
	822	double y = 0;
	823	double z = 0;
	824	double v = 0;
	825	double b0 = 0;
	826	double b1 = 0;
	827	double b2 = 0;
	828
	829	z = 0.5*x;
	830	System.Diagnostics.Debug.Assert((double)(z)>(double)(0), "Domain error in K1");
	831	if( (double)(x)<=(double)(2) )
	832	{
	833	y = x*x-2.0;
	834	besselm1firstcheb(-7.02386347938628759343E-18, ref b0, ref b1, ref b2);
	835	besselm1nextcheb(y, -2.42744985051936593393E-15, ref b0, ref b1, ref b2);
	836	besselm1nextcheb(y, -6.66690169419932900609E-13, ref b0, ref b1, ref b2);
	837	besselm1nextcheb(y, -1.41148839263352776110E-10, ref b0, ref b1, ref b2);
	838	besselm1nextcheb(y, -2.21338763073472585583E-8, ref b0, ref b1, ref b2);
	839	besselm1nextcheb(y, -2.43340614156596823496E-6, ref b0, ref b1, ref b2);
	840	besselm1nextcheb(y, -1.73028895751305206302E-4, ref b0, ref b1, ref b2);
	841	besselm1nextcheb(y, -6.97572385963986435018E-3, ref b0, ref b1, ref b2);
	842	besselm1nextcheb(y, -1.22611180822657148235E-1, ref b0, ref b1, ref b2);
	843	besselm1nextcheb(y, -3.53155960776544875667E-1, ref b0, ref b1, ref b2);
	844	besselm1nextcheb(y, 1.52530022733894777053E0, ref b0, ref b1, ref b2);
	845	v = 0.5*(b0-b2);
	846	result = Math.Log(z)*besseli1(x)+v/x;
	847	}
	848	else
	849	{
	850	y = 8.0/x-2.0;
	851	besselm1firstcheb(-5.75674448366501715755E-18, ref b0, ref b1, ref b2);
	852	besselm1nextcheb(y, 1.79405087314755922667E-17, ref b0, ref b1, ref b2);
	853	besselm1nextcheb(y, -5.68946255844285935196E-17, ref b0, ref b1, ref b2);
	854	besselm1nextcheb(y, 1.83809354436663880070E-16, ref b0, ref b1, ref b2);
	855	besselm1nextcheb(y, -6.05704724837331885336E-16, ref b0, ref b1, ref b2);
	856	besselm1nextcheb(y, 2.03870316562433424052E-15, ref b0, ref b1, ref b2);
	857	besselm1nextcheb(y, -7.01983709041831346144E-15, ref b0, ref b1, ref b2);
	858	besselm1nextcheb(y, 2.47715442448130437068E-14, ref b0, ref b1, ref b2);
	859	besselm1nextcheb(y, -8.97670518232499435011E-14, ref b0, ref b1, ref b2);
	860	besselm1nextcheb(y, 3.34841966607842919884E-13, ref b0, ref b1, ref b2);
	861	besselm1nextcheb(y, -1.28917396095102890680E-12, ref b0, ref b1, ref b2);
	862	besselm1nextcheb(y, 5.13963967348173025100E-12, ref b0, ref b1, ref b2);
	863	besselm1nextcheb(y, -2.12996783842756842877E-11, ref b0, ref b1, ref b2);
	864	besselm1nextcheb(y, 9.21831518760500529508E-11, ref b0, ref b1, ref b2);
	865	besselm1nextcheb(y, -4.19035475934189648750E-10, ref b0, ref b1, ref b2);
	866	besselm1nextcheb(y, 2.01504975519703286596E-9, ref b0, ref b1, ref b2);
	867	besselm1nextcheb(y, -1.03457624656780970260E-8, ref b0, ref b1, ref b2);
	868	besselm1nextcheb(y, 5.74108412545004946722E-8, ref b0, ref b1, ref b2);
	869	besselm1nextcheb(y, -3.50196060308781257119E-7, ref b0, ref b1, ref b2);
	870	besselm1nextcheb(y, 2.40648494783721712015E-6, ref b0, ref b1, ref b2);
	871	besselm1nextcheb(y, -1.93619797416608296024E-5, ref b0, ref b1, ref b2);
	872	besselm1nextcheb(y, 1.95215518471351631108E-4, ref b0, ref b1, ref b2);
	873	besselm1nextcheb(y, -2.85781685962277938680E-3, ref b0, ref b1, ref b2);
	874	besselm1nextcheb(y, 1.03923736576817238437E-1, ref b0, ref b1, ref b2);
	875	besselm1nextcheb(y, 2.72062619048444266945E0, ref b0, ref b1, ref b2);
	876	v = 0.5*(b0-b2);
	877	result = Math.Exp(-x)*v/Math.Sqrt(x);
	878	}
	879	return result;
	880	}
	881
	882
	883	/*************************************************************************
	884	Modified Bessel function, second kind, integer order
	885
	886	Returns modified Bessel function of the second kind
	887	of order n of the argument.
	888
	889	The range is partitioned into the two intervals [0,9.55] and
	890	(9.55, infinity). An ascending power series is used in the
	891	low range, and an asymptotic expansion in the high range.
	892
	893	ACCURACY:
	894
	895	Relative error:
	896	arithmetic domain # trials peak rms
	897	IEEE 0,30 90000 1.8e-8 3.0e-10
	898
	899	Error is high only near the crossover point x = 9.55
	900	between the two expansions used.
	901
	902	Cephes Math Library Release 2.8: June, 2000
	903	Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
	904	*************************************************************************/
	905	public static double besselkn(int nn,
	906	double x)
	907	{
	908	double result = 0;
	909	double k = 0;
	910	double kf = 0;
	911	double nk1f = 0;
	912	double nkf = 0;
	913	double zn = 0;
	914	double t = 0;
	915	double s = 0;
	916	double z0 = 0;
	917	double z = 0;
	918	double ans = 0;
	919	double fn = 0;
	920	double pn = 0;
	921	double pk = 0;
	922	double zmn = 0;
	923	double tlg = 0;
	924	double tox = 0;
	925	int i = 0;
	926	int n = 0;
	927	double eul = 0;
	928
	929	eul = 5.772156649015328606065e-1;
	930	if( nn<0 )
	931	{
	932	n = -nn;
	933	}
	934	else
	935	{
	936	n = nn;
	937	}
	938	System.Diagnostics.Debug.Assert(n<=31, "Overflow in BesselKN");
	939	System.Diagnostics.Debug.Assert((double)(x)>(double)(0), "Domain error in BesselKN");
	940	if( (double)(x)<=(double)(9.55) )
	941	{
	942	ans = 0.0;
	943	z0 = 0.25xx;
	944	fn = 1.0;
	945	pn = 0.0;
	946	zmn = 1.0;
	947	tox = 2.0/x;
	948	if( n>0 )
	949	{
	950	pn = -eul;
	951	k = 1.0;
	952	for(i=1; i<=n-1; i++)
	953	{
	954	pn = pn+1.0/k;
	955	k = k+1.0;
	956	fn = fn*k;
	957	}
	958	zmn = tox;
	959	if( n==1 )
	960	{
	961	ans = 1.0/x;
	962	}
	963	else
	964	{
	965	nk1f = fn/n;
	966	kf = 1.0;
	967	s = nk1f;
	968	z = -z0;
	969	zn = 1.0;
	970	for(i=1; i<=n-1; i++)
	971	{
	972	nk1f = nk1f/(n-i);
	973	kf = kf*i;
	974	zn = zn*z;
	975	t = nk1f*zn/kf;
	976	s = s+t;
	977	System.Diagnostics.Debug.Assert((double)(AP.Math.MaxRealNumber-Math.Abs(t))>(double)(Math.Abs(s)), "Overflow in BesselKN");
	978	System.Diagnostics.Debug.Assert(!((double)(tox)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/tox)<(double)(zmn)), "Overflow in BesselKN");
	979	zmn = zmn*tox;
	980	}
	981	s = s*0.5;
	982	t = Math.Abs(s);
	983	System.Diagnostics.Debug.Assert(!((double)(zmn)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/zmn)<(double)(t)), "Overflow in BesselKN");
	984	System.Diagnostics.Debug.Assert(!((double)(t)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/t)<(double)(zmn)), "Overflow in BesselKN");
	985	ans = s*zmn;
	986	}
	987	}
	988	tlg = 2.0Math.Log(0.5x);
	989	pk = -eul;
	990	if( n==0 )
	991	{
	992	pn = pk;
	993	t = 1.0;
	994	}
	995	else
	996	{
	997	pn = pn+1.0/n;
	998	t = 1.0/fn;
	999	}
	1000	s = (pk+pn-tlg)*t;
	1001	k = 1.0;
	1002	do
	1003	{
	1004	t = t(z0/(k(k+n)));
	1005	pk = pk+1.0/k;
	1006	pn = pn+1.0/(k+n);
	1007	s = s+(pk+pn-tlg)*t;
	1008	k = k+1.0;
	1009	}
	1010	while( (double)(Math.Abs(t/s))>(double)(AP.Math.MachineEpsilon) );
	1011	s = 0.5*s/zmn;
	1012	if( n%2!=0 )
	1013	{
	1014	s = -s;
	1015	}
	1016	ans = ans+s;
	1017	result = ans;
	1018	return result;
	1019	}
	1020	if( (double)(x)>(double)(Math.Log(AP.Math.MaxRealNumber)) )
	1021	{
	1022	result = 0;
	1023	return result;
	1024	}
	1025	k = n;
	1026	pn = 4.0kk;
	1027	pk = 1.0;
	1028	z0 = 8.0*x;
	1029	fn = 1.0;
	1030	t = 1.0;
	1031	s = t;
	1032	nkf = AP.Math.MaxRealNumber;
	1033	i = 0;
	1034	do
	1035	{
	1036	z = pn-pk*pk;
	1037	t = tz/(fnz0);
	1038	nk1f = Math.Abs(t);
	1039	if( i>=n & (double)(nk1f)>(double)(nkf) )
	1040	{
	1041	break;
	1042	}
	1043	nkf = nk1f;
	1044	s = s+t;
	1045	fn = fn+1.0;
	1046	pk = pk+2.0;
	1047	i = i+1;
	1048	}
	1049	while( (double)(Math.Abs(t/s))>(double)(AP.Math.MachineEpsilon) );
	1050	result = Math.Exp(-x)Math.Sqrt(Math.PI/(2.0x))*s;
	1051	return result;
	1052	}
	1053
	1054
	1055	/*************************************************************************
	1056	Internal subroutine
	1057
	1058	Cephes Math Library Release 2.8: June, 2000
	1059	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	1060	*************************************************************************/
	1061	private static void besselmfirstcheb(double c,
	1062	ref double b0,
	1063	ref double b1,
	1064	ref double b2)
	1065	{
	1066	b0 = c;
	1067	b1 = 0.0;
	1068	b2 = 0.0;
	1069	}
	1070
	1071
	1072	/*************************************************************************
	1073	Internal subroutine
	1074
	1075	Cephes Math Library Release 2.8: June, 2000
	1076	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	1077	*************************************************************************/
	1078	private static void besselmnextcheb(double x,
	1079	double c,
	1080	ref double b0,
	1081	ref double b1,
	1082	ref double b2)
	1083	{
	1084	b2 = b1;
	1085	b1 = b0;
	1086	b0 = x*b1-b2+c;
	1087	}
	1088
	1089
	1090	/*************************************************************************
	1091	Internal subroutine
	1092
	1093	Cephes Math Library Release 2.8: June, 2000
	1094	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	1095	*************************************************************************/
	1096	private static void besselm1firstcheb(double c,
	1097	ref double b0,
	1098	ref double b1,
	1099	ref double b2)
	1100	{
	1101	b0 = c;
	1102	b1 = 0.0;
	1103	b2 = 0.0;
	1104	}
	1105
	1106
	1107	/*************************************************************************
	1108	Internal subroutine
	1109
	1110	Cephes Math Library Release 2.8: June, 2000
	1111	Copyright 1984, 1987, 2000 by Stephen L. Moshier
	1112	*************************************************************************/
	1113	private static void besselm1nextcheb(double x,
	1114	double c,
	1115	ref double b0,
	1116	ref double b1,
	1117	ref double b2)
	1118	{
	1119	b2 = b1;
	1120	b1 = b0;
	1121	b0 = x*b1-b2+c;
	1122	}
	1123
	1124
	1125	private static void besselasympt0(double x,
	1126	ref double pzero,
	1127	ref double qzero)
	1128	{
	1129	double xsq = 0;
	1130	double p2 = 0;
	1131	double q2 = 0;
	1132	double p3 = 0;
	1133	double q3 = 0;
	1134
	1135	xsq = 64.0/(x*x);
	1136	p2 = 0.0;
	1137	p2 = 2485.271928957404011288128951+xsq*p2;
	1138	p2 = 153982.6532623911470917825993+xsq*p2;
	1139	p2 = 2016135.283049983642487182349+xsq*p2;
	1140	p2 = 8413041.456550439208464315611+xsq*p2;
	1141	p2 = 12332384.76817638145232406055+xsq*p2;
	1142	p2 = 5393485.083869438325262122897+xsq*p2;
	1143	q2 = 1.0;
	1144	q2 = 2615.700736920839685159081813+xsq*q2;
	1145	q2 = 156001.7276940030940592769933+xsq*q2;
	1146	q2 = 2025066.801570134013891035236+xsq*q2;
	1147	q2 = 8426449.050629797331554404810+xsq*q2;
	1148	q2 = 12338310.22786324960844856182+xsq*q2;
	1149	q2 = 5393485.083869438325560444960+xsq*q2;
	1150	p3 = -0.0;
	1151	p3 = -4.887199395841261531199129300+xsq*p3;
	1152	p3 = -226.2630641933704113967255053+xsq*p3;
	1153	p3 = -2365.956170779108192723612816+xsq*p3;
	1154	p3 = -8239.066313485606568803548860+xsq*p3;
	1155	p3 = -10381.41698748464093880530341+xsq*p3;
	1156	p3 = -3984.617357595222463506790588+xsq*p3;
	1157	q3 = 1.0;
	1158	q3 = 408.7714673983499223402830260+xsq*q3;
	1159	q3 = 15704.89191515395519392882766+xsq*q3;
	1160	q3 = 156021.3206679291652539287109+xsq*q3;
	1161	q3 = 533291.3634216897168722255057+xsq*q3;
	1162	q3 = 666745.4239319826986004038103+xsq*q3;
	1163	q3 = 255015.5108860942382983170882+xsq*q3;
	1164	pzero = p2/q2;
	1165	qzero = 8*p3/q3/x;
	1166	}
	1167
	1168
	1169	private static void besselasympt1(double x,
	1170	ref double pzero,
	1171	ref double qzero)
	1172	{
	1173	double xsq = 0;
	1174	double p2 = 0;
	1175	double q2 = 0;
	1176	double p3 = 0;
	1177	double q3 = 0;
	1178
	1179	xsq = 64.0/(x*x);
	1180	p2 = -1611.616644324610116477412898;
	1181	p2 = -109824.0554345934672737413139+xsq*p2;
	1182	p2 = -1523529.351181137383255105722+xsq*p2;
	1183	p2 = -6603373.248364939109255245434+xsq*p2;
	1184	p2 = -9942246.505077641195658377899+xsq*p2;
	1185	p2 = -4435757.816794127857114720794+xsq*p2;
	1186	q2 = 1.0;
	1187	q2 = -1455.009440190496182453565068+xsq*q2;
	1188	q2 = -107263.8599110382011903063867+xsq*q2;
	1189	q2 = -1511809.506634160881644546358+xsq*q2;
	1190	q2 = -6585339.479723087072826915069+xsq*q2;
	1191	q2 = -9934124.389934585658967556309+xsq*q2;
	1192	q2 = -4435757.816794127856828016962+xsq*q2;
	1193	p3 = 35.26513384663603218592175580;
	1194	p3 = 1706.375429020768002061283546+xsq*p3;
	1195	p3 = 18494.26287322386679652009819+xsq*p3;
	1196	p3 = 66178.83658127083517939992166+xsq*p3;
	1197	p3 = 85145.16067533570196555001171+xsq*p3;
	1198	p3 = 33220.91340985722351859704442+xsq*p3;
	1199	q3 = 1.0;
	1200	q3 = 863.8367769604990967475517183+xsq*q3;
	1201	q3 = 37890.22974577220264142952256+xsq*q3;
	1202	q3 = 400294.4358226697511708610813+xsq*q3;
	1203	q3 = 1419460.669603720892855755253+xsq*q3;
	1204	q3 = 1819458.042243997298924553839+xsq*q3;
	1205	q3 = 708712.8194102874357377502472+xsq*q3;
	1206	pzero = p2/q2;
	1207	qzero = 8*p3/q3/x;
	1208	}
	1209	}
	1210	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences