|
|
C - RELATE PLASMA Z-FUNCTION TO THE COMPLENETARY ERROR FUNCTION AND USE A
|
|
|
C LIBRARY TO CALCULATE THE COMPLEMENTARY ERROR FUNCTION
|
|
|
|
|
|
complex*16 FUNCTION Zz(ZARG)
|
|
|
IMPLICIT NONE
|
|
|
COMPLEX*16 I,ZARG
|
|
|
REAL*8 XI,YI,U,V,PI
|
|
|
PARAMETER(I=(0.0D0,1.0D0),PI=3.141592653589793d0)
|
|
|
LOGICAL FLAG
|
|
|
|
|
|
c write(*,*)"zeta",zarg
|
|
|
|
|
|
XI=dREAL(ZARG)
|
|
|
YI=diMAG(ZARG)
|
|
|
CALL WOFZ(XI, YI, U, V, FLAG)
|
|
|
Zz=DSQRT(PI)*I*(U+I*V)
|
|
|
|
|
|
RETURN
|
|
|
END
|
|
|
|
|
|
C ALGORITHM 680, COLLECTED ALGORITHMS FROM ACM.
|
|
|
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
|
|
|
C VOL. 16, NO. 1, PP. 47.
|
|
|
SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
|
|
|
C
|
|
|
C GIVEN A COMPLEX NUMBER Z = (XI,YI), THIS SUBROUTINE COMPUTES
|
|
|
C THE VALUE OF THE FADDEEVA-FUNCTION W(Z) = EXP(-Z**2)*ERFC(-I*Z),
|
|
|
C WHERE ERFC IS THE COMPLEX COMPLEMENTARY ERROR-FUNCTION AND I
|
|
|
C MEANS SQRT(-1).
|
|
|
C THE ACCURACY OF THE ALGORITHM FOR Z IN THE 1ST AND 2ND QUADRANT
|
|
|
C IS 14 SIGNIFICANT DIGITS; IN THE 3RD AND 4TH IT IS 13 SIGNIFICANT
|
|
|
C DIGITS OUTSIDE A CIRCULAR REGION WITH RADIUS 0.126 AROUND A ZERO
|
|
|
C OF THE FUNCTION.
|
|
|
C ALL REAL VARIABLES IN THE PROGRAM ARE DOUBLE PRECISION.
|
|
|
C
|
|
|
C
|
|
|
C THE CODE CONTAINS A FEW COMPILER-DEPENDENT PARAMETERS :
|
|
|
C RMAXREAL = THE MAXIMUM VALUE OF RMAXREAL EQUALS THE ROOT OF
|
|
|
C RMAX = THE LARGEST NUMBER WHICH CAN STILL BE
|
|
|
C IMPLEMENTED ON THE COMPUTER IN DOUBLE PRECISION
|
|
|
C FLOATING-POINT ARITHMETIC
|
|
|
C RMAXEXP = LN(RMAX) - LN(2)
|
|
|
C RMAXGONI = THE LARGEST POSSIBLE ARGUMENT OF A DOUBLE PRECISION
|
|
|
C GONIOMETRIC FUNCTION (DCOS, DSIN, ...)
|
|
|
C THE REASON WHY THESE PARAMETERS ARE NEEDED AS THEY ARE DEFINED WILL
|
|
|
C BE EXPLAINED IN THE CODE BY MEANS OF COMMENTS
|
|
|
C
|
|
|
C
|
|
|
C PARAMETER LIST
|
|
|
C XI = REAL PART OF Z
|
|
|
C YI = IMAGINARY PART OF Z
|
|
|
C U = REAL PART OF W(Z)
|
|
|
C V = IMAGINARY PART OF W(Z)
|
|
|
C FLAG = AN ERROR FLAG INDICATING WHETHER OVERFLOW WILL
|
|
|
C OCCUR OR NOT; TYPE LOGICAL;
|
|
|
C THE VALUES OF THIS VARIABLE HAVE THE FOLLOWING
|
|
|
C MEANING :
|
|
|
C FLAG=.FALSE. : NO ERROR CONDITION
|
|
|
C FLAG=.TRUE. : OVERFLOW WILL OCCUR, THE ROUTINE
|
|
|
C BECOMES INACTIVE
|
|
|
C XI, YI ARE THE INPUT-PARAMETERS
|
|
|
C U, V, FLAG ARE THE OUTPUT-PARAMETERS
|
|
|
C
|
|
|
C FURTHERMORE THE PARAMETER FACTOR EQUALS 2/SQRT(PI)
|
|
|
C
|
|
|
C THE ROUTINE IS NOT UNDERFLOW-PROTECTED BUT ANY VARIABLE CAN BE
|
|
|
C PUT TO 0 UPON UNDERFLOW;
|
|
|
C
|
|
|
C REFERENCE - GPM POPPE, CMJ WIJERS; MORE EFFICIENT COMPUTATION OF
|
|
|
C THE COMPLEX ERROR-FUNCTION, ACM TRANS. MATH. SOFTWARE.
|
|
|
C
|
|
|
*
|
|
|
*
|
|
|
*
|
|
|
*
|
|
|
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
|
|
|
*
|
|
|
LOGICAL A, B, FLAG
|
|
|
PARAMETER (FACTOR = 1.12837916709551257388D0,
|
|
|
* RMAXREAL = 0.5D+154,
|
|
|
* RMAXEXP = 708.503061461606D0,
|
|
|
* RMAXGONI = 3.53711887601422D+15)
|
|
|
*
|
|
|
FLAG = .FALSE.
|
|
|
*
|
|
|
XABS = DABS(XI)
|
|
|
YABS = DABS(YI)
|
|
|
X = XABS/6.3
|
|
|
Y = YABS/4.4
|
|
|
*
|
|
|
C
|
|
|
C THE FOLLOWING IF-STATEMENT PROTECTS
|
|
|
C QRHO = (X**2 + Y**2) AGAINST OVERFLOW
|
|
|
C
|
|
|
IF ((XABS.GT.RMAXREAL).OR.(YABS.GT.RMAXREAL)) GOTO 100
|
|
|
*
|
|
|
QRHO = X**2 + Y**2
|
|
|
*
|
|
|
XABSQ = XABS**2
|
|
|
XQUAD = XABSQ - YABS**2
|
|
|
YQUAD = 2*XABS*YABS
|
|
|
*
|
|
|
A = QRHO.LT.0.085264D0
|
|
|
*
|
|
|
IF (A) THEN
|
|
|
C
|
|
|
C IF (QRHO.LT.0.085264D0) THEN THE FADDEEVA-FUNCTION IS EVALUATED
|
|
|
C USING A POWER-SERIES (ABRAMOWITZ/STEGUN, EQUATION (7.1.5), P.297)
|
|
|
C N IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
|
|
|
C ACCURACY
|
|
|
C
|
|
|
QRHO = (1-0.85*Y)*DSQRT(QRHO)
|
|
|
N = IDNINT(6 + 72*QRHO)
|
|
|
J = 2*N+1
|
|
|
XSUM = 1.0/J
|
|
|
YSUM = 0.0D0
|
|
|
DO 10 I=N, 1, -1
|
|
|
J = J - 2
|
|
|
XAUX = (XSUM*XQUAD - YSUM*YQUAD)/I
|
|
|
YSUM = (XSUM*YQUAD + YSUM*XQUAD)/I
|
|
|
XSUM = XAUX + 1.0/J
|
|
|
10 CONTINUE
|
|
|
U1 = -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0
|
|
|
V1 = FACTOR*(XSUM*XABS - YSUM*YABS)
|
|
|
DAUX = DEXP(-XQUAD)
|
|
|
U2 = DAUX*DCOS(YQUAD)
|
|
|
V2 = -DAUX*DSIN(YQUAD)
|
|
|
*
|
|
|
U = U1*U2 - V1*V2
|
|
|
V = U1*V2 + V1*U2
|
|
|
*
|
|
|
ELSE
|
|
|
C
|
|
|
C IF (QRHO.GT.1.O) THEN W(Z) IS EVALUATED USING THE LAPLACE
|
|
|
C CONTINUED FRACTION
|
|
|
C NU IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
|
|
|
C ACCURACY
|
|
|
C
|
|
|
C IF ((QRHO.GT.0.085264D0).AND.(QRHO.LT.1.0)) THEN W(Z) IS EVALUATED
|
|
|
C BY A TRUNCATED TAYLOR EXPANSION, WHERE THE LAPLACE CONTINUED FRACTION
|
|
|
C IS USED TO CALCULATE THE DERIVATIVES OF W(Z)
|
|
|
C KAPN IS THE MINIMUM NUMBER OF TERMS IN THE TAYLOR EXPANSION NEEDED
|
|
|
C TO OBTAIN THE REQUIRED ACCURACY
|
|
|
C NU IS THE MINIMUM NUMBER OF TERMS OF THE CONTINUED FRACTION NEEDED
|
|
|
C TO CALCULATE THE DERIVATIVES WITH THE REQUIRED ACCURACY
|
|
|
C
|
|
|
*
|
|
|
IF (QRHO.GT.1.0) THEN
|
|
|
H = 0.0D0
|
|
|
KAPN = 0
|
|
|
QRHO = DSQRT(QRHO)
|
|
|
NU = IDINT(3 + (1442/(26*QRHO+77)))
|
|
|
ELSE
|
|
|
QRHO = (1-Y)*DSQRT(1-QRHO)
|
|
|
H = 1.88*QRHO
|
|
|
H2 = 2*H
|
|
|
KAPN = IDNINT(7 + 34*QRHO)
|
|
|
NU = IDNINT(16 + 26*QRHO)
|
|
|
ENDIF
|
|
|
*
|
|
|
B = (H.GT.0.0)
|
|
|
*
|
|
|
IF (B) QLAMBDA = H2**KAPN
|
|
|
*
|
|
|
RX = 0.0
|
|
|
RY = 0.0
|
|
|
SX = 0.0
|
|
|
SY = 0.0
|
|
|
*
|
|
|
DO 11 N=NU, 0, -1
|
|
|
NP1 = N + 1
|
|
|
TX = YABS + H + NP1*RX
|
|
|
TY = XABS - NP1*RY
|
|
|
C = 0.5/(TX**2 + TY**2)
|
|
|
RX = C*TX
|
|
|
RY = C*TY
|
|
|
IF ((B).AND.(N.LE.KAPN)) THEN
|
|
|
TX = QLAMBDA + SX
|
|
|
SX = RX*TX - RY*SY
|
|
|
SY = RY*TX + RX*SY
|
|
|
QLAMBDA = QLAMBDA/H2
|
|
|
ENDIF
|
|
|
11 CONTINUE
|
|
|
*
|
|
|
IF (H.EQ.0.0) THEN
|
|
|
U = FACTOR*RX
|
|
|
V = FACTOR*RY
|
|
|
ELSE
|
|
|
U = FACTOR*SX
|
|
|
V = FACTOR*SY
|
|
|
END IF
|
|
|
*
|
|
|
IF (YABS.EQ.0.0) U = DEXP(-XABS**2)
|
|
|
*
|
|
|
END IF
|
|
|
*
|
|
|
*
|
|
|
C
|
|
|
C EVALUATION OF W(Z) IN THE OTHER QUADRANTS
|
|
|
C
|
|
|
*
|
|
|
IF (YI.LT.0.0) THEN
|
|
|
*
|
|
|
IF (A) THEN
|
|
|
U2 = 2*U2
|
|
|
V2 = 2*V2
|
|
|
ELSE
|
|
|
XQUAD = -XQUAD
|
|
|
*
|
|
|
C
|
|
|
C THE FOLLOWING IF-STATEMENT PROTECTS 2*EXP(-Z**2)
|
|
|
C AGAINST OVERFLOW
|
|
|
C
|
|
|
IF ((YQUAD.GT.RMAXGONI).OR.
|
|
|
* (XQUAD.GT.RMAXEXP)) GOTO 100
|
|
|
*
|
|
|
W1 = 2*DEXP(XQUAD)
|
|
|
U2 = W1*DCOS(YQUAD)
|
|
|
V2 = -W1*DSIN(YQUAD)
|
|
|
END IF
|
|
|
*
|
|
|
U = U2 - U
|
|
|
V = V2 - V
|
|
|
IF (XI.GT.0.0) V = -V
|
|
|
ELSE
|
|
|
IF (XI.LT.0.0) V = -V
|
|
|
END IF
|
|
|
*
|
|
|
RETURN
|
|
|
*
|
|
|
100 FLAG = .TRUE.
|
|
|
RETURN
|
|
|
*
|
|
|
END
|
|
|
|
