integral Subroutine

SUBROUTINE integral(a, b, result, dk)
  !     Gaussian integration of exp(k.cosx) from a to b.
  REAL (dp), INTENT(IN) :: dk
  REAL, INTENT(IN)      :: a, b
  REAL, INTENT(OUT)     :: result
  !     Local variables
  REAL (dp)  :: xmid, range, x1, x2,                                    &
       x(3) = (/0.238619186083197_dp, 0.661209386466265_dp, 0.932469514203152_dp/), &
       w(3) = (/0.467913934572691_dp, 0.360761573048139_dp, 0.171324492379170_dp/)
  INTEGER    :: i
  xmid = (a + b)/2._dp
  range = (b - a)/2._dp
  result = 0._dp
  DO i = 1, 3
     x1 = xmid + x(i)*range
     x2 = xmid - x(i)*range
     result = result + w(i)*(EXP(dk*COS(x1)) + EXP(dk*COS(x2)))
  END DO
  result = result * range
  RETURN
END SUBROUTINE integral