airy.c

来自「该程序是用vc开发的对动态数组进行管理的DLL」· C语言 代码 · 共 647 行 · 第 1/2 页

			frac = ay-iy;
			g1 = (1.0f-frac)*gn1[iy]+frac*gn1[iy+1];
			g2 = (1.0f-frac)*gn2[iy]+frac*gn2[iy+1];
			oy = 1.0f/y;
			aip = (float)((g1*sin(oy)-g2*cos(oy))*sqrt(sx));
		}
	}
	return aip;
}


float airyb (float x)
/*****************************************************************************
Return approximation to the Airy function Bi(x)
******************************************************************************
Input:
x		value at which to evaluate Bi(x)

Returned:	Bi(x)
******************************************************************************
Notes:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 07/05/89
*****************************************************************************/
{
	static float bp[101] = {
		0.61492f, 0.61940f, 0.62389f, 0.62837f, 0.63286f,
		0.63735f, 0.64184f, 0.64634f, 0.65084f, 0.65534f,
		0.65986f, 0.66438f, 0.66890f, 0.67343f, 0.67798f,
		0.68253f, 0.68709f, 0.69167f, 0.69625f, 0.70085f,
		0.70546f, 0.71008f, 0.71472f, 0.71938f, 0.72405f,
		0.72874f, 0.73345f, 0.73818f, 0.74292f, 0.74769f,
		0.75248f, 0.75729f, 0.76213f, 0.76699f, 0.77187f,
		0.77678f, 0.78172f, 0.78669f, 0.79169f, 0.79671f,
		0.80177f, 0.80686f, 0.81198f, 0.81714f, 0.82233f,
		0.82755f, 0.83282f, 0.83812f, 0.84346f, 0.84885f,
		0.85427f, 0.85974f, 0.86525f, 0.87081f, 0.87641f,
		0.88206f, 0.88776f, 0.89350f, 0.89930f, 0.90515f,
		0.91106f, 0.91702f, 0.92303f, 0.92910f, 0.93524f,
		0.94143f, 0.94768f, 0.95399f, 0.96037f, 0.96681f,
		0.97332f, 0.97990f, 0.98655f, 0.99327f, 1.00006f,
		1.00693f, 1.01387f, 1.02088f, 1.02798f, 1.03516f,
		1.04242f, 1.04976f, 1.05719f, 1.06470f, 1.07230f,
		1.07999f, 1.08778f, 1.09566f, 1.10363f, 1.11170f,
		1.11987f, 1.12814f, 1.13651f, 1.14499f, 1.15357f,
		1.16226f, 1.17106f, 1.17998f, 1.18901f, 1.19815f, 1.20742f };
	static float fp1[11] = {
		0.619912f, 0.620335f, 0.620327f, 0.619799f, 0.618649f,
		0.616764f, 0.614022f, 0.610309f, 0.605543f, 0.599723f, 0.593015f };
	static float fp2[11] = {
		0.593015f, 0.589451f, 0.585855f, 0.582330f, 0.578985f,
		0.575908f, 0.573135f, 0.570636f, 0.568343f, 0.566204f, 0.564190f };
	static float bn1[101] = {
		0.61492f, 0.61044f, 0.60596f, 0.60147f, 0.59698f,
		0.59249f, 0.58800f, 0.58351f, 0.57901f, 0.57450f,
		0.56999f, 0.56548f, 0.56096f, 0.55643f, 0.55189f,
		0.54735f, 0.54280f, 0.53824f, 0.53367f, 0.52909f,
		0.52450f, 0.51990f, 0.51529f, 0.51067f, 0.50604f,
		0.50139f, 0.49674f, 0.49207f, 0.48738f, 0.48268f,
		0.47797f, 0.47325f, 0.46851f, 0.46375f, 0.45898f,
		0.45419f, 0.44939f, 0.44457f, 0.43974f, 0.43488f,
		0.43002f, 0.42513f, 0.42023f, 0.41531f, 0.41037f,
		0.40541f, 0.40043f, 0.39544f, 0.39043f, 0.38540f,
		0.38035f, 0.37528f, 0.37019f, 0.36508f, 0.35996f,
		0.35481f, 0.34965f, 0.34446f, 0.33936f, 0.33403f,
		0.32879f, 0.32352f, 0.31824f, 0.31294f, 0.30761f,
		0.30227f, 0.29691f, 0.29153f, 0.28612f, 0.28070f,
		0.27526f, 0.26980f, 0.26432f, 0.25883f, 0.25331f,
		0.24777f, 0.24222f, 0.23665f, 0.23106f, 0.22545f,
		0.21982f, 0.21418f, 0.20852f, 0.20284f, 0.19714f,
		0.19143f, 0.18570f, 0.17996f, 0.17420f, 0.16842f,
		0.16263f, 0.15683f, 0.15101f, 0.14518f, 0.13933f,
		0.13347f, 0.12760f, 0.12171f, 0.11582f, 0.10991f, 0.10399f };
	static float bn2[91] = {
		 0.10399f, 0.04432f,-0.01582f,-0.07576f,-0.13472f,
		-0.19178f,-0.24596f,-0.29620f,-0.34140f,-0.38046f,
		-0.41230f,-0.43590f,-0.45036f,-0.45492f,-0.44905f,
		-0.43242f,-0.40500f,-0.36709f,-0.31929f,-0.26258f,
		-0.19828f,-0.12807f,-0.05390f, 0.02196f, 0.09710f,
		 0.16893f, 0.23486f, 0.29235f, 0.33904f, 0.37289f,
		 0.39223f, 0.39593f, 0.38346f, 0.35494f, 0.31122f,
		 0.25387f, 0.18514f, 0.10794f, 0.02570f,-0.05774f,
		-0.13836f,-0.21208f,-0.27502f,-0.32371f,-0.35531f,
		-0.36781f,-0.36017f,-0.33245f,-0.28589f,-0.22282f,
		-0.14669f,-0.06182f, 0.02679f, 0.11373f, 0.19354f,
		 0.26101f, 0.31159f, 0.34172f, 0.34908f, 0.33283f,
		 0.29376f, 0.23425f, 0.15821f, 0.07087f,-0.02159f,
		-0.11246f,-0.19493f,-0.26267f,-0.31030f,-0.33387f,
		-0.33125f,-0.30230f,-0.24904f,-0.17550f,-0.08751f,
		 0.00775f, 0.10235f, 0.18820f, 0.25778f, 0.30483f,
		 0.32494f, 0.31603f, 0.27858f, 0.21570f, 0.13293f,
		 0.03778f,-0.06091f,-0.15379f,-0.23186f,-0.28738f,-0.31467f };
	static float fn1[6] = {
		0.39752f, 0.39781f, 0.39809f, 0.39838f, 0.39866f, 0.39894f };
	static float fn2[6] = {
		0.40028f, 0.40002f, 0.39975f, 0.39949f, 0.39921f, 0.39894f };

	int ix,iy;
	float ax,frac,bi,sx,y,ay,f,f1,f2,oy;

	if (x>=0.0) {
		if (x<1.0) {
			ax = x*100.0f;
			ix = (int)ax;
			frac = ax-ix;
			bi = (1.0f-frac)*bp[ix]+frac*bp[ix+1];
		} else {
			sx = (float)sqrt(x);
			y = 1.5f/(x*sx);
			if (y>0.5) {
				ay = 15.0f-y*10.0f;
				iy = (int)ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				f = (1.0f-frac)*fp1[iy]+frac*fp1[iy+1];
				bi = (float)(exp(1.0/y)*f/sqrt(sx));
			} else {
				ay = 10.0f-y*20.0f;
				iy = (int)ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				f = (1.0f-frac)*fp2[iy]+frac*fp2[iy+1];
				bi = (float)(exp(1.0/y)*f/sqrt(sx));
			}
		}
	} else {
		if (x>-10.0) {
			if (x>-1.0) {
				ax = -x*100.0f;
				ix = (int)ax;
				frac = ax-ix;
				bi = (1.0f-frac)*bn1[ix]+frac*bn1[ix+1];
			} else {
				ax = (-x-1.0f)*10.0f;
				ix = (int)ax;
				frac = ax-ix;
				bi = (1.0f-frac)*bn2[ix]+frac*bn2[ix+1];
			}
		} else {
			sx = (float)(sqrt(-x));
			y = 1.5f/(-x*sx);
			ay = 5.0f-y*100.0f;
			iy = (int)ay;
			frac = ay-iy;
			f1 = (1.0f-frac)*fn1[iy]+frac*fn1[iy+1];
			f2 = (1.0f-frac)*fn2[iy]+frac*fn2[iy+1];
			oy = 1.0f/y;
			bi = (float)((f2*cos(oy)-f1*sin(oy))/sqrt(sx));
		}
	}
	return bi;
}


float airybp (float x)
/*****************************************************************************
Return approximation to the derivative of the Airy function Bi'(x)
******************************************************************************
Input:
x		value at which to evaluate Bi'(x)

Returned:	Bi'(x)
******************************************************************************
Notes:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.
******************************************************************************
Authors:  Dave Hale, Craig Artley, Colorado School of Mines, 11/13/90
*****************************************************************************/
{
	static float bpp[101] = {
		0.44828f, 0.44831f, 0.44841f, 0.44856f, 0.44878f,
		0.44907f, 0.44942f, 0.44984f, 0.45033f, 0.45088f,
		0.45151f, 0.45220f, 0.45297f, 0.45381f, 0.45472f,
		0.45571f, 0.45677f, 0.45791f, 0.45912f, 0.46041f,
		0.46178f, 0.46324f, 0.46477f, 0.46638f, 0.46808f,
		0.46986f, 0.47172f, 0.47367f, 0.47571f, 0.47783f,
		0.48004f, 0.48235f, 0.48474f, 0.48722f, 0.48980f,
		0.49247f, 0.49524f, 0.49810f, 0.50106f, 0.50412f,
		0.50728f, 0.51053f, 0.51389f, 0.51736f, 0.52092f,
		0.52459f, 0.52837f, 0.53225f, 0.53625f, 0.54035f,
		0.54457f, 0.54890f, 0.55334f, 0.55789f, 0.56257f,
		0.56736f, 0.57227f, 0.57730f, 0.58246f, 0.58774f,
		0.59314f, 0.59867f, 0.60433f, 0.61012f, 0.61603f,
		0.62209f, 0.62827f, 0.63460f, 0.64106f, 0.64766f,
		0.65440f, 0.66129f, 0.66832f, 0.67549f, 0.68282f,
		0.69029f, 0.69792f, 0.70571f, 0.71365f, 0.72174f,
		0.73000f, 0.73842f, 0.74701f, 0.75576f, 0.76468f,
		0.77378f, 0.78304f, 0.79249f, 0.80211f, 0.81191f,
		0.82190f, 0.83207f, 0.84243f, 0.85298f, 0.86373f,
		0.87467f, 0.88581f, 0.89716f, 0.90871f, 0.92046f, 0.93243f };
	static float gp1[11] = {
		0.478728f, 0.479925f, 0.481658f, 0.484018f, 0.487107f,
		0.491037f, 0.495921f, 0.501859f, 0.508909f, 0.517032f, 0.526011f };
	static float gp2[11] = {
		0.526011f, 0.530678f, 0.535345f, 0.539902f, 0.544235f,
		0.548255f, 0.551930f, 0.555296f, 0.558428f, 0.561382f, 0.564190f };
	static float bpn1[101] = {
		0.44828f, 0.44831f, 0.44841f, 0.44856f, 0.44877f,
		0.44903f, 0.44936f, 0.44974f, 0.45017f, 0.45066f,
		0.45121f, 0.45180f, 0.45245f, 0.45315f, 0.45390f,
		0.45470f, 0.45554f, 0.45643f, 0.45737f, 0.45835f,
		0.45938f, 0.46045f, 0.46156f, 0.46272f, 0.46391f,
		0.46515f, 0.46642f, 0.46773f, 0.46908f, 0.47046f,
		0.47188f, 0.47333f, 0.47481f, 0.47632f, 0.47787f,
		0.47944f, 0.48105f, 0.48268f, 0.48434f, 0.48602f,
		0.48773f, 0.48946f, 0.49122f, 0.49299f, 0.49479f,
		0.49660f, 0.49844f, 0.50029f, 0.50215f, 0.50403f,
		0.50593f, 0.50784f, 0.50976f, 0.51169f, 0.51363f,
		0.51557f, 0.51753f, 0.51949f, 0.52145f, 0.52342f,
		0.52540f, 0.52737f, 0.52934f, 0.53132f, 0.53329f,
		0.53525f, 0.53721f, 0.53917f, 0.54112f, 0.54306f,
		0.54499f, 0.54692f, 0.54883f, 0.55072f, 0.55260f,
		0.55447f, 0.55632f, 0.55815f, 0.55996f, 0.56176f,
		0.56353f, 0.56527f, 0.56699f, 0.56869f, 0.57036f,
		0.57200f, 0.57362f, 0.57520f, 0.57675f, 0.57826f,
		0.57974f, 0.58119f, 0.58260f, 0.58397f, 0.58530f,
		0.58659f, 0.58783f, 0.58904f, 0.59019f, 0.59131f, 0.59237f };
	static float bpn2[91] = {
		 0.59237f, 0.60011f, 0.60171f, 0.59592f, 0.58165f,
		 0.55790f, 0.52389f, 0.47906f, 0.42315f, 0.35624f,
		 0.27879f, 0.19168f, 0.09622f,-0.00581f,-0.11223f,
		-0.22042f,-0.32739f,-0.42989f,-0.52445f,-0.60751f,
		-0.67561f,-0.72544f,-0.75412f,-0.75926f,-0.73920f,
		-0.69311f,-0.62117f,-0.52461f,-0.40581f,-0.26829f,
		-0.11667f, 0.04347f, 0.20575f, 0.36320f, 0.50858f,
		 0.63474f, 0.73494f, 0.80328f, 0.83508f, 0.82721f,
		 0.77841f, 0.68948f, 0.56345f, 0.40555f, 0.22307f,
		 0.02511f,-0.17783f,-0.37440f,-0.55300f,-0.70247f,
		-0.81289f,-0.87622f,-0.88697f,-0.84276f,-0.74461f,
		-0.59717f,-0.40856f,-0.19009f, 0.04437f, 0.27926f,
		 0.49824f, 0.68542f, 0.82650f, 0.90998f, 0.92812f,
		 0.87780f, 0.76095f, 0.58474f, 0.36122f, 0.10670f,
		-0.15945f,-0.41615f,-0.64232f,-0.81860f,-0.92910f,
		-0.96296f,-0.91547f,-0.78882f,-0.59221f,-0.34136f,
		-0.05740f, 0.23484f, 0.50894f, 0.73928f, 0.90348f,
		 0.98471f, 0.97349f, 0.86898f, 0.67936f, 0.42147f, 0.11941f };
	static float gn1[6] = {
		0.40092f, 0.40052f, 0.40012f, 0.39972f, 0.39933f, 0.39894f };
	static float gn2[6] = {
		0.39704f, 0.39741f, 0.39779f, 0.39817f, 0.39855f, 0.39894f };

	int ix,iy;
	float ax,frac,bip,sx,y,ay,g,g1,g2,oy;

	if (x>=0.0) {
		if (x<1.0) {
			ax = x*100.0f;
			ix = (int)ax;
			frac = ax-ix;
			bip = (1.0f-frac)*bpp[ix]+frac*bpp[ix+1];
		} else {
			sx = (float)sqrt(x);
			y = 1.5f/(x*sx);
			if (y>0.5) {
				ay = 15.0f-y*10.0f;
				iy = (int)ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				g = (1.0f-frac)*gp1[iy]+frac*gp1[iy+1];
				bip = (float)(exp(1.0/y)*g*sqrt(sx));
			} else {
				ay = 10.0f-y*20.0f;
				iy = (int)ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				g = (1.0f-frac)*gp2[iy]+frac*gp2[iy+1];
				bip = (float)(exp(1.0/y)*g*sqrt(sx));
			}
		}
	} else {
		if (x>-10.0) {
			if (x>-1.0) {
				ax = -x*100.0f;
				ix = (int)ax;
				frac = ax-ix;
				bip = (1.0f-frac)*bpn1[ix]+frac*bpn1[ix+1];
			} else {
				ax = (-x-1.0f)*10.0f;
				ix = (int)ax;
				frac = ax-ix;
				bip = (1.0f-frac)*bpn2[ix]+frac*bpn2[ix+1];
			}
		} else {
			sx = (float)sqrt(-x);
			y = 1.5f/(-x*sx);
			ay = 5.0f-y*100.0f;
			iy = (int)ay;
			frac = ay-iy;
			g1 = (1.0f-frac)*gn1[iy]+frac*gn1[iy+1];
			g2 = (1.0f-frac)*gn2[iy]+frac*gn2[iy+1];
			oy = 1.0f/y;
			bip = (float)((g1*cos(oy)+g2*sin(oy))*sqrt(sx));
		}
	}
	return bip;
}

#ifdef TEST
main()
{
	int nx,i;
	float xmin,xmax,x,dx;

	scanf("%d %*[^\n]", &nx);
	scanf("%f %*[^\n]", &xmin);
	scanf("%f %*[^\n]", &xmax);

	for (i=0,x=xmin,dx=(xmax-xmin)/(nx-1); i<nx; i++,x+=dx)
		printf("ai[%f] = %f\n", x, airya(x));
		printf("aip[%f] = %f\n", x, airyap(x));
		printf("bi[%f] = %f\n", x, airyb(x));
		printf("bip[%f] = %f\n", x, airybp(x));
}
#endif