!********************************************************* !* Program to demonstrate Bessel Coefficients Subroutine * !* ----------------------------------------------------- * !* 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: * !* * !* BESSEL COEFFICIENTS * !* What is the order of the Bessel function ? 0 * !* What degree is desired ? 41 * !* The coefficients are: *************** !* A( 0) = 0.100000E+01 A( 1) = 0.000000E+00 A( 2) = -0.250000E+00 * !* A( 3) = 0.000000E+00 A( 4) = 0.156250E-01 A( 5) = 0.000000E+00 * !* A( 6) = -0.434028E-03 A( 7) = 0.000000E+00 A( 8) = 0.678168E-05 * !* A( 9) = 0.000000E+00 A(10) = -0.678168E-07 A(11) = 0.000000E+00 * !* A(12) = 0.470950E-09 A(13) = 0.000000E+00 A(14) = -0.240281E-11 * !* A(15) = 0.000000E+00 A(16) = 0.938597E-14 A(17) = 0.000000E+00 * !* A(18) = -0.289690E-16 A(19) = 0.000000E+00 A(20) = 0.724226E-19 * !* A(21) = 0.000000E+00 A(22) = -0.149633E-21 A(23) = 0.000000E+00 * !* A(24) = 0.259780E-24 A(25) = 0.000000E+00 A(26) = -0.384290E-27 * !* A(27) = 0.000000E+00 A(28) = 0.490166E-30 A(29) = 0.000000E+00 * !* A(30) = -0.544629E-33 A(31) = 0.000000E+00 A(32) = 0.531864E-36 * !* A(33) = 0.000000E+00 A(34) = -0.460090E-39 A(35) = 0.000000E+00 * !* A(36) = 0.355008E-42 A(37) = 0.000000E+00 A(38) = -0.245850E-45 * !* A(39) = 0.000000E+00 A(40) = 0.153657E-48 A(41) = 0.000000E+00 * !* Argument ? 1 * !* Y= 0.765198 * !*********************************************************************** PROGRAM Bessel implicit none integer, parameter :: SIZE=100 real*8 A(0:SIZE) real*8 x,y integer i, m, n print *, 'BESSEL COEFFICIENTS' write(*,'(" What is the order of the Bessel function ? ")',advance='no') read *, n write(*,'(" What degree is desired ? ")',advance='no') read *, m print *, ' ' call Bessel_coeff(m,n,A) print *, 'The coefficients are:' write(*,50) ((i,A(i)),i=0,m) write(*,'(" Argument ? ")',advance='no') read *, x y = A(0) + A(1) * x; if (m > 1) y = y + A(2) * x * x; if (m > 2) y = y + A(3) * x * x * x; if (m > 3) y = y + A(4) * x**4 do i = 4, m if (m > i) y = y + A(i + 1) * x**(i+1) end do write(*,60) y stop 50 format(3(' A(',i2,') = ',e13.6)) 60 format(' Y=',f10.6/) end !************************************************************* !* Bessel function series coefficient evaluation subroutine * !* m+1 is the number of coefficients desired, n is the order * !* of the Bessel function. The coefficients are returned in * !* A(i). * !************************************************************* subroutine Bessel_coeff(m,n,A) implicit none integer i,m,n real*8 A(0:m+n), B(0:m+n) real*8 a1,b1; a1 = 1.d0; b1 = 1.d0; do i = 1, n B(i - 1) = 0.d0; b1 = b1 * i; a1 = a1 / 2.d0 end do b1 = a1 / b1; a1 = 1.d0 i = 0 do while (i <= m) A(i) = a1 * b1; A(i + 1) = 0 a1 = -a1 / ((i + 2) * (n + n + i + 2)) i = i + 2 end do a1 = a1 / 2.d0; do i = 0, m B(i + n) = A(i) end do A = B return end ! End file bessel.f90