~ --------------------------------------------------------------------- Bessel functions: This library defines the 0-th and 1st-order Bessel functions I_0(x), I_1(x), J_0(x), J_1(x), K_0(x), K_1(x). These are called as BesJ0(x), BesJ1(x), etc. and return the values of the evaluated function. Accuracy checked for 1<=x<=10 against Bronstein-Semendjajew's tables (4 digits). They agree within the displayed digits except BesI1, which seems somewhat inaccurate for x=4, 5, 6. Note: These are lifted straight out of "Numerical Recipes" (C version) and therefore assumed to be copyright by its authors and publisher. They are made available to the MathPad community as a service and for non-commercial use only. U. Wienands, Stanford Linear Accelerator Center, uli@slac.stanford.edu Parameters: x: the argument. Note: This file has to be saved with the same Degree/Radian setting as the "include"ing file (at least for MathPad 2.5.2), else BesJ will return erroneous values! ~ -- I_0(x) BesI0(x) = (I0_y:=(x/3.75)^2, 1.0+I0_y*(3.5156229+I0_y*(3.0899424+I0_y*(1.2067492 + I0_y*(0.2659732+I0_y*(0.360768e-1+I0_y*0.45813e-2))))) ) when abs(x) < 3.75, (I0_y:=3.75/abs(x), (exp(abs(x))/sqrt(abs(x)))*(0.39894228+I0_y*(0.1328592e-1 + I0_y*(0.225319e-2+I0_y*(-0.157565e-2+I0_y*(0.916281e-2 + I0_y*(-0.2057706e-1+I0_y*(0.2635537e-1+I0_y*(-0.1647633e-1 + I0_y*0.392377e-2)))))))) ) otherwise -- I_1(x) BesI1(x) = I1_ax:=abs(x), (I1_y:=(x/3.75)^2, (I1_ax*(0.5+I1_y*(0.87890594+I1_y*(0.51498869+I1_y*(0.15084934 + I1_y*(0.2658733e-1+I1_y*(0.301532e-2+I1_y*0.32411e-3))))))) * x/I1_ax ) when I1_ax<3.75, (I1_y:=3.75/I1_ax, I1_ans:=0.2282967e-1+I1_y*(-0.2895312e-1+I1_y*(0.1787654e-1, I1_y*0.420059e-2)), I1_ans:=0.39894228+I1_y*(-0.3988024e-1+I1_y*(-0.362018e-2 + I1_y*(0.163801e-2+I1_y*(-0.1031555e-1+I1_y*I1_ans)))), I1_ans*(exp(I1_ax)/sqrt(I1_ax))*x/I1_ax ) otherwise -- J_0(x) BesJ0(x) = (J0_y:=x*x, (57568490574.0+J0_y*(-13362590354.0+J0_y*(651619640.7 + J0_y*(-11214424.18+J0_y*(77392.33017+J0_y*(-184.9052456)))))) / (57568490411.0+J0_y*(1029532985.0+J0_y*(9494680.718 + J0_y*(59272.64853+J0_y*(267.8532712+J0_y*1.0))))) ) when abs(x)<8, (J0_ax:=abs(x), J0_z:=8/J0_ax, J0_y:=J0_z*J0_z, J0_xx:=J0_ax-0.785398164, sqrt(0.636619772/J0_ax) * (cos(J0_xx*Radians) * (1.0+J0_y*(-0.1098628627e-2+J0_y*(0.2734510407e-4 + J0_y*(-0.2073370639e-5+J0_y*0.2093887211e-6)))) - J0_z*sin(J0_xx*Radians)* (-0.1562499995e-1+J0_y*(0.1430488765e-3 + J0_y*(-0.6911147651e-5+J0_y*(0.7621095161e-6 - J0_y*0.934935152e-7))))) ) otherwise -- J_1(x) BesJ1(x) = (J1_y:=x*x, (x*(72362614232.0+J1_y*(-7895059235.0+J1_y*(242396853.1 + J1_y*(-2972611.439+J1_y*(15704.48260+J1_y*(-30.16036606))))))) / (144725228442.0+J1_y*(2300535178.0+J1_y*(18583304.74 + J1_y*(99447.43394+J1_y*(376.9991397+J1_y*1.0))))) ) when abs(x)<8, (J1_ax:=abs(x), J1_z:=8/J1_ax, J1_y:=J1_z*J1_z, J1_xx:=J1_ax-2.356194491, sqrt(0.636619772/J1_ax)*x/J1_ax * (cos(J1_xx*Radians) * (1.0+J1_y*(0.183105e-2+J1_y*(-0.3516396496e-4 + J1_y*(0.2457520174e-5+J1_y*(-0.240337019e-6))))) - J1_z*sin(J1_xx*Radians) * (0.04687499995+J1_y*(-0.2002690873e-3 + J1_y*(0.8449199096e-5+J1_y*(-0.88228987e-6 + J1_y*0.105787412e-6))))) ) otherwise -- K_0(x) BesK0(x) = (K0_y:=x*x/4, (-ln(x/2.0)*BesI0(x))+(-0.57721566+K0_y*(0.42278420 + K0_y*(0.23069756+K0_y*(0.3488590e-1+K0_y*(0.262698e-2 + K0_y*(0.10750e-3+K0_y*0.74e-5)))))) ) when x<=2, (K0_y:=2/x, (exp(-x)/sqrt(x))*(1.25331414+K0_y*(-0.7832358e-1 + K0_y*(0.2189568e-1+K0_y*(-0.1062446e-1+K0_y*(0.587872e-2 + K0_y*(-0.251540e-2+K0_y*0.53208e-3)))))) ) otherwise -- K_1(x) BesK1(x) = (K1_y:=x*x/4, (ln(x/2.0)*BesI1(x))+(1.0/x)*(1.0+K1_y*(0.15443144 + K1_y*(-0.67278579+K1_y*(-0.18156897+K1_y*(-0.1919402e-1 + K1_y*(-0.110404e-2+K1_y*(-0.4686e-4))))))) ) when x<=2, (K1_y:=2/x, (exp(-x)/sqrt(x))*(1.25331414+K1_y*(0.23498619 + K1_y*(-0.3655620e-1+K1_y*(0.1504268e-1+K1_y*(-0.780353e-2 + K1_y*(0.325614e-2+K1_y*(-0.68245e-3))))))) ) otherwise ~----------------------------------------------------------