airy.c

su 的源代码库

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

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

	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.0;
			ix = ax;
			frac = ax-ix;
			bip = (1.0-frac)*bpp[ix]+frac*bpp[ix+1];
		} else {
			sx = sqrt(x);
			y = 1.5/(x*sx);
			if (y>0.5) {
				ay = 15.0-y*10.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				g = (1.0-frac)*gp1[iy]+frac*gp1[iy+1];
				bip = exp(1.0/y)*g*sqrt(sx);
			} else {
				ay = 10.0-y*20.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				g = (1.0-frac)*gp2[iy]+frac*gp2[iy+1];
				bip = exp(1.0/y)*g*sqrt(sx);
			}
		}
	} else {
		if (x>-10.0) {
			if (x>-1.0) {
				ax = -x*100.0;
				ix = ax;
				frac = ax-ix;
				bip = (1.0-frac)*bpn1[ix]+frac*bpn1[ix+1];
			} else {
				ax = (-x-1.0)*10.0;
				ix = ax;
				frac = ax-ix;
				bip = (1.0-frac)*bpn2[ix]+frac*bpn2[ix+1];
			}
		} else {
			sx = sqrt(-x);
			y = 1.5/(-x*sx);
			ay = 5.0-y*100.0;
			iy = ay;
			frac = ay-iy;
			g1 = (1.0-frac)*gn1[iy]+frac*gn1[iy+1];
			g2 = (1.0-frac)*gn2[iy]+frac*gn2[iy+1];
			oy = 1.0/y;
			bip = (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