DataMuseum.dk

top - download

    Length: 2202 (0x89a)
    Types: TextFile
    Notes: UNIX file
    Names: »j1.c«

└─⟦f27320a65⟧ Bits:30001972 Commodore 900 hard disk image with partial source code
    └─⟦f4b8d8c84⟧ UNIX V7 Filesystem
        └─ ⟦this⟧ »libm/j1.c« 

TextFile

/*
 * Evaluate the Bessel function of the first kind for the
 * 0'th order.
 */
#include <math.h>

/*
 * (Hart 6046, 17.13)
 */
double	j1sntab[] ={
	 0.2214887880421963139207647803e+18,
	-0.2512374214703212789513276482e+17,
	 0.8482420744781272654092270714e+15,
	-0.1249820367262024853059386404e+14,
	 0.9316565296724673204941569090e+11,
	-0.3686668987022981626057360734e+09,
	 0.7437023817119996441033971568e+06,
	-0.6079530179607413599422162589e+03
};
double	j1smtab[] ={
	 0.4429775760843926213068606431e+18,
	 0.5124712716484872112355190833e+16,
	 0.2989836307725487159974863805e+14,
	 0.1158192127466889329393351964e+12,
	 0.3281940344534196444499723799e+09,
	 0.6988586184485075744033771749e+06,
	 0.1077741289433304357312995618e+04,
	 0.1000000000000000000000000000e+01
};

/*
 * (Hart 6748, 16.96)
 */
double	j1pntab[] ={
	-0.1746957669440928588993119874e+07,
	-0.2366943914052142834823993894e+07,
	-0.8285036373872377453364829379e+06,
	-0.7957456871350595920771321737e+05,
	-0.1427906628827098109925825113e+04
};
double	j1pmtab[] ={
	-0.1746957669440928570099951216e+07,
	-0.2363745139022653985015374700e+07,
	-0.8242369906562818850190419776e+06,
	-0.7814405008939111083484764129e+05,
	-0.1308452983339079687508155034e+04,
	 0.1000000000000000000000000000e+01
};

/*
 * (Hart 7148, 17.16)
 */
double	j1qntab[] ={
	 0.1437540512394930951685647580e+05,
	 0.2288127654437247527255486540e+05,
	 0.9758129726644802407557480600e+04,
	 0.1212097049668108451079609000e+04,
	 0.3134055238623227003160600000e+02
};
double	j1qmtab[] ={
	 0.3066753093109186472663730271e+06,
	 0.4894441578927938446024906527e+06,
	 0.2102091447442188749228749142e+06,
	 0.2667395055960220714890160450e+05,
	 0.7531715785291537658631247000e+03,
	 0.1000000000000000000000000000e+01
};

double
j1(x)
double x;
{
	double r, p, q, z, xn;
	register int m;

	m = 0;
	if (x < 0.0) {
		x = -x;
		m = 1;
	}
	if (x <= 8.0) {
		r = x*x;
		r = x*_pol(r, j1sntab, 8)/_pol(r, j1smtab, 8);
	} else {
		z = 8.0 / x;
		xn = x - (3.0*PI)/4.0;
		r = x*x;
		p = _pol(r, j1pntab, 5) / _pol(r, j1pmtab, 6);
		q = z*_pol(r, j1qntab, 5) / _pol(r, j1qmtab, 6);
		r = sqrt(2.0/(PI*x)) * (p*cos(xn) - q*sin(xn));
	}
	if (m)
		r = -r;
	return (r);
}