求助-复变量的贝塞尔第一类和第二类函数的Fortran子程序

风平老涡

[Fortran] 纯文本查看 复制代码

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

风平老涡

jlg1206 发表于 2021-1-31 18:57
翻墙搜了很多，但都不能运行，编译不能通过。我也查了很多书带的子程序，但是可能代码太老也是没法编译 ...

[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

风平老涡

jlg1206 发表于 2021-1-31 18:57
翻墙搜了很多，但都不能运行，编译不能通过。我也查了很多书带的子程序，但是可能代码太老也是没法编译 ...

由于字数的限置，我把整个程序分成几个回帖。

jlg1206

风平老涡 发表于 2021-1-26 00:34
去网上搜一搜，应该有很多的。

翻墙搜了很多，但都不能运行，编译不能通过。我也查了很多书带的子程序，但是可能代码太老也是没法编译

风平老涡

去网上搜一搜，应该有很多的。