[通用算法] 求助-复变量的贝塞尔第一类和第二类函数的Fortran子程序

SUBROUTINE cbdb(cz,cnu,fn,w)
!-----------------------------------------------------------------------

!         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)
DO  k = 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)
DO  i = 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) - 1  FOR -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) FOR  X >= 10  WHERE
!     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