Trailing-Edge
-
PDP-10 Archives
-
decus_20tap1_198111
-
decus/20-0011/w.des
There are 2 other files named w.des in the archive. Click here to see a list.
UNIVERSITY OF WESTERN AUSTRALIA
COMPUTING CENTRE
PROGRAM INFORMATION SHEET NO. 62
W
PROGRAM NO: UWA6-07.0.004
AUTHOR: I.D.PUGSLEY, PHYSIOLOGY
SOURCE LANGUAGE: FORTRAN IV
DESCRIPTION: THIS FORTRAN IV FUNCTION WILL CALCULATE THE PROBABILITY
INTEGRAL OR ERROR FUNCTION WHICH ARISES IN PROBLEMS OF
DIFFUSION, HEAT FLOW AND DISTRIBUTED ELECTRICAL NETWORKS.
THE ERROR FUNCTION IS COMPUTED FOR A COMPLEX ARGUMENT;
COMPLEX NUMBERS ARISE WHERE AN ELEMENT OF THE SYSTEM
CONTAINS A TIME-DEPENDENT FUNCTION.
SPECIAL CASES OF THIS FUNCTION MAY BE USED TO COMPUTE
ERROR FUNCTION AND ITS COMPLEMENT FOR REAL ARGUMENT,
INTEGRAL ERROR FUNCTIONS, PROBABILITY FUNCTIONS,
DAWSON'S INTEGRAL, FRESNEL INTEGRALS, CONFLUENT HYPER-
GEOMETRIC FUNCTION, PARABOLIC CYLINDER FUNCTIONS AND
SPHERICAL BESSEL FUNCTIONS.
FOR THESE AND FURTHER APPLICATIONS AND EQUATIONS SEE
'HANDBOOK OF MATHEMATICAL FUNCTIONS'
ED. ABRAMOWITZ AND STEGUN, N.B.S., A.M.SERIES 55, 1964.
CALLING SEQUENCE: W(Z) WHERE W IS A COMPLEX FUNCTION
Z IS A COMPLEX ARGUMENT
IF OVERFLOW OCCURS NEAR SINGULARITIES, A MESSAGE IS
TYPED AND COMPUTATION CONTINUES.
ERROR FUNCTION IS
Z
E**(-Z**2)(1+(2*I/SQRT(PI))*INT(E**(T**2)DT) = W(Z)
Z 0
WHERE INT DENOTES INTEGRAL FROM 0 TO Z
0
I IS SQRT(-1)
E IS 2.71828...
PI IS 3.14159...
PAGE 2
SOURCE LISTING:
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
PAGE 3
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
************************************************************************