mairya.f90

数值计算和数值分析在Fortran下的特殊函数库,是数值计算的必备

MODULE airya_func
 
! From the book "Computation of Special Functions"
!      by Shanjie Zhang and Jianming Jin
!   Copyright 1996 by John Wiley & Sons, Inc.
! The authors state:
!   "However, we give permission to the reader who purchases this book
!    to incorporate any of these programs into his or her programs
!    provided that the copyright is acknowledged."
 
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
 
CONTAINS


SUBROUTINE airya(x, ai, bi, ad, bd)

!    ======================================================
!    Purpose: Compute Airy functions and their derivatives
!    Input:   x  --- Argument of Airy function
!    Output:  AI --- Ai(x)
!             BI --- Bi(x)
!             AD --- Ai'(x)
!             BD --- Bi'(x)
!    Routine called:
!             AJYIK for computing Jv(x), Yv(x), Iv(x) and
!             Kv(x) with v=1/3 and 2/3
!    ======================================================

REAL (dp), INTENT(IN)      :: x
REAL (dp), INTENT(OUT)     :: ai
REAL (dp), INTENT(OUT)     :: bi
REAL (dp), INTENT(OUT)     :: ad
REAL (dp), INTENT(OUT)     :: bd

REAL (dp)  :: xa, pir, c1, c2, sr3, z, xq, vj1, vj2, vy1, vy2, vi1, vi2,  &
              vk1, vk2

xa = ABS(x)
pir = 0.318309886183891_dp
c1 = 0.355028053887817_dp
c2 = 0.258819403792807_dp
sr3 = 1.732050807568877_dp
z = xa ** 1.5 / 1.5_dp
xq = SQRT(xa)
CALL ajyik(z, vj1, vj2, vy1, vy2, vi1, vi2, vk1, vk2)
IF (x == 0.0_dp) THEN
  ai = c1
  bi = sr3 * c1
  ad = -c2
  bd = sr3 * c2
ELSE IF (x > 0.0_dp) THEN
  ai = pir * xq / sr3 * vk1
  bi = xq * (pir*vk1 + 2.0_dp/sr3*vi1)
  ad = -xa / sr3 * pir * vk2
  bd = xa * (pir*vk2 + 2.0_dp/sr3*vi2)
ELSE
  ai = 0.5_dp * xq * (vj1 - vy1/sr3)
  bi = -0.5_dp * xq * (vj1/sr3 + vy1)
  ad = 0.5_dp * xa * (vj2 + vy2/sr3)
  bd = 0.5_dp * xa * (vj2/sr3-vy2)
END IF
RETURN
END SUBROUTINE airya


SUBROUTINE ajyik(x, vj1, vj2, vy1, vy2, vi1, vi2, vk1, vk2)

!       =======================================================
!       Purpose: Compute Bessel functions Jv(x) and Yv(x),
!                and modified Bessel functions Iv(x) and
!                Kv(x), and their derivatives with v=1/3,2/3
!       Input :  x --- Argument of Jv(x),Yv(x),Iv(x) and
!                      Kv(x) ( x