airyb.c

seismic software,very useful

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

/*
FUNCTION:  return approximation to the Airy function Bi(x)

PARAMETERS:
x			i value at which to evaluate 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
*/

#include <math.h>

float airyb (float x)
{
	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;
}