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

MODULE Complex_Bessel
! SUBROUTINE CBSSLJ evaluates Bessel functions of complex ORDER
! SUBROUTINE CBSSLJ(x, nu, w)
! x = argument, nu = order, w = function value

! Latest revision - 20 July 2001
! amiller @ bigpond.net.au

! WARNING (from the NSWC Math. Library Manual)
! Precision:
! The real and imaginary parts are normally accurate to 11-12 significant
! digits when the function value is not near zero.   The exception is when
! Real(z) is small, Imag(nu) is small, and 14 < |z| < 17.5 + 0.5|nu|^2
! In this region, the real part of the result is still accurate, but all
! accuracy may be lost in the imaginary part.

USE constants_NSWC
IMPLICIT NONE

CONTAINS


FUNCTION cpabs(x, y) RESULT(fn_val)
!     --------------------------------------
!     EVALUATION OF SQRT(X*X + Y*Y)
!     --------------------------------------
REAL (dp), INTENT(IN)  :: x, y
REAL (dp)              :: fn_val

! Local variable
REAL (dp)  :: a

IF (ABS(x) > ABS(y)) THEN
  a = y / x
  fn_val = ABS(x) * SQRT(1.0_dp + a*a)
  RETURN
END IF
IF (y /= 0.0_dp) THEN
  a = x / y
  fn_val = ABS(y) * SQRT(1.0_dp + a*a)
  RETURN
END IF
fn_val = 0.0_dp
RETURN
END FUNCTION cpabs

SUBROUTINE dcrec(x, y, u, v)
!-----------------------------------------------------------------------
!             COMPLEX RECIPROCAL U + I*V = 1/(X + I*Y)
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN)   :: x, y
REAL (dp), INTENT(OUT)  :: u, v

! Local variables
REAL (dp)  :: t, d

IF (ABS(x) <= ABS(y)) THEN
  t = x / y
  d = y + t * x
  u = t / d
  v = -1.0_dp / d
  RETURN
END IF
t = y / x
d = x + t * y
u = 1.0_dp / d
v = -t / d
RETURN
END SUBROUTINE dcrec



FUNCTION cdiv(a,b) RESULT(fn_val)
!-----------------------------------------------------------------------
!              COMPLEX DIVISION A/B WHERE B IS NONZERO
!-----------------------------------------------------------------------

COMPLEX (dp), INTENT(IN)  :: a, b
COMPLEX (dp)              :: fn_val


REAL (dp)  :: ai, ar, bi, br, d, t
REAL (dp)  :: u, v

ar = REAL(a, KIND=dp)
ai = AIMAG(a)
br = REAL(b, KIND=dp)
bi = AIMAG(b)

IF (ABS(br) >= ABS(bi)) THEN
  t = bi / br
  d = br + t * bi
  u = (ar+ai*t) / d
  v = (ai-ar*t) / d
  fn_val = CMPLX(u,v, KIND=dp)
  RETURN
END IF
t = br / bi
d = bi + t * br
u = (ar*t+ai) / d
v = (ai*t-ar) / d
fn_val = CMPLX(u,v, KIND=dp)
RETURN
END FUNCTION cdiv

FUNCTION dsin1(x) RESULT(fn_val)
!-----------------------------------------------------------------------

!                REAL (dp) EVALUATION OF SIN(PI*X)

!                             --------------

!     THE EXPANSION FOR SIN(PI*A) (ABS(A) <= PI/4) USING A1,...,A13
!     IS ACCURATE TO WITHIN 2 UNITS OF THE 40-TH SIGNIFICANT DIGIT, AND
!     THE EXPANSION FOR COS(PI*A) (ABS(A) <= PI/4) USING B1,...,B13
!     IS ACCURATE TO WITHIN 4 UNITS OF THE 40-TH SIGNIFICANT DIGIT.

!-----------------------------------------------------------------------
REAL (dp), INTENT(IN)  :: x
REAL (dp)              :: fn_val

! Local variables
REAL (dp)            :: a, t, w
REAL (dp), PARAMETER :: pi = 3.141592653589793238462643383279502884197D+00
REAL (dp), PARAMETER :: a1 = -.1028083791780141522795259479153765743002D+00,   &
      a2  = .3170868848763100170457042079710451905600D-02,   &
      a3  = -.4657026956105571623449026167864697920000D-04,  &
      a4  = .3989844942879455643410226655783424000000D-06,   &
      a5  = -.2237397227721999776371894030796800000000D-08,  &
      a6  = .8847045483056962709715066675200000000000D-11,   &
      a7  = -.2598715447506450292885585920000000000000D-13,  &
      a8  = .5893449774331011070033920000000000000000D-16 ,  &
      a9  = -.1062975472045522550784000000000000000000D-18,   &
      a10 = .1561182648301780992000000000000000000000D-21,    &
      a11 = -.1903193516670976000000000000000000000000D-24,   &
      a12 = .1956617650176000000000000000000000000000D-27,    &
      a13 = -.1711276032000000000000000000000000000000D-30
REAL (dp), PARAMETER :: b1 = -.3084251375340424568385778437461297229882D+00, &
      b2  = .1585434424381550085228521039855226435920D-01,   &
      b3  = -.3259918869273900136414318317506279360000D-03,  &
      b4  = .3590860448591510079069203991239232000000D-05,   &
      b5  = -.2461136950494199754009084061808640000000D-07,  &
      b6  = .1150115912797405152263195572224000000000D-09,   &
      b7  = -.3898073171259675439899172864000000000000D-12,  &
      b8  = .1001886461636271969091584000000000000000D-14,   &
      b9  = -.2019653396886572027084800000000000000000D-17,  &
      b10 = .3278483561466560512000000000000000000000D-20,   &
      b11 = -.4377345082051788800000000000000000000000D-23,  &
      b12 = .4891532381388800000000000000000000000000D-26,   &
      b13 = -.4617089843200000000000000000000000000000D-29
INTEGER  :: max, n
!------------------------

!     ****** MAX IS A MACHINE DEPENDENT CONSTANT. MAX IS THE
!            LARGEST POSITIVE INTEGER THAT MAY BE USED.

!                       MAX = IPMPAR(3)
max = HUGE(3)

!------------------------
a = ABS(x)
t = MAX
IF (a >= t) THEN
  fn_val = 0.0_dp
  RETURN
END IF

n = a
t = n
a = a - t
IF (a <= 0.75_dp) THEN
  IF (a < 0.25_dp) GO TO 10

!                    0.25 <= A <= 0.75

  a = 0.25_dp + (0.25_dp-a)
  t = 16._dp * a * a
  fn_val = (((((((((((((b13*t + b12)*t + b11)*t + b10)*t + b9)*t + b8)*t  &
           + b7)*t + b6)*t + b5)*t + b4)*t + b3)*t + b2)*t + b1)*t +  &
           0.5_dp) + 0.5_dp
  GO TO 20
END IF

!                 A < 0.25  OR  A > 0.75

a = 0.25_dp + (0.75_dp-a)
10 t = 16._dp * a * a
w = (((((((((((((a13*t + a12)*t + a11)*t + a10)*t + a9)*t + a8)*t + a7)*t  &
    + a6)*t + a5)*t + a4)*t + a3)*t + a2)*t + a1)*t + 0.5_dp) + 0.5_dp
fn_val = pi * a * w

!                        TERMINATION

20 IF (x < 0.0) fn_val = -fn_val
IF (MOD(n,2) /= 0) fn_val = -fn_val
RETURN
END FUNCTION dsin1



FUNCTION dcos1 (x) RESULT(fn_val)
 
!-----------------------------------------------------------------------

!                REAL (dp) EVALUATION OF COS(PI*X)

!                             --------------

!     THE EXPANSION FOR SIN(PI*A) (ABS(A) .LE. PI/4) USING A1,...,A13
!     IS ACCURATE TO WITHIN 2 UNITS OF THE 40-TH SIGNIFICANT DIGIT, AND
!     THE EXPANSION FOR COS(PI*A) (ABS(A) .LE. PI/4) USING B1,...,B13
!     IS ACCURATE TO WITHIN 4 UNITS OF THE 40-TH SIGNIFICANT DIGIT.

!-----------------------------------------------------------------------

REAL (dp), INTENT(IN)  :: x
REAL (dp)              :: fn_val

REAL (dp)  :: a, t, w
INTEGER    :: MAX, n
!------------------------
REAL (dp), PARAMETER  :: pi = 3.141592653589793238462643383279502884197_dp
!------------------------
REAL (dp), PARAMETER  :: &
    a1  = -.1028083791780141522795259479153765743002D+00,  &
    a2  =  .3170868848763100170457042079710451905600D-02,  &
    a3  = -.4657026956105571623449026167864697920000D-04,  &
    a4  =  .3989844942879455643410226655783424000000D-06,  &
    a5  = -.2237397227721999776371894030796800000000D-08,  &
    a6  =  .8847045483056962709715066675200000000000D-11,  &
    a7  = -.2598715447506450292885585920000000000000D-13,  &
    a8  =  .5893449774331011070033920000000000000000D-16,  &
    a9  = -.1062975472045522550784000000000000000000D-18,  &
    a10 =  .1561182648301780992000000000000000000000D-21,  &
    a11 = -.1903193516670976000000000000000000000000D-24,  &
    a12 =  .1956617650176000000000000000000000000000D-27,  &
    a13 = -.1711276032000000000000000000000000000000D-30
!------------------------
REAL (dp), PARAMETER  :: &
    b1  = -.3084251375340424568385778437461297229882D+00,  &
    b2  =  .1585434424381550085228521039855226435920D-01,  &
    b3  = -.3259918869273900136414318317506279360000D-03,  &
    b4  =  .3590860448591510079069203991239232000000D-05,  &
    b5  = -.2461136950494199754009084061808640000000D-07,  &
    b6  =  .1150115912797405152263195572224000000000D-09,  &
    b7  = -.3898073171259675439899172864000000000000D-12,  &
    b8  =  .1001886461636271969091584000000000000000D-14,  &
    b9  = -.2019653396886572027084800000000000000000D-17,  &
    b10 =  .3278483561466560512000000000000000000000D-20,  &
    b11 = -.4377345082051788800000000000000000000000D-23,  &
    b12 =  .4891532381388800000000000000000000000000D-26,  &
    b13 = -.4617089843200000000000000000000000000000D-29
!------------------------

!     ****** MAX IS A MACHINE DEPENDENT CONSTANT. MAX IS THE
!            LARGEST POSITIVE INTEGER THAT MAY BE USED.

MAX = HUGE(0)

!------------------------
a = ABS(x)
t = MAX
IF (a < t) GO TO 10
fn_val = 1.d0
RETURN

10 n = a
t = n
a = a - t
IF (a > 0.75D0) GO TO 20
IF (a < 0.25D0) GO TO 21

!                    0.25 .LE. A .LE. 0.75

a = 0.25D0 + (0.25D0 - a)
t = 16.d0*a*a
w = (((((((((((((a13*t + a12)*t + a11)*t + a10)*t + a9)*t +  &
    a8)*t + a7)*t + a6)*t + a5)*t + a4)*t + a3)*t + a2)*t + a1)*t + 0.5D0) + 0.5D0
fn_val = pi*a*w
GO TO 30

!                 A .LT. 0.25  OR  A .GT. 0.75

20 a = 0.25D0 + (0.75D0 - a)
n = n - 1
21 t = 16.d0*a*a
fn_val = (((((((((((((b13*t + b12)*t + b11)*t + b10)*t + b9)*t + b8)*t + &
         b7)*t + b6)*t + b5)*t + b4)*t + b3)*t + b2)*t + b1)*t + 0.5D0) + 0.5D0

!                        TERMINATION

30 IF (MOD(n,2) /= 0) fn_val = -fn_val
RETURN
END FUNCTION dcos1

FUNCTION drexp(x) RESULT(fn_val)
!-----------------------------------------------------------------------
!            EVALUATION OF THE FUNCTION EXP(X) - 1
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN)  :: x
REAL (dp)              :: fn_val

! Local variables
REAL (dp) :: e, w, z
REAL (dp) :: a0 = .248015873015873015873016D-04,   &
    a1 = -.344452080605731005808147D-05, a2 = .206664230430046597475413D-06,  &
    a3 = -.447300111094328162971036D-08, a4 = .114734027080634968083920D-11,  &
    b1 = -.249994190011341852652396D+00, b2 = .249987228833107957725728D-01,  &
    b3 = -.119037506846942249362528D-02, b4 = .228908693387350391768682D-04
REAL (dp) :: c1 = .1666666666666666666666666666666667D+00,   &
             c2 = .4166666666666666666666666666666667D-01,   &
             c3 = .8333333333333333333333333333333333D-02,   &
             c4 = .1388888888888888888888888888888889D-02,   &
             c5 = .1984126984126984126984126984126984D-03
!---------------------------
IF (ABS(x) <= 0.15_dp) THEN

!        Z IS A MINIMAX APPROXIMATION OF THE SERIES

!                C6 + C7*X + C8*X**2 + ....

!        THIS APPROXIMATION IS ACCURATE TO WITHIN
!        1 UNIT OF THE 23-RD SIGNIFICANT DIGIT.
!        THE RESULTING VALUE FOR W IS ACCURATE TO
!        WITHIN 1 UNIT OF THE 33-RD SIGNIFICANT DIGIT.

  z = ((((a4*x + a3)*x + a2)*x + a1)*x + a0) /  &
      ((((b4*x + b3)*x + b2)*x + b1)*x + 1._dp)
  w = ((((((z*x + c5)*x + c4)*x + c3)*x + c2)*x + c1)*x + 0.5_dp)*x + 1._dp
  fn_val = x * w
  RETURN
END IF

IF (x >= 0.0_dp) THEN
  e = EXP(x)
  fn_val = e * (0.5_dp + (0.5_dp - 1.0_dp/e))
  RETURN
END IF
IF (x >= -77._dp) THEN
  fn_val = (EXP(x)-0.5_dp) - 0.5_dp
  RETURN
END IF
fn_val = -1._dp
RETURN
END FUNCTION drexp



SUBROUTINE dcgama(mo, z, w)
!-----------------------------------------------------------------------

!        EVALUATION OF THE COMPLEX GAMMA AND LOGGAMMA FUNCTIONS

!                        ---------------

!     MO IS AN INTEGER.  Z AND W ARE INTERPRETED AS REAL (dp)
!     COMPLEX NUMBERS.  IT IS ASSUMED THAT Z(1) AND Z(2) ARE THE REAL
!     AND IMAGINARY PARTS OF THE COMPLEX NUMBER Z, AND THAT W(1) AND
!     W(2) ARE THE REAL AND IMAGINARY PARTS OF W.

!                 W = GAMMA(Z)       IF MO = 0
!                 W = LN(GAMMA(Z))   OTHERWISE

!-----------------------------------------------------------------------
INTEGER, INTENT(IN)        :: mo
COMPLEX (dp), INTENT(IN)   :: z
COMPLEX (dp), INTENT(OUT)  :: w

! Local variables
REAL (dp), PARAMETER :: c0(30)  &
        = (/ .8333333333333333333333333333333333333333D-01,  &
        -.2777777777777777777777777777777777777778D-02,  &
         .7936507936507936507936507936507936507937D-03,  &
        -.5952380952380952380952380952380952380952D-03,  &
         .8417508417508417508417508417508417508418D-03,  &
        -.1917526917526917526917526917526917526918D-02,  &
         .6410256410256410256410256410256410256410D-02,  &
        -.2955065359477124183006535947712418300654D-01,  &
         .1796443723688305731649384900158893966944D+00,  &
        -.1392432216905901116427432216905901116427D+01,  &
         .1340286404416839199447895100069013112491D+02,  &
        -.1568482846260020173063651324520889738281D+03,  &
         .2193103333333333333333333333333333333333D+04,  &
        -.3610877125372498935717326521924223073648D+05,  &
         .6914722688513130671083952507756734675533D+06,  &
        -.1523822153940741619228336495888678051866D+08,  &
         .3829007513914141414141414141414141414141D+09,  &
        -.1088226603578439108901514916552510537473D+11,  &
         .3473202837650022522522522522522522522523D+12,  &
        -.1236960214226927445425171034927132488108D+14,  &
         .4887880647930793350758151625180229021085D+15,  &
        -.2132033396091937389697505898213683855747D+17,  &
         .1021775296525700077565287628053585500394D+19,  &
        -.5357547217330020361082770919196920448485D+20,  &
         .3061578263704883415043151051329622758194D+22,  &
        -.1899991742639920405029371429306942902947D+24,  &
         .1276337403382883414923495137769782597654D+26,  &
        -.9252847176120416307230242348347622779519D+27,  &
         .7218822595185610297836050187301637922490D+29,  &
        -.6045183405995856967743148238754547286066D+31 /),  &
        dlpi = 1.144729885849400174143427351353058711647_dp,  &
        hl2p =  .9189385332046727417803297364056176398614_dp,  &
        pi = 3.141592653589793238462643383279502884197_dp,  &
        pi2 = 6.283185307179586476925286766559005768394_dp
REAL (dp) :: a, a1, a2, c, cn, cut, d, eps, et, e2t, h1, h2, q1, q2, s, sn,  &
             s1, s2, t, t1, t2, u, u1, u2, v1, v2, w1, w2, x, y, y2
INTEGER   :: j, k, max, n, nm1
!---------------------------
!     DLPI = LOG(PI)
!     HL2P = 0.5 * LOG(2*PI)
!---------------------------

!     ****** MAX AND EPS ARE MACHINE DEPENDENT CONSTANTS.
!            MAX IS THE LARGEST POSITIVE INTEGER THAT MAY
!            BE USED, AND EPS IS THE SMALLEST NUMBER SUCH
!            THAT  1._dp + EPS > 1._dp.

!                      MAX = IPMPAR(3)
max = HUGE(3)
eps = EPSILON(1.0_dp)

!---------------------------
x = REAL(z, KIND=dp)
y = AIMAG(z)
IF (x < 0._dp) THEN
!-----------------------------------------------------------------------
!            CASE WHEN THE REAL PART OF Z IS NEGATIVE
!-----------------------------------------------------------------------
  y = ABS(y)
  t = -pi * y
  et = EXP(t)
  e2t = et * et

!     SET  A1 = (1 + E2T)/2  AND  A2 = (1 - E2T)/2

  a1 = 0.5_dp * (1._dp+e2t)
  t2 = t + t
  IF (t2 >= -0.15_dp) THEN
    a2 = -0.5_dp * drexp(t2)
  ELSE
    a2 = 0.5_dp * (0.5_dp+(0.5_dp-e2t))
  END IF

!     COMPUTE SIN(PI*X) AND COS(PI*X)

  u = MAX
  IF (ABS(x) >= MIN(u,1._dp/eps)) GO TO 80
  k = ABS(x)
  u = x + k
  k = MOD(k,2)
  IF (u <= -0.5_dp) THEN
    u = 0.5_dp + (0.5_dp+u)
    k = k + 1
  END IF
  u = pi * u
  sn = SIN(u)
  cn = COS(u)
  IF (k == 1) THEN
    sn = -sn
    cn = -cn
  END IF

!     SET  H1 + H2*I  TO  PI/SIN(PI*Z)  OR  LOG(PI/SIN(PI*Z))

  a1 = sn * a1
  a2 = cn * a2
  a = a1 * a1 + a2 * a2
  IF (a == 0._dp) GO TO 80
  IF (mo == 0) THEN

    h1 = a1 / a
    h2 = -a2 / a
    c = pi * et
    h1 = c * h1
    h2 = c * h2
  ELSE

    h1 = (dlpi+t) - 0.5_dp * LOG(a)
    h2 = -ATAN2(a2,a1)
  END IF
  IF (AIMAG(z) >= 0._dp) THEN
    x = 1.0 - x
    y = -y
  ELSE
    h2 = -h2
    x = 1.0 - x
  END IF
END IF
!-----------------------------------------------------------------------
!           CASE WHEN THE REAL PART OF Z IS NONNEGATIVE
!-----------------------------------------------------------------------
w1 = 0._dp
w2 = 0._dp
n = 0
t = x
y2 = y * y
a = t * t + y2
cut = 225._dp
IF (eps > 1.d-30) cut = 144._dp
IF (eps > 1.d-20) cut = 64._dp
IF (a < cut) THEN
  IF (a == 0._dp) GO TO 80
  10 n = n + 1
  t = t + 1._dp
  a = t * t + y2
  IF (a < cut) GO TO 10

!     LET S1 + S2*I BE THE PRODUCT OF THE TERMS (Z+J)/(Z+N)

  u1 = (x*t+y2) / a
  u2 = y / a
  s1 = u1
  s2 = n * u2
  IF (n >= 2) THEN
    u = t / a
    nm1 = n - 1
    DO j = 1, nm1
      v1 = u1 + j * u
      v2 = (n-j) * u2
      c = s1 * v1 - s2 * v2
      d = s1 * v2 + s2 * v1
      s1 = c
      s2 = d
    END DO
  END IF

!     SET  W1 + W2*I = LOG(S1 + S2*I)  WHEN MO IS NONZERO

  s = s1 * s1 + s2 * s2
  IF (mo /= 0) THEN
    w1 = 0.5_dp * LOG(s)
    w2 = ATAN2(s2,s1)
  END IF
END IF

!     SET  V1 + V2*I = (Z - 0.5) * LOG(Z + N) - Z

t1 = 0.5_dp * LOG(a) - 1._dp
t2 = ATAN2(y,t)
u = x - 0.5_dp
v1 = (u*t1-0.5_dp) - y * t2
v2 = u * t2 + y * t1

!     LET A1 + A2*I BE THE ASYMPTOTIC SUM

u1 = t / a
u2 = -y / a
q1 = u1 * u1 - u2 * u2
q2 = 2._dp * u1 * u2
a1 = 0._dp
a2 = 0._dp
DO j = 1, 30
  t1 = a1
  t2 = a2
  a1 = a1 + c0(j) * u1
  a2 = a2 + c0(j) * u2
  IF (a1 == t1) THEN
    IF (a2 == t2) GO TO 40
  END IF
  t1 = u1 * q1 - u2 * q2
  t2 = u1 * q2 + u2 * q1
  u1 = t1
  u2 = t2
END DO
!-----------------------------------------------------------------------
!                 GATHERING TOGETHER THE RESULTS
!-----------------------------------------------------------------------
40 w1 = (((a1+hl2p)-w1)+v1) - n
w2 = (a2-w2) + v2
IF (REAL(z, KIND=dp) < 0.0_dp) GO TO 60
IF (mo == 0) THEN

!     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO = 0

  a = EXP(w1)
  w1 = a * COS(w2)
  w2 = a * SIN(w2)
  IF (n == 0) GO TO 70
  c = (s1*w1+s2*w2) / s
  d = (s1*w2-s2*w1) / s
  w1 = c
  w2 = d
  GO TO 70
END IF

!     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO IS NONZERO.
!     THE ANGLE W2 IS REDUCED TO THE INTERVAL -PI < W2 <= PI.

50 IF (w2 <= pi) THEN
  k = 0.5_dp - w2 / pi2
  w2 = w2 + pi2 * k
  GO TO 70
END IF
k = w2 / pi2 - 0.5_dp
u = k + 1
w2 = w2 - pi2 * u
IF (w2 <= -pi) w2 = pi
GO TO 70

!     CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO IS NONZERO

60 IF (mo /= 0) THEN
  w1 = h1 - w1
  w2 = h2 - w2
  GO TO 50
END IF

!     CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO = 0

a = EXP(-w1)
t1 = a * COS(-w2)
t2 = a * SIN(-w2)
w1 = h1 * t1 - h2 * t2
w2 = h1 * t2 + h2 * t1
IF (n /= 0) THEN
  c = w1 * s1 - w2 * s2
  d = w1 * s2 + w2 * s1
  w1 = c
  w2 = d
END IF

!     TERMINATION

70 w = CMPLX(w1, w2, KIND=dp)
RETURN
!-----------------------------------------------------------------------
!             THE REQUESTED VALUE CANNOT BE COMPUTED
!-----------------------------------------------------------------------
80 w = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
RETURN
END SUBROUTINE dcgama

FUNCTION cgam0(z) RESULT(fn_val)
!-----------------------------------------------------------------------
!          EVALUATION OF 1/GAMMA(1 + Z)  FOR ABS(Z) < 1.0
!-----------------------------------------------------------------------

COMPLEX (dp), INTENT(IN)  :: z
COMPLEX (dp)              :: fn_val

COMPLEX (dp)  :: w
INTEGER       :: i, k, n
!-----------------------
REAL (dp)  :: x, y
REAL, PARAMETER :: a(25) = (/ .577215664901533_dp, -.655878071520254_dp,  &
        -.420026350340952D-01, .166538611382291_dp, -.421977345555443D-01,  &
        -.962197152787697D-02, .721894324666310D-02, -.116516759185907D-02,  &
        -.215241674114951D-03, .128050282388116D-03, -.201348547807882D-04,  &
        -.125049348214267D-0, .113302723198170D-05, -.205633841697761D-0,  &
        .611609510448142D-08, .500200764446922D-08, -.118127457048702D-08,  &
        .104342671169110D-09, .778226343990507D-11, -.369680561864221D-11,  &
        .510037028745448D-12, -.205832605356651D-13, -.534812253942302D-14,  &
        .122677862823826D-14, -.118125930169746D-15 /)
!-----------------------
n = 25
x = REAL(z, KIND=dp)
y = AIMAG(z)
IF (x*x+y*y <= 0.04D0) n = 14

k = n
w = a(n)
DO  i = 2, n
  k = k - 1
  w = a(k) + z * w
END DO
fn_val = 1.0D0 + z * w
RETURN
END FUNCTION cgam0



FUNCTION dgamma(a) RESULT(fn_val)
!-----------------------------------------------------------------------

!                EVALUATION OF THE GAMMA FUNCTION FOR
!                     REAL (dp) ARGUMENTS

!                           -----------

!     DGAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
!     BE COMPUTED.

!-----------------------------------------------------------------------
!     WRITTEN BY ALFRED H. MORRIS, JR.
!          NAVAL SURFACE WEAPONS CENTER
!          DAHLGREN, VIRGINIA
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN) :: a
REAL (dp)             :: fn_val

! Local variables
REAL (dp), PARAMETER :: d = 0.41893853320467274178032973640562_dp,  &
                        pi = 3.14159265358979323846264338327950_dp
REAL (dp) :: s, t, x, w
INTEGER   :: j, n
!-----------------------------------------------------------------------
!     D = 0.5*(LN(2*PI) - 1)
!-----------------------------------------------------------------------
fn_val = 0.0_dp
x = a
IF (ABS(a) <= 20._dp) THEN
!-----------------------------------------------------------------------
!             EVALUATION OF DGAMMA(A) FOR ABS(A) <= 20
!-----------------------------------------------------------------------
  t = 1.0_dp
  n = x
  n = n - 1

!     LET T BE THE PRODUCT OF A-J WHEN A >= 2

  IF (n < 0) THEN
    GO TO 40
  ELSE IF (n == 0) THEN
    GO TO 30
  END IF

  DO j = 1, n
    x = x - 1._dp
    t = x * t
  END DO
  30 x = x - 1._dp
  GO TO 60

!     LET T BE THE PRODUCT OF A+J WHEN A < 1

  40 t = a
  IF (a <= 0._dp) THEN
    n = -n - 1
    IF (n /= 0) THEN
      DO j = 1, n
        x = x + 1._dp
        t = x * t
      END DO
    END IF
    x = (x+0.5_dp) + 0.5_dp
    t = x * t
    IF (t == 0._dp) RETURN
  END IF

!     THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
!     CODE MAY BE OMITTED IF DESIRED.

  IF (ABS(t) < 1.d-33) THEN
    IF (ABS(t)*HUGE(1.0_dp) <= 1.000000001_dp) RETURN
    fn_val = 1._dp / t
    RETURN
  END IF

!     COMPUTE DGAMMA(1 + X) FOR 0 <= X < 1

  60 fn_val = 1._dp / (1._dp + dgam1(x))

!     TERMINATION

  IF (a >= 1._dp) THEN
    fn_val = fn_val * t
    RETURN
  END IF
  fn_val = fn_val / t
  RETURN
END IF
!-----------------------------------------------------------------------
!           EVALUATION OF DGAMMA(A) FOR ABS(A) > 20
!-----------------------------------------------------------------------
IF (ABS(a) >= 1.d3) RETURN
IF (a <= 0._dp) THEN
  s = dsin1(a) / pi
  IF (s == 0._dp) RETURN
  x = -a
END IF

!     COMPUTE THE MODIFIED ASYMPTOTIC SUM

w = dpdel(x)

!     FINAL ASSEMBLY

w = (d+w) + (x-0.5_dp) * (LOG(x)-1._dp)
IF (w > dxparg(0)) RETURN
fn_val = EXP(w)
IF (a < 0._dp) fn_val = (1._dp/(fn_val*s)) / x

RETURN
END FUNCTION dgamma



FUNCTION glog(x) RESULT(fn_val)
!     -------------------
!     EVALUATION OF LN(X) FOR X >= 15
!     -------------------
REAL (dp), INTENT(IN) :: x
REAL (dp)             :: fn_val

! Local variables
REAL (dp), PARAMETER :: c1 = .286228750476730_dp, c2 = .399999628131494+dp,  &
                        c3 = .666666666752663_dp
REAL (dp), PARAMETER :: w(163) = (/ .270805020110221007D+01, .277258872223978124D+01,  &
  .283321334405621608D+01, .289037175789616469D+01, .294443897916644046D+01, .299573227355399099D+01, .304452243772342300D+01, &
  .309104245335831585D+01, .313549421592914969D+01, .317805383034794562D+01, .321887582486820075D+01, .325809653802148205D+01, &
  .329583686600432907D+01, .333220451017520392D+01, .336729582998647403D+01, .340119738166215538D+01, .343398720448514625D+01, &
  .346573590279972655D+01, .349650756146648024D+01, .352636052461616139D+01, .355534806148941368D+01, .358351893845611000D+01, &
  .361091791264422444D+01, .363758615972638577D+01, .366356164612964643D+01, .368887945411393630D+01, .371357206670430780D+01, &
  .373766961828336831D+01, .376120011569356242D+01, .378418963391826116D+01, .380666248977031976D+01, .382864139648909500D+01, &
  .385014760171005859D+01, .387120101090789093D+01, .389182029811062661D+01, .391202300542814606D+01, .393182563272432577D+01, &
  .395124371858142735D+01, .397029191355212183D+01, .398898404656427438D+01, .400733318523247092D+01, .402535169073514923D+01, &
  .404305126783455015D+01, .406044301054641934D+01, .407753744390571945D+01, .409434456222210068D+01, .411087386417331125D+01, &
  .412713438504509156D+01, .414313472639153269D+01, .415888308335967186D+01, .417438726989563711D+01, .418965474202642554D+01, &
  .420469261939096606D+01, .421950770517610670D+01, .423410650459725938D+01, .424849524204935899D+01, .426267987704131542D+01, &
  .427666611901605531D+01, .429045944114839113D+01, .430406509320416975D+01, .431748811353631044D+01, .433073334028633108D+01, &
  .434380542185368385D+01, .435670882668959174D+01, .436944785246702149D+01, .438202663467388161D+01, .439444915467243877D+01, &
  .440671924726425311D+01, .441884060779659792D+01, .443081679884331362D+01, .444265125649031645D+01, .445434729625350773D+01, &
  .446590811865458372D+01, .447733681447820647D+01, .448863636973213984D+01, .449980967033026507D+01, .451085950651685004D+01, &
  .452178857704904031D+01, .453259949315325594D+01, .454329478227000390D+01, .455387689160054083D+01, .456434819146783624D+01, &
  .457471097850338282D+01, .458496747867057192D+01, .459511985013458993D+01, .460517018598809137D+01, .461512051684125945D+01, &
  .462497281328427108D+01, .463472898822963577D+01, .464439089914137266D+01, .465396035015752337D+01, .466343909411206714D+01, &
  .467282883446190617D+01, .468213122712421969D+01, .469134788222914370D+01, .470048036579241623D+01, .470953020131233414D+01, &
  .471849887129509454D+01, .472738781871234057D+01, .473619844839449546D+01, .474493212836325007D+01, .475359019110636465D+01, &
  .476217393479775612D+01, .477068462446566476D+01, .477912349311152939D+01, .478749174278204599D+01, .479579054559674109D+01, &
  .480402104473325656D+01, .481218435537241750D+01, .482028156560503686D+01, .482831373730230112D+01, .483628190695147800D+01, &
  .484418708645859127D+01, .485203026391961717D+01, .485981240436167211D+01, .486753445045558242D+01, .487519732320115154D+01, &
  .488280192258637085D+01, .489034912822175377D+01, .489783979995091137D+01, .490527477843842945D+01, .491265488573605201D+01, &
  .491998092582812492D+01, .492725368515720469D+01, .493447393313069176D+01, .494164242260930430D+01, .494875989037816828D+01, &
  .495582705760126073D+01, .496284463025990728D+01, .496981329957600062D+01, .497673374242057440D+01, .498360662170833644D+01, &
  .499043258677873630D+01, .499721227376411506D+01, .500394630594545914D+01, .501063529409625575D+01, .501727983681492433D+01, &
  .502388052084627639D+01, .503043792139243546D+01, .503695260241362916D+01, .504342511691924662D+01, .504985600724953705D+01, &
  .505624580534830806D+01, .506259503302696680D+01, .506890420222023153D+01, .507517381523382692D+01, .508140436498446300D+01, &
  .508759633523238407D+01, .509375020080676233D+01, .509986642782419842D+01, .510594547390058061D+01, .511198778835654323D+01, &
  .511799381241675511D+01, .512396397940325892D+01, .512989871492307347D+01, .513579843705026176D+01, .514166355650265984D+01, &
  .514749447681345304D+01, .515329159449777895D+01, .515905529921452903D+01, .516478597392351405D+01, .517048399503815178D+01, &
  .517614973257382914D+01 /)
REAL (dp) :: t, t2, z
INTEGER   :: n
!     -------------------

IF (x < 178.0_dp) THEN
  n = x
  t = (x-n) / (x+n)
  t2 = t * t
  z = (((c1*t2 + c2)*t2 + c3)*t2 + 2.0) * t
  fn_val = w(n-14) + z
  RETURN
END IF

fn_val = LOG(x)
RETURN
END FUNCTION glog

SUBROUTINE cbsslj(z,cnu,w)
!-----------------------------------------------------------------------

!         EVALUATION OF THE COMPLEX BESSEL FUNCTION J   (Z)
!                                                    CNU
!-----------------------------------------------------------------------

!     WRITTEN BY
!         ANDREW H. VAN TUYL AND ALFRED H. MORRIS, JR.
!         NAVAL SURFACE WARFARE CENTER
!         OCTOBER, 1991

!     A MODIFICATION OF THE PROCEDURE DEVELOPED BY ALLEN V. HERSHEY
!     (NAVAL SURFACE WARFARE CENTER) IN 1978 FOR HANDLING THE DEBYE
!     APPROXIMATION IS EMPLOYED.

!-----------------------------------------------------------------------

COMPLEX (dp), INTENT(IN)   :: z
COMPLEX (dp), INTENT(IN)   :: cnu
COMPLEX (dp), INTENT(OUT)  :: w

COMPLEX (dp)  :: c, nu, s, sm1, sm2, t, tsc, w0, w1, zn, zz
!-----------------------
REAL (dp) :: a, cn1, cn2, e, fn
REAL (dp) :: pn, qm, qn, qnp1
REAL (dp) :: r, rn2, r2, sn, t1, t2
REAL (dp) :: u, v, x, y
INTEGER   :: i, k, m, n
REAL (dp), PARAMETER  :: pi = 3.141592653589793238462643383279502884197_dp
!-----------------------
x = REAL(z, KIND=dp)
y = AIMAG(z)
r = cpabs(x,y)
cn1 = REAL(cnu, KIND=dp)
cn2 = AIMAG(cnu)
rn2 = cn1 * cn1 + cn2 * cn2
pn = INT(cn1)
fn = cn1 - pn
sn = 1.0_dp

!          CALCULATION WHEN ORDER IS AN INTEGER

IF (fn == 0.0_dp .AND. cn2 == 0.0_dp) THEN
  n = pn
  pn = ABS(pn)
  cn1 = pn
  IF (n < 0 .AND. n /= (n/2)*2) sn = -1.0_dp
END IF

!          SELECTION OF METHOD

IF (r > 17.5D0) THEN
  IF (r > 17.5D0 + 0.5D0*rn2) GO TO 10
  GO TO 20
END IF

!          USE MACLAURIN EXPANSION AND RECURSION

IF (cn1 < 0.0D0) THEN
  qn = -1.25D0 * (r + 0.5D0*ABS(cn2) - ABS(y-0.5D0*cn2))
  IF (cn1 < qn) THEN
    qn = 1.25D0 * (r - MAX(1.2D0*r,ABS(y-cn2)))
    IF (cn1 < qn) THEN
      qn = MIN(pn, REAL(-INT(1.25D0*(r-ABS(cn2))), KIND=dp))
      GO TO 60
    END IF
  END IF
END IF

r2 = r * r
qm = 0.0625D0 * r2 * r2 - cn2 * cn2
qn = MAX(pn, REAL(INT(SQRT(MAX(0.0D0,qm))), KIND=dp))
GO TO 60

!          USE ASYMPTOTIC EXPANSION

10 CALL cbja(z,cnu,w)
RETURN

!          CALCULATION FOR 17.5 < ABS(Z) <= 17.5 + 0.5*ABS(CNU)**2

20 n = 0
IF (ABS(cn2) < 0.8D0*ABS(y)) THEN
  qm = -1.25D0 * (r + 0.5D0*ABS(cn2) - ABS(y-0.5D0*cn2))
  IF (cn1 < qm) THEN
    qm = 1.25D0 * (r - MAX(1.2D0*r, ABS(y-cn2)))
    IF (cn1 < qm) n = 1
  END IF
END IF

qn = pn
a = 4.d-3 * r * r
zz = z
IF (x < 0.0D0) zz = -z

!          CALCULATION OF ZONE OF EXCLUSION OF DEBYE APPROXIMATION

30 nu = CMPLX(qn+fn,cn2, KIND=dp)
zn = nu / z
t2 = AIMAG(zn) * AIMAG(zn)
u = 1.0D0 - REAL(zn, KIND=dp)
t1 = u * u + t2
u = 1.0D0 + DBLE(zn)
t2 = u * u + t2
u = t1 * t2
v = a * u / (t1*t1 + t2*t2)
IF (u*v*v <= 1.0D0) THEN
  
!          THE ARGUMENT LIES INSIDE THE ZONE OF EXCLUSION
  
  qn = qn + 1.0D0
  IF (n == 0) GO TO 30
  
!          USE MACLAURIN EXPANSION WITH FORWARD RECURRENCE
  
  qn = MIN(pn, REAL(-INT(1.25D0*(r-ABS(cn2))), KIND=dp))
ELSE
  
!          USE BACKWARD RECURRENCE STARTING FROM THE ASYMPTOTIC EXPANSION
  
  qnp1 = qn + 1.0D0
  IF (ABS(qn) < ABS(pn)) THEN
    IF (r >= 17.5D0 + 0.5D0*(qnp1*qnp1 + cn2*cn2)) THEN
      
      nu = CMPLX(qn+fn,cn2, KIND=dp)
      CALL cbja(zz,nu,sm1)
      nu = CMPLX(qnp1+fn,cn2, KIND=dp)
      CALL cbja(zz,nu,sm2)
      GO TO 40
    END IF
  END IF
  
!          USE BACKWARD RECURRENCE STARTING FROM THE DEBYE APPROXIMATION
  
  nu = CMPLX(qn+fn,cn2, KIND=dp)
  CALL cbdb(zz,nu,fn,sm1)
  IF (qn == pn) GO TO 50
  nu = CMPLX(qnp1+fn,cn2, KIND=dp)
  CALL cbdb(zz,nu,fn,sm2)
  
  40 nu = CMPLX(qn+fn,cn2, KIND=dp)
  tsc = 2.0D0 * nu * sm1 / zz - sm2
  sm2 = sm1
  sm1 = tsc
  qn = qn - 1.0D0
  IF (qn /= pn) GO TO 40
  
  50 w = sm1
  IF (sn < 0.0D0) w = -w
  IF (x >= 0.0D0) RETURN
  
  nu = pi * CMPLX(-cn2,cn1, KIND=dp)
  IF (y < 0.0D0) nu = -nu
  w = EXP(nu) * w
  RETURN
END IF

!          USE MACLAURIN EXPANSION WITH FORWARD OR BACKWARD RECURRENCE.

60 m = qn - pn
IF (ABS(m) <= 1) THEN
  nu = CMPLX(cn1,cn2, KIND=dp)
  CALL cbjm(z,nu,w)
ELSE
  nu = CMPLX(qn+fn,cn2, KIND=dp)
  CALL cbjm(z,nu,w1)
  w0 = 0.25D0 * z * z
  IF (m <= 0) THEN
    
!          FORWARD RECURRENCE
    
    m = ABS(m)
    nu = nu + 1.0D0
    CALL cbjm(z,nu,w)
    DO  i = 2, m
      c = nu * (nu+1.0D0)
      t = (c/w0) * (w-w1)
      w1 = w
      w = t
      nu = nu + 1.0D0
    END DO
  ELSE
    
!          BACKWARD RECURRENCE
    
    nu = nu - 1.0D0
    CALL cbjm(z,nu,w)
    DO  i = 2, m
      c = nu * (nu+1.0D0)
      t = (w0/c) * w1
      w1 = w
      w = w - t
      nu = nu - 1.0D0
    END DO
  END IF
END IF

!          FINAL ASSEMBLY

IF (fn == 0.0D0 .AND. cn2 == 0.0D0) THEN
  k = pn
  IF (k == 0) RETURN
  e = sn / dgamma(pn+1.0D0)
  w = e * w * (0.5D0*z) ** k
  RETURN
END IF

s = cnu * LOG(0.5D0*z)
w = EXP(s) * w
IF (rn2 <= 0.81D0) THEN
  w = w * cgam0(cnu)
  RETURN
END IF
CALL dcgama(0,cnu,t)
w = cdiv(w,cnu*t)
RETURN
END SUBROUTINE cbsslj



SUBROUTINE cbjm(z,cnu,w)
!-----------------------------------------------------------------------

!       COMPUTATION OF  (Z/2)**(-CNU) * GAMMA(CNU + 1) * J(CNU,Z)

!                           -----------------

!     THE MACLAURIN EXPANSION IS USED.  IT IS ASSUMED THAT CNU IS NOT
!     A NEGATIVE INTEGER.

!-----------------------------------------------------------------------

COMPLEX (dp), INTENT(IN)   :: z
COMPLEX (dp), INTENT(IN)   :: cnu
COMPLEX (dp), INTENT(OUT)  :: w

COMPLEX (dp)  :: nu, nup1, p, s, sn, t, ti
!--------------------------
REAL (dp)  :: a, a0, eps, inu, m, rnu
INTEGER    :: i, imin, k, km1, km2

!--------------------------

!     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
!            SMALLEST NUMBER SUCH THAT 1.0 + EPS .GT. 1.0 .

eps = EPSILON(0.0_dp)

!--------------------------
s = -0.25D0 * (z*z)
nu = cnu
rnu = REAL(nu, KIND=dp)
inu = AIMAG(nu)
a = 0.5D0 + (0.5D0+rnu)
nup1 = CMPLX(a,inu, KIND=dp)

IF (a > 0.0D0) THEN
  m = 1.0D0
  t = s / nup1
  w = 1.0D0 + t
ELSE
  
!     ADD 1.0 AND THE FIRST K-1 TERMS
  
  k = INT(-a) + 2
  km1 = k - 1
  w = (1.0D0,0.0D0)
  t = w
  DO  i = 1, km1
    m = i
    t = t * (s/(m*(nu+m)))
    w = w + t
    IF (anorm(t) <= eps*anorm(w)) GO TO 20
  END DO
  GO TO 50
  
!     CHECK IF THE (K-1)-ST AND K-TH TERMS CAN BE IGNORED.
!     IF SO THEN THE SUMMATION IS COMPLETE.
  
  20 IF (i /= km1) THEN
    imin = i + 1
    IF (imin < k-5) THEN
      ti = t
      
      m = km1
      t = s / (nu+m)
      a0 = anorm(t) / m
      t = t * (s/(nu+(m+1.0D0)))
      a = anorm(t) / (m*(m+1.0D0))
      a = MAX(a,a0)
      
      t = (1.0D0,0.0D0)
      km2 = k - 2
      DO  i = imin, km2
        m = i
        t = t * (s/(m*(nu+m)))
        IF (a*anorm(t) < 0.5D0) RETURN
      END DO
      t = t * ti
      imin = km2
    END IF
    
!     ADD THE (K-1)-ST TERM
    
    a = 1.0D0
    p = (1.0D0,0.0D0)
    sn = p
    DO  i = imin, km1
      m = i
      a = a * m
      p = p * (nu+m)
      sn = s * sn
    END DO
    t = t * (cdiv(sn,p)/a)
    w = w + t
  END IF
END IF

!     ADD THE REMAINING TERMS

50 m = m + 1.0D0
t = t * (s/(m*(nu+m)))
w = w + t
IF (anorm(t) > eps*anorm(w)) GO TO 50

RETURN
END SUBROUTINE cbjm