random_von_Mises Function

FUNCTION random_von_Mises(k, first) RESULT(fn_val)
  !     Algorithm VMD from:
  !     Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a
  !     comparison of random numbers', J. of Appl. Statist., 17, 165-168.
  !     Fortran 90 code by Alan Miller
  !     CSIRO Division of Mathematical & Information Sciences
  !     Arguments:
  !     k (real)        parameter of the von Mises distribution.
  !     first (logical) set to .TRUE. the first time that the function
  !                     is called, or the first time with a new value
  !                     for k.   When first = .TRUE., the function sets
  !                     up starting values and may be very much slower.
  REAL, INTENT(IN)     :: k
  LOGICAL, INTENT(IN)  :: first
  REAL                 :: fn_val
  !     Local variables
  INTEGER          :: j, n
  INTEGER, SAVE    :: nk
  REAL, PARAMETER  :: pi = 3.14159265
  REAL, SAVE       :: p(20), theta(0:20)
  REAL             :: sump, r, th, lambda, rlast
  REAL (dp)        :: dk
  IF (first) THEN                        ! Initialization, if necessary
     IF (k < zero) THEN
        WRITE(*, *) '** Error: argument k for random_von_Mises = ', k
        RETURN
     END IF
     nk = k + k + one
     IF (nk > 20) THEN
        WRITE(*, *) '** Error: argument k for random_von_Mises = ', k
        RETURN
     END IF
     dk = k
     theta(0) = zero
     IF (k > half) THEN
        !     Set up array p of probabilities.
        sump = zero
        DO j = 1, nk
           IF (j < nk) THEN
              theta(j) = ACOS(one - j/k)
           ELSE
              theta(nk) = pi
           END IF
           !     Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j)
           CALL integral(theta(j-1), theta(j), p(j), dk)
           sump = sump + p(j)
        END DO
        p(1:nk) = p(1:nk) / sump
     ELSE
        p(1) = one
        theta(1) = pi
     END IF                         ! if k > 0.5
  END IF                           ! if first
  CALL RANDOM_NUMBER(r)
  DO j = 1, nk
     r = r - p(j)
     IF (r < zero) EXIT
  END DO
  r = -r/p(j)
  DO
     th = theta(j-1) + r*(theta(j) - theta(j-1))
     lambda = k - j + one - k*COS(th)
     n = 1
     rlast = lambda
     DO
        CALL RANDOM_NUMBER(r)
        IF (r > rlast) EXIT
        n = n + 1
        rlast = r
     END DO
     IF (n .NE. 2*(n/2)) EXIT         ! is n even?
     CALL RANDOM_NUMBER(r)
  END DO
  fn_val = SIGN(th, (r - rlast)/(one - rlast) - half)
  RETURN
END FUNCTION random_von_Mises