!****************************************************************** !* Complex Bessel Function of the 1st Kind of integer order * !* -------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Complex Bessel Function of the 1st Kind of integer order * !* * !* Input complex argument (real imaginary): 1 2 * !* Input integer order: 1 * !* * !* Function value: (1.291847,1.010488) * !* * !* * !* F90 Release 1.1 By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !* * !* Release 1.1: MAXK=20 (15 before). * !****************************************************************** Program Test_CBESSJ integer nu complex z,z1 real*8 x,y print *,' ' print *,' Complex Bessel Function of the 1st Kind of integer order' print *,' ' write(*,10,advance='no'); read *, x, y z = CMPLX(x,y) write(*,20,advance='no'); read *, nu print *,' ' call CBESSJ(z,nu,z1) print *,' Function value: ', z1 print *,' ' 10 format(' Input complex argument (real imaginary): ') 20 format(' Input integer order: ') END real*8 Function Fact(K) Integer i Real*8 f F=1.d0 do i=2, k f=f*dfloat(i) end do Fact=f return End !******************************************* !* FUNCTION GAMMA(X) * !* --------------------------------------- * !* Returns the value of Gamma(x) in double * !* precision as EXP(LN(GAMMA(X))) for X>0. * !******************************************* real*8 Function Gamma(xx) parameter(ONE=1.d0,FPF=5.5d0,HALF=0.5d0) real*8 xx real*8 cof(6) real*8 stp,x,tmp,ser integer j cof(1)=76.18009173d0 cof(2)=-86.50532033d0 cof(3)=24.01409822d0 cof(4)=-1.231739516d0 cof(5)=0.120858003d-2 cof(6)=-0.536382d-5 stp=2.50662827465d0 x=xx-ONE tmp=x+FPF tmp=(x+HALF)*LOG(tmp)-tmp ser=ONE do j=1, 6 x=x+ONE ser=ser+cof(j)/x end do Gamma = EXP(tmp+LOG(stp*ser)) return End Subroutine CBESSJ(z, nu, z1) !--------------------------------------------------- ! inf. (-z^2/4)^k ! Jnu(z) = (z/2)^nu x Sum ------------------ ! k=0 k! x Gamma(nu+k+1) ! (nu must be >= 0). !--------------------------------------------------- Parameter(MAXK=20,ZERO=0.d0) Complex z,z1 Integer k Complex sum,tmp Real*8 Fact, Gamma sum = CMPLX(ZERO,ZERO) do k=0, MAXK !calculate (-z**2/4)**k tmp = (-z*z/4.d0)**k !divide by k! tmp = tmp / Fact(k) !divide by Gamma(nu+k+1) tmp = tmp / Gamma(dfloat(nu+k+1)) !actualize sum sum = sum + tmp end do !calculate (z/2)**nu tmp = (z/2)**nu !multiply (z/2)**nu by sum z1 = tmp*sum return End !end of file cbessj.f90