FUNCTION VMISES(T, VK)                                            VMI   10
C VMISES RETURNS THE LEFT TAIL AREA OF THE VON MISES DISTRIBUTION,
C EQUAL TO THE INCOMPLETE MODIFIED BESSEL FUNCTION,  OF THE FIRST
C KIND AND ZERO-TH ORDER.
C     T  = ANGLE IN RADIANS, TREATED AS DEVIATION FROM ZERO (MEAN)
C            REDUCED MODULO 2PI TO THE RANGE (-PI,+PI).
C     VK = CONCENTRATION PARAMETER, KAPPA.    NEGATIVE VALUES ARE
C            TREATED AS ZERO.
C NEEDS AS AUXILIARY ROUTINE EITHER
C FUNCTION GAUSS(X),  RETURNING THE NORMAL DISTRIBUTION TAIL AREA
C                         TO THE LEFT OF THE BOUNDING ORDINATE, X,
C OR  FUNCTION ERF(X),  RETURNING THE ERROR FUNCTION AREA BETWEEN
C      ZERO AND THE BOUNDING ORDINATE, X, WHERE ERF(-X) = -ERF(X).
C      IF ERF IS USED THEN EACH OF THE LAST 3 COMMENTS (COLUMNS 7
C      TO 50 ONLY) MUST REPLACE THE STATEMEMT PRECEDING IT.
      DATA PI /3.1415926535898/, TPI /6.2831853071760/
C CONSTANTS APPROPRIATE FOR 8 DECIMAL DIGIT ACCURACY.
      DATA A1, A2, A3, A4, CK, C1 /12.0,0.8,8.0,1.0,10.5,56.0/
      Z = VK
C CONVERT ANGLE T MODULO 2PI TO RANGE (-PI,+PI).
      U = AMOD(T+PI,TPI)
      IF (U.LT.0.0) U = U + TPI
      Y = U - PI
      IF (Z.GT.CK) GO TO 30
      V = 0.0
      IF (Z.LE.0.0) GO TO 20
C FOR SMALL VK SUM IP TERMS BY BACKWARDS RECURSION.
      IP = Z*A2 - A3/(Z+A4) + A1
      P = FLOAT(IP)
      S = SIN(Y)
      C = COS(Y)
      Y = P*Y
      SN = SIN(Y)
      CN = COS(Y)
      R = 0.0
      Z = 2.0/Z
      DO 10 N=2,IP
        P = P - 1.0
        Y = SN
        SN = SN*C - CN*S
        CN = CN*C + Y*S
        R = 1.0/(P*Z+R)
        V = (SN/P+V)*R
   10 CONTINUE
   20 VMISES = (U*0.5+V)/PI
      GO TO 40
C FOR LARGE VK COMPUTE THE NORMAL APPROXIMATION AND LEFT TAIL.
   30 C = 24.0*Z
      V = C - C1
      R = SQRT((54.0/(347.0/V+26.0-C)-6.0+C)/6.0)
C     R=SQRT((54.0/(347.0/V+26.0-C)-6.0+C)/12.0)           IF ERF
      Z = SIN(Y*0.5)*R
      S = Z*Z
C     S=Z*Z*2.0                                            IF ERF
      V = V - S + 3.0
      Y = (C-S-S-16.0)/3.0
      Y = ((S+1.75)*S+83.5)/V - Y
      VMISES = GAUSS(Z-S/(Y*Y)*Z)
C     VMISES=ERF(Z-S/(Y*Y)*Z)*0.5+0.5                      IF ERF
   40 IF (VMISES.LT.0.0) VMISES = 0.0
      IF (VMISES.GT.1.0) VMISES = 1.0
      RETURN
      END

.