#include <math.h>

void bessy01(float x, float *y0, float *y1)
{
	if (x < 8.0) {
		int i;
		float bessj0(float);
		float bessj1(float);
		float z,z2,c,lnx,b0,b1,b2;
		static float ar1[15]={0.164349e-14, -0.8747341e-13,
			0.402633082e-11, -0.15837552542e-9, 0.524879478733e-8,
			-0.14407233274019e-6, 0.32065325376548e-5,
			-0.563207914105699e-4, 0.753113593257774e-3,
			-0.72879624795521e-2, 0.471966895957634e-1,
			-0.177302012781143, 0.261567346255047,
			0.179034314077182, -0.274474305529745};
		static float ar2[15]={0.42773e-15, -0.2440949e-13,
			0.121143321e-11, -0.5172121473e-10, 0.187547032473e-8,
			-0.5688440039919e-7, 0.141662436449235e-5,
			-0.283046401495148e-4, 0.440478629867099e-3,
			-0.51316411610611e-2, 0.423191803533369e-1,
			-0.226624991556755, 0.675615780772188,
			-0.767296362886646, -0.128697384381350};
		c=0.636619772367581;
		lnx=c*log(x);
		c /= x;
		x /= 8.0;
		z=2.0*x*x-1.0;
		z2=z+z;
		b1=b2=0.0;
		for (i=0; i<=14; i++) {
			b0=z2*b1-b2+ar1[i];
			b2=b1;
			b1=b0;
		}
		*y0 = lnx*bessj0(8.0*x)+z*b1-b2-0.33146113203285e-1;
		b1=b2=0.0;
		for (i=0; i<=14; i++) {
			b0=z2*b1-b2+ar2[i];
			b2=b1;
			b1=b0;
		}
		*y1 = lnx*bessj1(8.0*x)-c+x*(z*b1-b2+0.2030410588593425e-1);
	} else {
		void besspq0(float, float *, float *);
		void besspq1(float, float *, float *);
		float c,cosx,sinx,p0,q0,p1,q1;
		c=0.797884560802865/sqrt(x);
		besspq0(x,&p0,&q0);
		besspq1(x,&p1,&q1);
		x -= 0.706858347057703e1;
		cosx=cos(x);
		sinx=sin(x);
		*y0 = c*(p0*sinx+q0*cosx);
		*y1 = c*(q1*sinx-p1*cosx);
	}
}