Trailing-Edge
-
PDP-10 Archives
-
decuslib10-01
-
43,50041/w.f4
There are no other files named w.f4 in the archive.
C THE PROBABILITY INTEGRAL FOR COMPLEX ARGUMENT. W(Z) 01
C (ERROR FUNCTION FOR COMPLEX ARGUMENT.) W(Z) 02
COMPLEX FUNCTION W(Z) W(Z) 03
C W(Z) = EXP(-Z*Z) * ERFC(-I*Z) W(Z) 04
C = EXP(-Z**2)*(1.+2.*J/ROOT PI*(INTEGRAL, 0 TO Z, (EXP(T**2).DT)))W(Z) 05
C = J/PI*(INTEGRAL, -INFINITY TO +INFINITY, (EXP(-T**2)/(Z-T).DT)) W(Z) 06
C W(Z) 07
C THIS FUNCTION HAS BEEN TABULATED TO SIX DECIMALS BY W(Z) 08
C FADDEYEVA AND TERENT@EV. W(Z) 09
C ( USSR ACADEMY OF SCIENCES. (1961) PERGAMON. ED. BY V.A. FOX.) W(Z) 10
C ALSO U.S. NATIONAL BUREAU OF STANDARDS. (1964) W(Z) 11
C @HANDBOOK OF MATHEMATICAL FUNCTIONS@ W(Z) 12
C W(Z) 13
C***********************************************************************W(Z) 14
COMPLEX Z,MZSQU W(Z) 15
DOUBLE PRECISION A,B,C,PART1,PART2,PART3,PART4,PART1S,PART3S W(Z) 16
DATA C, D / 1.1283791670955126D0, 4.35 / W(Z)
DATA CLOWER,C25,CUPPER/O145440672724,O232400000000,O221652211200/ W(Z) 18
ABS Z = CABS(Z) W(Z) 19
IF(ABS Z.LE.C LOWER) GO TO 10 W(Z) 20
IF(ABS Z.GE.C UPPER) GO TO 30 W(Z) 21
MZSQU = -Z*Z W(Z) 22
E = ( ( CABS(Z+D) + CABS(Z-D) )/2. )**2 - D*D W(Z) 23
IF(E.LT.1.) GO TO 20 W(Z) 24
C***********************************************************************W(Z) 25
C CONTINUED FRACTION METHOD. W(Z) 26
M = 4. + 40./E W(Z) 27
W = FLOAT(M) W(Z) 28
DO 5 N = 1,M W(Z) 29
RN = M - N + 1 W(Z) 30
5 W = RN*(RN-.5)/(MZSQU + 2.*RN + .5 - W) W(Z) 31
W = (0.,-.564189584)*Z/(.5+MZSQU-W) W(Z) 32
IF(REAL(MZSQU).LE.(-89.4159863)) RETURN W(Z) 35
IF(AIMAG(Z).EQ.0.) W = CMPLX(EXP(REAL(MZSQU)),AIMAG(W)) W(Z) 33
IF(AIMAG(Z).GE.0.) RETURN W(Z) 34
IF(ABS(AIMAG(MZSQU)).GE.C25) GO TO 201 W(Z) 36
IF(REAL(MZSQU).LT.(+87.336544 )) GO TO 8 W(Z) 37
WRITE(6,9200) Z W(Z) 38
MZSQU = CMPLX(87.336544 ,AIMAG(MZSQU)) W(Z) 39
8 W = W + 2.*CEXP(MZSQU) W(Z) 40
RETURN W(Z) 41
C RELATIVE INACCURACY AT COMPLEX ZEROS OF W IS IGNORED. W(Z) 42
C END OF CONTINUED FRACTION METHOD. W(Z) 43
C***********************************************************************W(Z) 44
10 W = CMPLX(1.,REAL(Z)*SNGL(C)) W(Z) 45
RETURN W(Z) 46
C***********************************************************************W(Z) 47
C DOUBLE-PRECISION CONVERGENT SERIES METHOD. W(Z) 48
20 PART 1 = 1. W(Z) 49
PART 2 = 0. W(Z) 50
PART 3 = 1. W(Z) 51
PART 4 = 0. W(Z) 52
M = 4. + ABS Z * (7.5 + 2.0 * ABS Z) W(Z) 53
A = REAL(MZSQU) W(Z) 54
B = AIMAG(MZSQU) W(Z) 55
DO 25 N = 1,M W(Z) 56
PART 1S = PART 1 W(Z) 57
PART 3S = PART 3 W(Z) 58
RN = M - N + 1 W(Z) 59
RNPLUS = RN + 0.5 W(Z) 60
PART 1 = 1. + (A*PART 1 - B*PART 2)/RN W(Z) 61
PART 2 = (B*PART 1S + A*PART 2)/RN W(Z) 62
PART 3 = 1. + (A*PART 3 - B*PART 4) /RNPLUS W(Z) 63
25 PART 4 = (A*PART 4 + B*PART 3S)/RNPLUS W(Z) 64
PART 3 = PART 3 * C W(Z) 65
PART 4 = PART 4 * C W(Z) 67
PART 1 = PART 1 - PART 3 * AIMAG(Z) - PART 4 * REAL(Z) W(Z) 68
PART 2 = PART 2 + PART 3 * REAL(Z) - PART 4 * AIMAG(Z) W(Z) 69
W = CMPLX(SNGL(PART 1), SNGL(PART 2)) W(Z) 70
RETURN W(Z) 71
C END OF CONVERGENT SERIES. W(Z) 72
C***********************************************************************W(Z) 73
30 W = (0.,.564189584)/Z W(Z) 74
IF(ABS(REAL(Z)).GT.(-AIMAG(Z))) RETURN W(Z) 75
201 TYPE 9201,Z W(Z) 76
RETURN W(Z) 77
9200 FORMAT(31H SINCE W(Z) OVERFLOWS FOR Z = (,1PE17.9,1H,,E17.9,2H),,/W(Z) 78
1,34H), AN APPROXIMATION HAS BEEN MADE.) W(Z) 79
9201 FORMAT(38H SINCE W(Z) IS INDETERMINATE FOR Z = (,1PE17.9,1H,,E17.9W(Z) 80
1,2H),,/,35H -W(-Z) HAS BEEN INSERTED FOR W(Z).) W(Z) 81
END W(Z) 82