| Re: calculation of bessel functions j_l(x) |
|
 |
|
 |
|
 |
|
 |
Group: comp.lang.fortran · Group Profile
Author: FatemehFatemeh Date: Sep 13, 2008 04:16
Dear all ;
There is a command at the first line of bessj.f90 that when I comment
it, the program runs without any error. (this command is : call
assert(n >= 2, 'bessj_s args') )
bessj.f90:
--------------------
FUNCTION bessj_s(n,x)
USE nrtype; USE nrutil, ONLY : assert
USE nr, ONLY : bessj0,bessj1
IMPLICIT NONE
INTEGER(I4B), INTENT(IN) :: n
REAL(SP), INTENT(IN) :: x
REAL(SP) :: bessj_s
INTEGER(I4B), PARAMETER :: IACC=40,IEXP=maxexponent(x)/2
INTEGER(I4B) :: j,jsum,m
REAL(SP) :: ax,bj,bjm,bjp,summ,tox
call assert(n >= 2, 'bessj_s args')
ax=abs(x)
if (ax*ax <= 8.0_sp*tiny(x)) then
bessj_s=0.0
else if (ax > real(n,sp)) then
tox=2.0_sp/ax
bjm=bessj0(ax)
bj=bessj1(ax)
do j=1,n-1
bjp=j*tox*bj-bjm
bjm=bj
bj=bjp
end do
bessj_s=bj
else
tox=2.0_sp/ax
m=2*((n+int(sqrt(real(IACC*n,sp))))/2)
bessj_s=0.0
jsum=0
summ=0.0
bjp=0.0
bj=1.0
do j=m,1,-1
bjm=j*tox*bj-bjp
bjp=bj
bj=bjm
if (exponent(bj) > IEXP) then
bj=scale(bj,-IEXP)
bjp=scale(bjp,-IEXP)
.
.
.
----------------------------------------------
INTERFACE assert
MODULE PROCEDURE
assert1,assert2,assert3,assert4,assert_v
END INTERFACE
---------------------------------
SUBROUTINE assert1(n1,string)
CHARACTER(LEN=*), INTENT(IN) :: string
LOGICAL, INTENT(IN) :: n1
if (.not. n1) then
write (*,*) 'nrerror: an assertion failed with this
tag:', &
string
STOP 'program terminated by assert1'
end if
END SUBROUTINE assert1
-------------------------------------------
Would you please tell me that is it true to comment this line:
call assert(n >= 2, 'bessj_s args')
Best regards,
Fatemeh
|