#include <math.h> float bessj1(float x) { if (x == 0.0) return 1.0; if (fabs(x) < 8.0) { int i; float z,z2,b0,b1,b2; static float ar[15]={-0.19554e-15, 0.1138572e-13, -0.57774042e-12, 0.2528123664e-10, -0.94242129816e-9, 0.2949707007278e-7, -0.76175878054003e-6, 0.158870192399321e-4, -0.260444389348581e-3, 0.324027018268386e-2, -0.291755248061542e-1, 0.177709117239728e0, -0.661443934134543e0, 0.128799409885768e1, -0.119180116054122e1}; 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+ar[i]; b2=b1; b1=b0; } return x*(z*b1-b2+0.648358770605265); } else { void besspq1(float, float *, float *); int sgnx; float c,cosx,sinx,p1,q1; sgnx = (x > 0.0) ? 1 : -1; x=fabs(x); c=0.797884560802865/sqrt(x); cosx=cos(x-0.706858347057703e1); sinx=sin(x-0.706858347057703e1); besspq1(x,&p1,&q1); return sgnx*c*(p1*sinx+q1*cosx); } }