airy.c

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

/* Copyright (c) Colorado School of Mines, 2003.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/*****************************************************************************
AIRY - Approximate the Airy functions  Ai(x), Bi(x) and their respective
       derivatives  Ai'(x), Bi'(x)

airya		return approximation of Ai(x)
airypa		return approximation of Ai'(x)
airyb		return approximation of Bi(x)
airybp		return approximation of Bi'(x)

******************************************************************************
Function Prototypes:
float airya (float x);
float airyap (float x);
float airyb (float x);
float airybp (float x);

******************************************************************************
Input:
x		value at which to evaluate Ai(x)

Returned:
airya		Ai(x)
airypa		Ai'(x)
airyb		Bi(x)
airybp		Bi'(x)

******************************************************************************
Reference:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/06/89
*****************************************************************************/
/**************** end self doc ********************************/
#include "cwp.h"

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

Returned:	Ai(x)
******************************************************************************
Reference:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/06/89
*****************************************************************************/
{
	static float ap[101] = {
		0.35502f, 0.35243f, 0.34985f, 0.34726f, 0.34467f,
		0.34209f, 0.33951f, 0.33693f, 0.33435f, 0.33177f,
		0.32920f, 0.32663f, 0.32406f, 0.32150f, 0.31894f,
		0.31639f, 0.31384f, 0.31130f, 0.30876f, 0.30623f,
		0.30370f, 0.30118f, 0.29866f, 0.29615f, 0.29365f,
		0.29116f, 0.28867f, 0.28619f, 0.28372f, 0.28126f,
		0.27880f, 0.27635f, 0.27392f, 0.27149f, 0.26906f,
		0.26665f, 0.26425f, 0.26186f, 0.25947f, 0.25710f,
		0.25474f, 0.25238f, 0.25004f, 0.24771f, 0.24539f,
		0.24308f, 0.24078f, 0.23849f, 0.23621f, 0.23394f,
		0.23169f, 0.22945f, 0.22721f, 0.22499f, 0.22279f,
		0.22059f, 0.21841f, 0.21624f, 0.21408f, 0.21193f,
		0.20980f, 0.20767f, 0.20556f, 0.20347f, 0.20138f,
		0.19931f, 0.19726f, 0.19521f, 0.19318f, 0.19116f,
		0.18916f, 0.18717f, 0.18519f, 0.18322f, 0.18127f,
		0.17933f, 0.17741f, 0.17549f, 0.17360f, 0.17171f,
		0.16984f, 0.16798f, 0.16614f, 0.16431f, 0.16249f,
		0.16069f, 0.15890f, 0.15713f, 0.15536f, 0.15362f,
		0.15188f, 0.15016f, 0.14845f, 0.14676f, 0.14508f,
		0.14341f, 0.14176f, 0.14012f, 0.13850f, 0.13689f, 0.13529f };
	static float fp1[11] = {
		0.527027f, 0.528783f, 0.530601f, 0.532488f, 0.534448f,
		0.536489f, 0.538618f, 0.540844f, 0.543180f, 0.545636f, 0.548230f};
	static float fp2[11] = {
		0.548230f, 0.549584f, 0.550980f, 0.552421f, 0.553912f,
		0.555456f, 0.557058f, 0.558724f, 0.560462f, 0.562280f, 0.564190f };
	static float an1[101] = {
		0.35502f, 0.35761f, 0.36020f, 0.36279f, 0.36537f,
		0.36796f, 0.37054f, 0.37312f, 0.37560f, 0.37827f,
		0.38084f, 0.38341f, 0.38597f, 0.38853f, 0.39109f,
		0.39364f, 0.39618f, 0.39871f, 0.40124f, 0.40376f,
		0.40628f, 0.40879f, 0.41128f, 0.41377f, 0.41625f,
		0.41872f, 0.42118f, 0.42363f, 0.42606f, 0.42849f,
		0.43090f, 0.43330f, 0.43568f, 0.43805f, 0.44041f,
		0.44275f, 0.44508f, 0.44739f, 0.44968f, 0.45196f,
		0.45422f, 0.45646f, 0.45868f, 0.46089f, 0.46307f,
		0.46523f, 0.46738f, 0.46950f, 0.47159f, 0.47367f,
		0.47572f, 0.47775f, 0.47976f, 0.48174f, 0.48369f,
		0.48562f, 0.48752f, 0.48939f, 0.49124f, 0.49306f,
		0.49484f, 0.49660f, 0.49833f, 0.50003f, 0.50170f,
		0.50333f, 0.50493f, 0.50650f, 0.50803f, 0.50953f,
		0.51100f, 0.51242f, 0.51382f, 0.51517f, 0.51649f,
		0.51777f, 0.51901f, 0.52021f, 0.52137f, 0.52249f,
		0.52357f, 0.52461f, 0.52560f, 0.52655f, 0.52746f,
		0.52832f, 0.52914f, 0.52991f, 0.53064f, 0.53132f,
		0.53195f, 0.53254f, 0.53308f, 0.53356f, 0.53400f,
		0.53439f, 0.53473f, 0.53501f, 0.53525f, 0.53543f, 0.53556f };
	static float an2[91] = {
		 0.53556f, 0.53381f, 0.52619f, 0.51227f, 0.49170f,
		 0.46425f, 0.42986f, 0.38860f, 0.34076f, 0.28680f,
		 0.22740f, 0.16348f, 0.09614f, 0.02670f,-0.04333f,
		-0.11232f,-0.17850f,-0.24003f,-0.29509f,-0.34190f,
		-0.37881f,-0.40438f,-0.41744f,-0.41718f,-0.40319f,
		-0.37553f,-0.33477f,-0.28201f,-0.21885f,-0.14741f,
		-0.07026f, 0.00967f, 0.08921f, 0.16499f, 0.23370f,
		 0.29215f, 0.33749f, 0.36736f, 0.38003f, 0.37453f,
		 0.35076f, 0.30952f, 0.25258f, 0.18256f, 0.10293f,
		 0.01778f,-0.06833f,-0.15062f,-0.22435f,-0.28512f,
		-0.32914f,-0.35351f,-0.35652f,-0.33734f,-0.29713f,
		-0.23802f,-0.16352f,-0.07831f, 0.01210f, 0.10168f,
		 0.18428f, 0.25403f, 0.30585f, 0.33577f, 0.34132f,
		 0.32177f, 0.27825f, 0.21372f, 0.13285f, 0.04170f,
		-0.05270f,-0.14290f,-0.22159f,-0.28223f,-0.31959f,
		-0.33029f,-0.31311f,-0.26920f,-0.20205f,-0.11726f,
		-0.02213f, 0.07495f, 0.16526f, 0.24047f, 0.29347f,
		 0.31910f, 0.31465f, 0.28023f, 0.21886f, 0.13623f, 0.04024f};
	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,ai,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;
			ai = (1.0f-frac)*ap[ix]+frac*ap[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];
				ai = (float)(0.5*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];
				ai = (float)(0.5*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;
				ai = (1.0f-frac)*an1[ix]+frac*an1[ix+1];
			} else {
				ax = (-x-1.0f)*10.0f;
				ix = (int)ax;
				frac = ax-ix;
				ai = (1.0f-frac)*an2[ix]+frac*an2[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;
			ai = (float)((f1*cos(oy)+f2*sin(oy))/sqrt(sx));
		}
	}
	return ai;
}


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

Returned:	Ai'(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 app[101] = {
		-0.25881F, -0.25880F, -0.25874F, -0.25866F, -0.25854F,
		-0.25838F, -0.25819F, -0.25797F, -0.25772F, -0.25744F,
		-0.25713F, -0.25678F, -0.25641F, -0.25600F, -0.25557F,
		-0.25511F, -0.25462F, -0.25411F, -0.25356F, -0.25300F,
		-0.25240F, -0.25178F, -0.25114F, -0.25047F, -0.24977F,
		-0.24906F, -0.24832F, -0.24756F, -0.24677F, -0.24597F,
		-0.24514F, -0.24429F, -0.24343F, -0.24254F, -0.24164F,
		-0.24071F, -0.23977F, -0.23881F, -0.23783F, -0.23684F,
		-0.23583F, -0.23480F, -0.23376F, -0.23270F, -0.23163F,
		-0.23054F, -0.22944F, -0.22833F, -0.22720F, -0.22606F,
		-0.22491F, -0.22374F, -0.22257F, -0.22138F, -0.22018F,
		-0.21897F, -0.21775F, -0.21653F, -0.21529F, -0.21404F,
		-0.21279F, -0.21153F, -0.21025F, -0.20898F, -0.20769F,
		-0.20640F, -0.20510F, -0.20380F, -0.20248F, -0.20117F,
		-0.19985F, -0.19852F, -0.19719F, -0.19585F, -0.19451F,
		-0.19317F, -0.19182F, -0.19047F, -0.18912F, -0.18777F,
		-0.18641f, -0.18505f, -0.18369f, -0.18232f, -0.18096f,
		-0.17959f, -0.17823f, -0.17686f, -0.17549f, -0.17413f,
		-0.17276f, -0.17139f, -0.17003f, -0.16866f, -0.16730f,
		-0.16593f, -0.16457f, -0.16321f, -0.16185f, -0.16050f, -0.15914f };
	static float gp1[11] = {
		0.619954f, 0.617156f, 0.614275f, 0.611305f, 0.608239f,
		0.605068f, 0.601782f, 0.598372f, 0.594823f, 0.591120f, 0.587245f };
	static float gp2[11] = {
		0.587245f, 0.585235f, 0.583174f, 0.581056f, 0.578878f,
		0.576635f, 0.574320f, 0.571927f, 0.569448f, 0.566873f, 0.564190f };
	static float apn1[101] = {
		-0.25881f, -0.25880f, -0.25874f, -0.25865f, -0.25852f,
		-0.25836f, -0.25816f, -0.25792f, -0.25763f, -0.25731f,
		-0.25695f, -0.25655f, -0.25661f, -0.25563f, -0.25510f,
		-0.25453f, -0.25392f, -0.25326f, -0.25256f, -0.25182f,
		-0.25103f, -0.25019f, -0.24931f, -0.24838f, -0.24741f,
		-0.24638f, -0.24531f, -0.24419f, -0.24303f, -0.24181f,
		-0.24054f, -0.23922f, -0.23785f, -0.23643f, -0.23496f,
		-0.23344f, -0.23186f, -0.23023f, -0.22855f, -0.22682f,
		-0.22503f, -0.22318f, -0.22128f, -0.21933f, -0.21732f,
		-0.21525f, -0.21313f, -0.21095f, -0.20872f, -0.20643f,
		-0.20408f, -0.20167f, -0.19920f, -0.19668f, -0.19410f,
		-0.19146f, -0.18875f, -0.18600f, -0.18318f, -0.18030f,
		-0.17736f, -0.17436f, -0.17130f, -0.16818f, -0.16500f,
		-0.16176f, -0.15846f, -0.15509f, -0.15167f, -0.14818f,
		-0.14464f, -0.14103f, -0.13736f, -0.13363f, -0.12984f,
		-0.12599f, -0.12207f, -0.11810f, -0.11406f, -0.10996f,
		-0.10580f, -0.10159f, -0.09731f, -0.09297f, -0.08857f,
		-0.08410f, -0.07958f, -0.07500f, -0.07036f, -0.06566f,
		-0.06091f, -0.05609f, -0.05121f, -0.04628f, -0.04129f,
		-0.03624f, -0.03114f, -0.02597f, -0.02076f, -0.01548f, -0.01016f };
	static float apn2[91] = {
		-0.01016f, 0.04602f, 0.10703f, 0.17199f, 0.23981f,
		 0.30918f, 0.37854f, 0.44612f, 0.50999f, 0.56809f,
		 0.61825f, 0.65834f, 0.68624f, 0.70003f, 0.69801f,
		 0.67885f, 0.64163f, 0.58600f, 0.51221f, 0.42118f,
		 0.31458f, 0.19482f, 0.06503f,-0.07096f,-0.20874f,
		-0.34344f,-0.46986f,-0.58272f,-0.67688f,-0.74755f,
		-0.79062f,-0.80287f,-0.78221f,-0.72794f,-0.64085f,
		-0.52336f,-0.37953f,-0.21499f,-0.03676f, 0.14695f,
		 0.32719f, 0.49458f, 0.63990f, 0.75457f, 0.83122f,
		 0.86419f, 0.85003f, 0.78781f, 0.67943f, 0.52962f,
		 0.34593f, 0.13836f,-0.08106f,-0.29899f,-0.50147f,
		-0.67495f,-0.80711f,-0.88790f,-0.91030f,-0.87103f,
		-0.77100f,-0.61552f,-0.41412f,-0.18009f, 0.07027f,
		 0.31880f, 0.54671f, 0.73605f, 0.87115f, 0.94004f,
		 0.93556f, 0.85621f, 0.70659f, 0.49727f, 0.24422f,
		-0.03231f,-0.30933f,-0.56297f,-0.77061f,-0.91289f,
		-0.97566f,-0.95149f,-0.84067f,-0.65149f,-0.39986f,
		-0.10809f, 0.19695f, 0.48628f, 0.73154f, 0.90781f, 0.99626f };
	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,aip,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;
			aip = (1.0f-frac)*app[ix]+frac*app[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];
				aip = (float)(-0.5f*exp(-1.0f/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];
				aip = (float)(-0.5*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;
				aip = (1.0f-frac)*apn1[ix]+frac*apn1[ix+1];
			} else {
				ax = (-x-1.0f)*10.0f;
				ix = (int)ax;
				frac = ax-ix;
				aip = (1.0f-frac)*apn2[ix]+frac*apn2[ix+1];
			}
		} else {
			sx = (float)sqrt(-x);
			y = 1.5f/(-x*sx);
			ay = 5.0f-y*100.0f;
			iy = (int)ay;