!-----------------------------------------------------------------------
! CALCULATION OF J (CZ) BY THE DEBYE APPROXIMATION
! CNU
! ------------------
! IT IS ASSUMED THAT REAL(CZ) .GE. 0 AND THAT REAL(CNU) = FN + K
! WHERE K IS AN INTEGER.
!-----------------------------------------------------------------------
COMPLEX (dp), INTENT(IN) :: cz, cnu
REAL (dp), INTENT(IN) :: fn
COMPLEX (dp), INTENT(OUT):: w
! Local variables
REAL (dp) :: is, inu, izn
COMPLEX (dp):: c1, c2, eta, nu, p, p1, q, r, s, s1, s2, sm, t, z, zn
!----------------------
! C = 1/SQRT(2*PI)
! BND = PI/3
!----------------------
REAL (dp), PARAMETER:: c = .398942280401433_dp, pi = 3.14159265358979_dp,&
pi2 = 6.28318530717959_dp, bnd = 1.04719755119660_dp
COMPLEX (dp), PARAMETER :: j = (0.0, 1.0)
REAL (dp):: alpha, am, aq, ar
REAL (dp):: phi, sgn, theta
REAL (dp):: u, v, x, y
INTEGER :: ind, k, l, m
!----------------------
! COEFFICIENTS OF THE FIRST 16 POLYNOMIALS
! IN THE DEBYE APPROXIMATION
!----------------------
REAL (dp):: a(136) = (/ 1.0_dp, -.208333333333333_dp, .125000000000000_dp, .334201388888889_dp, &
-.401041666666667_dp, .703125000000000D-01,-.102581259645062D+01, .184646267361111D+01, &
-.891210937500000_dp, .732421875000000D-01, .466958442342625D+01,-.112070026162230D+02, &
.878912353515625D+01,-.236408691406250D+01, .112152099609375_dp,-.282120725582002D+02, &
.846362176746007D+02,-.918182415432400D+02, .425349987453885D+02,-.736879435947963D+01, &
.227108001708984_dp, .212570130039217D+03,-.765252468141182D+03, .105999045252800D+04, &
-.699579627376133D+03, .218190511744212D+03,-.264914304869516D+02, .572501420974731_dp, &
-.191945766231841D+04, .806172218173731D+04,-.135865500064341D+05, .116553933368645D+05, &
-.530564697861340D+04, .120090291321635D+04,-.108090919788395D+03, .172772750258446D+01, &
.202042913309661D+05,-.969805983886375D+05, .192547001232532D+06,-.203400177280416D+06, &
.122200464983017D+06,-.411926549688976D+05, .710951430248936D+04,-.493915304773088D+03, &
.607404200127348D+01,-.242919187900551D+06, .131176361466298D+07,-.299801591853811D+07, &
.376327129765640D+07,-.281356322658653D+07, .126836527332162D+07,-.331645172484564D+06, &
.452187689813627D+05,-.249983048181121D+04, .243805296995561D+02, .328446985307204D+07, &
-.197068191184322D+08, .509526024926646D+08,-.741051482115327D+08, .663445122747290D+08, &
-.375671766607634D+08, .132887671664218D+08,-.278561812808645D+07, .308186404612662D+06, &
-.138860897537170D+05, .110017140269247D+03,-.493292536645100D+08, .325573074185766D+09, &
-.939462359681578D+09, .155359689957058D+10,-.162108055210834D+10, .110684281682301D+10, &
-.495889784275030D+09, .142062907797533D+09,-.244740627257387D+08, .224376817792245D+07, &
-.840054336030241D+05, .551335896122021D+03, .814789096118312D+09,-.586648149205185D+10, &
.186882075092958D+11,-.346320433881588D+11, .412801855797540D+11,-.330265997498007D+11, &
.179542137311556D+11,-.656329379261928D+10, .155927986487926D+10,-.225105661889415D+09, &
.173951075539782D+08,-.549842327572289D+06, .303809051092238D+04,-.146792612476956D+11, &
.114498237732026D+12,-.399096175224466D+12, .819218669548577D+12,-.109837515608122D+13, &
.100815810686538D+13,-.645364869245377D+12, .287900649906151D+12,-.878670721780233D+11, &
.176347306068350D+11,-.216716498322380D+10, .143157876718889D+09,-.387183344257261D+07, &
.182577554742932D+05, .286464035717679D+12,-.240629790002850D+13, .910934118523990D+13, &
-.205168994109344D+14, .305651255199353D+14,-.316670885847852D+14, .233483640445818D+14, &
-.123204913055983D+14, .461272578084913D+13,-.119655288019618D+13, .205914503232410D+12, &
-.218229277575292D+11, .124700929351271D+10,-.291883881222208D+08, .118838426256783D+06, &
-.601972341723401D+13, .541775107551060D+14,-.221349638702525D+15, .542739664987660D+15, &
-.889496939881026D+15, .102695519608276D+16,-.857461032982895D+15, .523054882578445D+15, &
-.232604831188940D+15, .743731229086791D+14,-.166348247248925D+14, .248500092803409D+13, &
-.229619372968246D+12, .114657548994482D+11,-.234557963522252D+09, .832859304016289D+06 /)
z = cz
nu = cnu
inu = AIMAG(cnu)
IF (inu < 0.0D0) THEN
z = CONJG(z)
nu = CONJG(nu)
END IF
x = REAL(z, KIND=dp)
y = AIMAG(z)
! TANH(GAMMA) = SQRT(1 - (Z/NU)**2) = W/NU
! T = EXP(NU*(TANH(GAMMA) - GAMMA))
zn = z / nu
izn = AIMAG(zn)
IF (ABS(izn) <= 0.1D0*ABS(REAL(zn, KIND=dp))) THEN
s = (1.0D0-zn) * (1.0D0+zn)
eta = 1.0D0 / s
q = SQRT(s)
s = 1.0D0 / (nu*q)
t = zn / (1.0D0 + q)
t = EXP(nu*(q + LOG(t)))
ELSE
s = (nu-z) * (nu+z)
eta = (nu*nu) / s
w = SQRT(s)
q = w / nu
IF (REAL(q, KIND=dp) < 0.0D0) w = -w
s = 1.0D0 / w
t = z / (nu+w)
t = EXP(w + nu*LOG(t))
END IF
is = AIMAG(s)
r = SQRT(s)
c1 = r * t
ar = REAL(r, KIND=dp) * REAL(r, KIND=dp) + AIMAG(r) * AIMAG(r)
aq = -1.0D0 / (REAL(q, KIND=dp)*REAL(q, KIND=dp) + AIMAG(q)*AIMAG(q))
phi = ATAN2(y,x) / 3.0D0
q = nu - z
theta = ATAN2(AIMAG(q),REAL(q, KIND=dp)) - phi
ind = 0
IF (ABS(theta) >= 2.0D0*bnd) THEN
ind = 1
CALL dcrec(REAL(t, KIND=dp), AIMAG(t),u,v)
c2 = -j * r * CMPLX(u,v, KIND=dp)
IF (is >= 0.0D0) THEN
IF (is > 0.0D0) GO TO 10
IF (REAL(s, KIND=dp) <= 0.0D0) GO TO 10
END IF
c2 = -c2
END IF
! SUMMATION OF THE SERIES S1 AND S2
10 sm = s * s
p = (a(2)*eta + a(3)) * s
p1 = ((a(4)*eta + a(5))*eta + a(6)) * sm
s1 = (1.0D0 + p) + p1
IF (ind /= 0) s2 = (1.0D0-p) + p1
sgn = 1.0D0
am = ar * ar
m = 4
l = 6
! P = VALUE OF THE M-TH POLYNOMIAL
20 l = l + 1
alpha = a(l)
p = CMPLX(a(l),0.0D0, KIND=dp)
DOk = 2, m
l = l + 1
alpha = a(l) + aq * alpha
p = a(l) + eta * p
END DO
! ONLY THE S1 SUM IS FORMED WHEN IND = 0
sm = s * sm
p = p * sm
s1 = s1 + p
IF (ind /= 0) THEN
sgn = -sgn
s2 = s2 + sgn * p
END IF
am = ar * am
IF (1.0D0 + alpha*am /= 1.0D0) THEN
m = m + 1
IF (m <= 16) GO TO 20
END IF
! FINAL ASSEMBLY
s1 = c * c1 * s1
IF (ind == 0) THEN
w = s1
ELSE
s2 = c * c2 * s2
q = nu + z
theta = ATAN2(AIMAG(q),REAL(q, KIND=dp)) - phi
IF (ABS(theta) <= bnd) THEN
w = s1 + s2
ELSE
alpha = pi2
IF (izn < 0.0D0) alpha = -alpha
t = alpha * CMPLX(ABS(inu),-fn, KIND=dp)
alpha = EXP(REAL(t))
u = AIMAG(t)
r = CMPLX(COS(u),SIN(u), KIND=dp)
t = s1 - (alpha*r) * s1
IF (x == 0.0D0 .AND. inu == 0.0D0) t = -t
IF (y < 0.0D0) THEN
IF (izn >= 0.0D0 .AND. theta <= SIGN(pi,theta)) s2 =&
s2 * ( CONJG(r)/alpha)
IF (x == 0.0D0) GO TO 40
IF (izn >= 0.0D0) THEN
IF (is < 0.0D0) GO TO 40
END IF
END IF
w = s2 + t
GO TO 50
40 w = s2 - t
END IF
END IF
50 IF (inu < 0.0D0) w = CONJG(w)
RETURN
END SUBROUTINE cbdb
SUBROUTINE cbja(cz,cnu,w)
!-----------------------------------------------------------------------
! COMPUTATION OF J(NU,Z) BY THE ASYMPTOTIC EXPANSION
!-----------------------------------------------------------------------
COMPLEX (dp), INTENT(IN) :: cz
COMPLEX (dp), INTENT(IN) :: cnu
COMPLEX (dp), INTENT(OUT):: w
! Local variables
REAL (dp) :: eps, inu, m
COMPLEX (dp):: a, a1, arg, e, eta, nu, p, q, t, z, zr, zz
!--------------------------
REAL (dp) :: r, rnu, tol, u, v
REAL (dp) :: x, y
INTEGER :: i, ind
!--------------------------
! PIHALF = PI/2
! C = 2*PI**(-1/2)
!--------------------------
REAL (dp), PARAMETER :: pihalf = 1.5707963267949_dp, c = 1.12837916709551_dp
COMPLEX (dp), PARAMETER :: j = (0.0_dp, 1.0_dp)
!--------------------------
! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
! SMALLEST NUMBER SUCH THAT 1.0 + EPS .GT. 1.0 .
eps = EPSILON(0.0_dp)
!--------------------------
z = cz
x = REAL(z, KIND=dp)
y = AIMAG(z)
nu = cnu
ind = 0
IF (ABS(x) <= 1.d-2*ABS(y)) THEN
IF (AIMAG(nu) < 0.0D0 .AND. ABS(REAL(nu)) < 1.d-2*ABS(AIMAG(nu))) THEN
ind = 1
nu = CONJG(nu)
z = CONJG(z)
y = -y
END IF
END IF
IF (x < -1.d-2*y) z = -z
zz = z + z
CALL dcrec(REAL(zz, KIND=dp),AIMAG(zz),u,v)
zr = CMPLX(u,v, KIND=dp)
eta = -zr * zr
p = (0.0D0,0.0D0)
q = (0.0D0,0.0D0)
a1 = nu * nu - 0.25D0
a = a1
t = a1
m = 1.0D0
tol = eps * anorm(a1)
DOi = 1, 16
a = a - 2.0D0 * m
m = m + 1.0D0
t = t * a * eta / m
p = p + t
a = a - 2.0D0 * m
m = m + 1.0D0
t = t * a / m
q = q + t
IF (anorm(t) <= tol) GO TO 20
END DO
20 p = p + 1.0D0
q = (q+a1) * zr
w = z - pihalf * nu
IF (ABS(AIMAG(w)) <= 1.0D0) THEN
arg = w - 0.5D0 * pihalf
w = c * SQRT(zr) * (p*COS(arg) - q*SIN(arg))
ELSE
e = EXP(-j*w)
t = q - j * p
IF (AIMAG(z) > 0.0D0 .AND. REAL(z, KIND=dp) <= 1.d-2*AIMAG(z).AND.&
ABS(REAL(nu, KIND=dp)) < 1.d-2*AIMAG(nu)) t = 0.5D0 * t
CALL dcrec(REAL(e, KIND=dp),AIMAG(e),u,v)
w = 0.5D0 * c * SQRT(j*zr) * ((p-j*q)*e + t*CMPLX(u,v, KIND=dp))
END IF
IF (x < -1.d-2*y) THEN
IF (y < 0.0D0) nu = -nu
! COMPUTATION OF EXP(I*PI*NU)
rnu = REAL(nu, KIND=dp)
inu = AIMAG(nu)
r = EXP(-2.0D0*pihalf*inu)
u = r * dcos1(rnu)
v = r * dsin1(rnu)
w = w * CMPLX(u,v, KIND=dp)
END IF
IF (ind /= 0) w = CONJG(w)
RETURN
END SUBROUTINE cbja
FUNCTION anorm(z) RESULT(fn_val)
! Replaces the statement function anorm in the F77 code.
COMPLEX (dp), INTENT(IN):: z
REAL (dp) :: fn_val
fn_val = MAX( ABS( REAL(z, KIND=dp)), ABS(AIMAG(z) ) )
RETURN
END FUNCTION anorm
FUNCTION dgam1(x) RESULT(fn_val)
!-----------------------------------------------------------------------
! EVALUATION OF 1/GAMMA(1 + X) - 1FOR -0.5 <= X <= 1.5
!-----------------------------------------------------------------------
! THE FOLLOWING ARE THE FIRST 49 COEFFICIENTS OF THE MACLAURIN
! EXPANSION FOR 1/GAMMA(1 + X) - 1. THE COEFFICIENTS ARE
! CORRECT TO 40 DIGITS. THE COEFFICIENTS WERE OBTAINED BY
! ALFRED H. MORRIS JR. (NAVAL SURFACE WARFARE CENTER) AND ARE
! GIVEN HERE FOR REFERENCE. ONLY THE FIRST 14 COEFFICIENTS ARE
! USED IN THIS CODE.
! -----------
! DATA A(1)/ .5772156649015328606065120900824024310422D+00/,
! * A(2)/-.6558780715202538810770195151453904812798D+00/,
! * A(3)/-.4200263503409523552900393487542981871139D-01/,
! * A(4)/ .1665386113822914895017007951021052357178D+00/,
! * A(5)/-.4219773455554433674820830128918739130165D-01/,
! * A(6)/-.9621971527876973562114921672348198975363D-02/,
! * A(7)/ .7218943246663099542395010340446572709905D-02/,
! * A(8)/-.1165167591859065112113971084018388666809D-02/,
! * A(9)/-.2152416741149509728157299630536478064782D-03/,
! * A(10) / .1280502823881161861531986263281643233949D-03/
! DATA A(11) /-.2013485478078823865568939142102181838229D-04/,
! * A(12) /-.1250493482142670657345359473833092242323D-05/,
! * A(13) / .1133027231981695882374129620330744943324D-05/,
! * A(14) /-.2056338416977607103450154130020572836513D-06/,
! * A(15) / .6116095104481415817862498682855342867276D-08/,
! * A(16) / .5002007644469222930055665048059991303045D-08/,
! * A(17) /-.1181274570487020144588126565436505577739D-08/,
! * A(18) / .1043426711691100510491540332312250191401D-09/,
! * A(19) / .7782263439905071254049937311360777226068D-11/,
! * A(20) /-.3696805618642205708187815878085766236571D-11/
! DATA A(21) / .5100370287454475979015481322863231802727D-12/,
! * A(22) /-.2058326053566506783222429544855237419746D-13/,
! * A(23) /-.5348122539423017982370017318727939948990D-14/,
! * A(24) / .1226778628238260790158893846622422428165D-14/,
! * A(25) /-.1181259301697458769513764586842297831212D-15/,
! * A(26) / .1186692254751600332579777242928674071088D-17/,
! * A(27) / .1412380655318031781555803947566709037086D-17/,
! * A(28) /-.2298745684435370206592478580633699260285D-18/,
! * A(29) / .1714406321927337433383963370267257066813D-19/,
! * A(30) / .1337351730493693114864781395122268022875D-21/
! DATA A(31) /-.2054233551766672789325025351355733796682D-21/,
! * A(32) / .2736030048607999844831509904330982014865D-22/,
! * A(33) /-.1732356445910516639057428451564779799070D-23/,
! * A(34) /-.2360619024499287287343450735427531007926D-25/,
! * A(35) / .1864982941717294430718413161878666898946D-25/,
! * A(36) /-.2218095624207197204399716913626860379732D-26/,
! * A(37) / .1297781974947993668824414486330594165619D-27/,
! * A(38) / .1180697474966528406222745415509971518560D-29/,
! * A(39) /-.1124584349277088090293654674261439512119D-29/,
! * A(40) / .1277085175140866203990206677751124647749D-30/
! DATA A(41) /-.7391451169615140823461289330108552823711D-32/,
! * A(42) / .1134750257554215760954165259469306393009D-34/,
! * A(43) / .4639134641058722029944804907952228463058D-34/,
! * A(44) /-.5347336818439198875077418196709893320905D-35/,
! * A(45) / .3207995923613352622861237279082794391090D-36/,
! * A(46) /-.4445829736550756882101590352124643637401D-38/,
! * A(47) /-.1311174518881988712901058494389922190237D-38/,
! * A(48) / .1647033352543813886818259327906394145400D-39/,
! * A(49) /-.1056233178503581218600561071538285049997D-40/
! -----------
! C = A(1) - 1 IS ALSO FREQUENTLY NEEDED. C HAS THE VALUE ...
! DATA C /-.4227843350984671393934879099175975689578D+00/
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN) :: x
REAL (dp) :: fn_val
! Local variables
REAL (dp) :: d, t, w, z
REAL (dp), PARAMETER :: a0 = .611609510448141581788D-08, a1&
= .624730830116465516210D-08, b1 = .203610414066806987300D+00, b2&
= .266205348428949217746D-01, b3 = .493944979382446875238D-03, b4&
= -.851419432440314906588D-05, b5 = -.643045481779353022248D-05, b6&
= .992641840672773722196D-06, b7 = -.607761895722825260739D-07, b8&
= .195755836614639731882D-09
REAL (dp), PARAMETER :: p0 = .6116095104481415817861D-08, p1&
= .6871674113067198736152D-08, p2 = .6820161668496170657, p3&
= .4686843322948848031080D-10, p4 = .1572833027710446286995D-11, p5&
= -.1249441572276366213222D-12, p6 = .4343529937408594255178D-14, q1&
= .3056961078365221025009D+00, q2 = .5464213086042296536016D-01, q3&
= .4956830093825887312, q4 = .2692369466186361192876D-03
REAL (dp), PARAMETER :: c = -.422784335098467139393487909917598D+00, c0&
= .577215664901532860606512090082402D+00, c1&
= -.655878071520253881077019515145390D+00, c2&
= -.420026350340952355290039348754298D-01, c3&
= .166538611382291489501700795102105D+00, c4&
= -.421977345555443367482083012891874D-01, c5&
= -.962197152787697356211492167234820D-02, c6&
= .721894324666309954239501034044657D-02, c7&
= -.116516759185906511211397108401839D-02, c8&
= -.215241674114950972815729963053648D-03, c9&
= .128050282388116186153198626328164D-03, c10&
= -.201348547807882386556893914210218D-04, c11&
= -.125049348214267065734535947383309D-05, c12&
= .113302723198169588237412962033074D-05, c13&
= -.205633841697760710345015413002057D-06
!----------------------------
t = x
d = x - 0.5_dp
IF (d > 0._dp) t = d - 0.5_dp
IF (t < 0.0_dp) THEN
GO TO 30
ELSE IF (t > 0.0_dp) THEN
GO TO 20
END IF
fn_val = 0._dp
RETURN
!------------
! CASE WHEN 0 < T <= 0.5
! W IS A MINIMAX APPROXIMATION FOR
! THE SERIES A(15) + A(16)*T + ...
!------------
20 w = ((((((p6*t + p5)*t + p4)*t + p3)*t + p2)*t + p1)*t + p0) / &
((((q4*t+q3)*t + q2)*t + q1)*t + 1._dp)
z = (((((((((((((w*t + c13)*t + c12)*t + c11)*t + c10)*t + c9)*t + c8)*t + c7)*t&
+ c6)*t + c5)*t + c4)*t + c3)*t + c2)*t + c1) * t + c0
IF (d <= 0._dp) THEN
fn_val = x * z
RETURN
END IF
fn_val = (t/x) * ((z-0.5_dp)-0.5_dp)
RETURN
!------------
! CASE WHEN -0.5 <= T < 0
! W IS A MINIMAX APPROXIMATION FOR
! THE SERIES A(15) + A(16)*T + ...
!------------
30 w = (a1*t + a0) / ((((((((b8*t + b7)*t + b6)*t + b5)*t + b4)*t + b3)*t + b2)*t + b1)*t+1._dp)
z = (((((((((((((w*t + c13)*t + c12)*t + c11)*t + c10)*t + c9)*t + c8)*t + c7)*t&
+ c6)*t + c5)*t + c4)*t + c3)*t + c2)*t + c1) * t + c
IF (d <= 0._dp) THEN
fn_val = x * ((z+0.5_dp)+0.5_dp)
RETURN
END IF
fn_val = t * z / x
RETURN
END FUNCTION dgam1
FUNCTION dpdel(x) RESULT(fn_val)
!-----------------------------------------------------------------------
! COMPUTATION OF THE FUNCTION DEL(X) FORX >= 10WHERE
! LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X)
! --------
! THE SERIES FOR DPDEL ON THE INTERVAL 0.0 TO 1.0 DERIVED BY
! A.H. MORRIS FROM THE CHEBYSHEV SERIES IN THE SLATEC LIBRARY
! OBTAINED BY WAYNE FULLERTON (LOS ALAMOS).
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN) :: x
REAL (dp) :: fn_val
! Local variables
REAL (dp), PARAMETER :: a(15) = (/ .833333333333333333333333333333D-01,&
-.277777777777777777777777752282D-04,&
.793650793650793650791732130419D-07,&
-.595238095238095232389839236182D-09,&
.841750841750832853294451671990D-11,&
-.191752691751854612334149171243D-12,&
.641025640510325475730918472625D-14,&
-.295506514125338232839867823991D-15,&
.179643716359402238723287696452D-16,&
-.139228964661627791231203060395D-17,&
.133802855014020915603275339093D-18,&
-.154246009867966094273710216533D-19,&
.197701992980957427278370133333D-20,&
-.234065664793997056856992426667D-21,&
.171348014966398575409015466667D-22 /)
REAL (dp) :: t, w
INTEGER :: i, k
!-----------------------------------------------------------------------
t = (10._dp/x) ** 2
w = a(15)
DO i = 1, 14
k = 15 - i
w = t * w + a(k)
END DO
fn_val = w / x
RETURN
END FUNCTION dpdel
END MODULE Complex_Bessel
风平老涡 发表于 2021-2-9 00:05
SUBROUTINE cbsslj(z,cnu,w)
!---------------------------------------------- ...
好久没登陆了，非常感谢您的帮助。但是我运行了一下您提供的这个程序，编译的时候还是出现很多比如变量未声明这样的问题。 jlg1206 发表于 2021-4-8 19:16
好久没登陆了，非常感谢您的帮助。但是我运行了一下您提供的这个程序，编译的时候还是出现很多比如变量未 ...
附件应该有俩个文件。 搜netlib里面的amos，里面有
页:
1
[2]