!****************************************************************** !* Program to demonstrate the Bessel function asymptotic series * !* -------------------------------------------------------------- * !* Reference: BASIC Scientific Subroutines, Vol. II * !* by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981. * !* * !* F90 Version by J.-P. Moreau, Paris. * !* (www.jpmoreau.fr) * !* -------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* X J0(X) J1(X) N E * !* -------------------------------------------------------------- * !* 1 0.733562 0.402234 1 0.1121521 * !* 2 0.221488 0.578634 1 0.0070095 * !* 3 -0.259956 0.340699 2 0.0007853 * !* 4 -0.396826 -0.065886 2 0.0001398 * !* 5 -0.177486 -0.327696 3 0.0000155 * !* 6 0.150635 -0.276760 3 0.0000036 * !* 7 0.300051 -0.004696 4 0.0000004 * !* 8 0.171638 0.234650 5 0.0000000 * !* 9 -0.090332 0.245324 5 0.0000000 * !* 10 -0.245930 0.043476 6 0.0000000 * !* 11 -0.171187 -0.176788 5 0.0000000 * !* 12 0.047689 -0.223450 5 0.0000000 * !* 13 0.206925 -0.070319 4 0.0000000 * !* 14 0.171072 0.133376 4 0.0000000 * !* 15 -0.014225 0.205105 4 0.0000000 * !* -------------------------------------------------------------- * !* * !****************************************************************** PROGRAM Bessel1 implicit none real*8 e, e3, J0, J1, x integer i, n e3=1.d-8 print *, ' ' print *,' X J0(X) J1(X) N E ' print *,'------------------------------------------------------------------' do i = 1, 15 x = dfloat(i) call Calcul_Bessel(n,e,e3,x,J0,J1) write(*,100) i,J0,J1,n,e end do print *,'------------------------------------------------------------------' stop 100 format(' ',i2,' ',f9.6,' ',f9.6,' ',i1,' ',f9.7) end !********************************************************** !* Bessel function asymptotic series subroutine. This * !* program calculates the zeroth and first order Bessel * !* functions using an asymptotic series expansion. * !* The required input are X and a convergence factor e. * !* Returned are the two Bessel functions, J0(X) and J1(X) * !* ------------------------------------------------------ * !* Reference; Algorithms for RPN calculators, by Ball, * !* L.A. Wiley and sons. * !********************************************************** Subroutine Calcul_Bessel(n,e,e3,x,J0,J1) ! Labels: 100, 200 implicit none real*8 a,a1,a2,b,c,e,e1,e2,e3,e4,J0,J1,pi,x,x1 real*8 m1, m2, n1, n2 integer m, n pi = 4.d0*atan(1.d0) !Calculate P and Q polynomials; P0(x)=m1, P1(x)=m2, !Q0(x)=n1, Q1(x)=n2 a = 1.d0; a1 = 1.d0; a2 = 1.d0; b = 1.d0; c = 1.d0 e1 = 1.d6 m = -1 x1 = 1.d0 / (8.d0 * x) x1 = x1 * x1 m1 = 1.d0; m2 = 1.d0 n1 = -1.d0 / (8.d0 * x) n2 = -3.d0 * n1 n = 0 100 m = m + 2; a = a * m * m m = m + 2; a = a * m * m c = c * x1 a1 = a1 * a2; a2 = a2 + 1.d0; a1 = a1 * a2 a2 = a2 + 1.d0; e2 = a * c / a1 e4 = 1.d0 + (m + 2) / m + (m + 2) * (m + 2) / (a2 * 8.d0 * x) + (m + 2) * (m + 4) / (a2 * 8.d0 * x); e4 = e4 * e2 !Test for divergence if (dabs(e4) > e1) goto 200 e1 = dabs(e2) m1 = m1 - e2 m2 = m2 + e2 * (m + 2) / m n1 = n1 + e2 * (m + 2) * (m + 2) / (a2 * 8.d0 * x) n2 = n2 - e2 * (m + 2) * (m + 4) / (a2 * 8.d0 * x) n = n + 1 !Test for convergence criterion if (e1 < e3) goto 200 goto 100 200 a = pi e = e2 b = dsqrt(2.d0 / (a * x)) J0 = b * (m1 * dcos(x - a / 4.d0) - n1 * dsin(x - a / 4.d0)) J1 = b * (m2 * dcos(x - 3.d0 * a / 4.d0) - n2 * dsin(x - 3.d0 * a / 4.d0)) return end ! End of file bessel1.f90