Web pdp-10.trailing-edge.com

Trailing-Edge - PDP-10 Archives - decuslib20-01 - decus/20-0011/w.for
There is 1 other file named w.for in the archive. Click here to see a list.
```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/"145440672724,"232400000000,"221652211200/ 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

```