Google
 

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


************************************************************************