!INTERFACE CovMat
!  SUBROUTINE MaternCov2d (nu, kappa, x, covariance, num_elem, theta_opt)
!  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
!    REAL ( kind = dp), INTENT( IN)  :: nu, kappa
!    REAL ( kind = dp), INTENT( IN), OPTIONAL  :: theta_opt
!    INTEGER, INTENT( IN)            :: num_elem
!    REAL ( kind = dp), INTENT( IN)  :: x(1:num_elem)
!    REAL ( kind = dp), INTENT(OUT)  :: covariance(1:num_elem)
!  END SUBROUTINE MaternCov2d  
!  SUBROUTINE HomogenExp ( phi, x, covariance, num_elem)
!  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
!    REAL ( kind = dp), INTENT( IN)  :: phi
!    INTEGER, INTENT( IN)            :: num_elem
!    REAL ( kind = dp), INTENT( IN)  :: x(1:num_elem)
!    REAL ( kind = dp), INTENT(OUT)  :: covariance(1:num_elem)
!  END SUBROUTINE
!END INTERFACE CovMat

SUBROUTINE HomogenExp ( phi, x, covariance, num_elem)
  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
  REAL ( kind = dp), INTENT( IN)  :: phi
  INTEGER, INTENT( IN)            :: num_elem
  REAL ( kind = dp), INTENT( IN)  :: x(1:num_elem)
  REAL ( kind = dp), INTENT(OUT)  :: covariance(1:num_elem)

  covariance = exp( -phi * x)

END SUBROUTINE

SUBROUTINE MaternCov (nu, kappa, x, covariance, num_elem, dim_opt, theta_opt)

  IMPLICIT INTEGER( I-K)
  
  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)

  INTEGER, PARAMETER          :: freq_n = 10000
  REAL (kind = dp), PARAMETER :: PI = 4*ATAN( 1.0_dp)

  REAL ( kind = dp), INTENT( IN)  :: nu, kappa
  REAL ( kind = dp), INTENT( IN), OPTIONAL  :: theta_opt, dim_opt
!  REAL ( kind = dp), INTENT( IN)  :: theta_opt, dim_opt
  INTEGER, INTENT( IN)            :: num_elem
  REAL ( kind = dp), INTENT( IN)  :: x(1:num_elem)
  REAL ( kind = dp), INTENT(OUT)  :: covariance(1:num_elem)

  REAL (kind = dp) :: theta, dimen, expon, vari, corr_eps, dw
  REAL (kind = dp) :: w(1:freq_n), spec(1:freq_n)
  LOGICAL          :: norm
  
  IF ( PRESENT( theta_opt)) THEN
    theta = theta_opt
  ELSE 
    theta = 0.0_dp
  ENDIF

  IF ( PRESENT( dim_opt)) THEN 
    dimen = dim_opt
  ELSE 
    dimen = 1
  ENDIF

  norm = .TRUE.

  IF ( theta .EQ. 0.0_dp ) THEN
    covariance = kappa * ABS(x)
    IF ( nu .GT. 0.0_dp) THEN
      vari = gamma(nu)/gamma(nu + dimen/2)/(4 * pi)**(dimen/2)/kappa**(2 * nu)
      IF (nu < 1) THEN
        expon = 2.0_dp * nu
      ELSE
        expon = 2.0_dp
      ENDIF
      corr_eps = 2.0_dp**(1-nu)/gamma(nu) * EPSILON( nu)**nu * BESSEL_K( EPSILON( nu), nu)
      WHERE ( covariance .GE. EPSILON ( nu))
        covariance = 2.0_dp**( 1-nu) / GAMMA( nu + dimen/2) / (4* PI)**(dimen/2) / &
          kappa**(2.0_dp * nu) * covariance**nu * BESSEL_K( covariance, nu)
      ELSEWHERE
        covariance = vari * (1.0_dp - (1 - corr_eps)/EPSILON( nu)**expon * covariance**expon)
      ENDWHERE
    ELSEIF ( nu .EQ. 0.0_dp ) THEN
      WHERE ( covariance .GE. EPSILON ( nu))
        covariance = 2.0_dp**( 1-nu) / GAMMA( nu + dimen/2) / (4* PI)**(dimen/2) / &
          kappa**(2.0_dp * nu) * covariance**nu * BESSEL_K( covariance, nu)
      ELSEWHERE
        covariance = 2/gamma( dimen/2.0_dp)/(4 * pi)**(dimen/2)* (-log(covariance/2) - 0.577215664901484)
      ENDWHERE
    ELSE
      WHERE ( covariance .GE. EPSILON ( nu))
        covariance = 2.0_dp**( 1-nu) / GAMMA( nu + dimen/2) / (4* PI)**(dimen/2) / &
          kappa**(2.0_dp * nu) * covariance**nu * BESSEL_K( covariance, nu)
      ELSEWHERE
        covariance = ( 2**(1-nu)/GAMMA( nu + dimen/2) / (4.0_dp* PI)**(dimen/2)/ &
          kappa**(2.0_dp * nu) * GAMMA( nu) * 2.0_dp**( nu - 1.0_dp)) * &
          (1.0_dp - (covariance / EPSILON( nu)) )+ 2**(1.0_dp - nu) / GAMMA(nu + dimen/2) /&
          (4.0_dp * PI)**(dimen/2)/kappa**(2.0_dp*nu) * EPSILON( nu)**nu * &
          BESSEL_K( EPSILON( nu), nu) * covariance/EPSILON( nu)
       ENDWHERE
    ENDIF
  ELSE
    covariance = ABS(x)
    IF (dimen .GT. 2.0_dp) THEN 
      write (*,*) "Dimension > 2 not implemented for oscillating models."
    ENDIF
    w(1) = 0
    DO I = 2, freq_n
      w( I) = W( I-1) + 1000/ MAXVAL( covariance) / freq_n
    ENDDO
    dw = w(2) - w(1)
    spec = 1/(2 * pi)**dimen/(kappa**4 + 2.0_dp * kappa**2 * cos(PI * theta) * w**2 + w**4)**((nu + dimen/2)/2)
    IF( dimen .EQ. 1.0_dp) THEN
      covariance = 0 + spec(1) * dw
      DO K = 2, freq_n
        covariance = covariance + 2 * cos( ABS(x) * w(k)) * spec(k) * dw
      ENDDO
    ELSE
      covariance = 0
      DO K = 2, freq_n
        covariance = covariance + w(k) * BESSEL_J0( ABS( x) * w(k)) * spec(k) * dw
      ENDDO
    ENDIF
    IF (norm) THEN
      covariance = covariance / covariance( 1)
    ENDIF
  ENDIF

END SUBROUTINE MaternCov

SUBROUTINE MaternCov2d_theta (nu, kappa, x, covariance, num_elem, theta)

  IMPLICIT INTEGER( I-K)
  
  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)

  INTEGER        :: freq_n
  REAL (kind = dp) :: freq_max
  REAL (kind = dp), PARAMETER :: PI = 4*ATAN( 1.0_dp)

  REAL ( kind = dp), INTENT( IN)  :: nu, kappa
  INTEGER, INTENT( IN)            :: num_elem
  REAL ( kind = dp), INTENT( IN)  :: x(1:num_elem,1:num_elem)
  REAL ( kind = dp), INTENT(OUT)  :: covariance(1:num_elem,1:num_elem)

  REAL (kind = dp) :: theta, expon, vari, corr_eps, dw, y(1:num_elem,1:num_elem)
  REAL (kind = dp), ALLOCATABLE :: w(:), spec(:)

  y = ABS( x)
  freq_max = 4/ MINVAL(y, MASK=y .GT. 0)
  freq_delta = 0.25_dp* 1/ MAXVAL( y)
  freq_n = CEILING( freq_max / freq_delta)

  ALLOCATE( w(1:freq_n))
  ALLOCATE( spec(1:freq_n))

  w(1) = 0
  DO I = 2, freq_n
    w( I) = (I-1) * freq_delta
  ENDDO
  dw = w(2) - w(1)
  spec = 1/(2* pi)**2/(kappa**4 + 2.0_dp * kappa**2 * cos(PI * theta) * w**2 + w**4)**((nu + 1)/2)
  covariance = 0
  DO K = 2, freq_n
    covariance = covariance + w(k) * BESSEL_J0( y * w(k)) * spec(k) * dw
  ENDDO
  covariance = covariance / covariance( 1,1)

END SUBROUTINE MaternCov2d_theta


SUBROUTINE MaternCov2d (nu, kappa, x, covariance, num_elem, theta_opt)

  IMPLICIT INTEGER( I-K)
  
  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)

  INTEGER, PARAMETER          :: freq_n = 10000, freq_max = 1000
  REAL (kind = dp), PARAMETER :: PI = 4*ATAN( 1.0_dp)

  REAL ( kind = dp), INTENT( IN)  :: nu, kappa
  REAL ( kind = dp), INTENT( IN), OPTIONAL  :: theta_opt
  INTEGER, INTENT( IN)            :: num_elem
  REAL ( kind = dp), INTENT( IN)  :: x(1:num_elem,1:num_elem)
  REAL ( kind = dp), INTENT(OUT)  :: covariance(1:num_elem,1:num_elem)

  REAL (kind = dp) :: theta, expon, vari, corr_eps, dw, y(1:num_elem,1:num_elem)
  REAL (kind = dp) :: w(1:freq_n), spec(1:freq_n)
  LOGICAL          :: norm
  
  IF ( PRESENT( theta_opt)) THEN
    theta = theta_opt
  ELSE 
    theta = 0.0_dp
  ENDIF

  norm = .TRUE.
  
  y = ABS( x)
  IF ( theta .EQ. 0.0_dp ) THEN
    covariance = kappa * y
    IF ( nu .GT. 0.0_dp) THEN
      vari = gamma(nu)/gamma(nu + 1)/(4 * pi)/kappa**(2 * nu)
      IF (nu < 1) THEN
        expon = 2.0_dp * nu
      ELSE
        expon = 2.0_dp
      ENDIF
      corr_eps = 2.0_dp**(1-nu)/gamma(nu) * EPSILON( nu)**nu * BESSEL_K( EPSILON( nu), nu)
      WHERE ( covariance .GE. EPSILON ( nu))
        covariance = 2.0_dp**( 1-nu) / GAMMA( nu + 1) / (4* PI) / &
          kappa**(2.0_dp * nu) * covariance**nu * BESSEL_K( covariance, nu)
      ELSEWHERE
        covariance = vari * (1.0_dp - (1 - corr_eps)/EPSILON( nu)**expon * covariance**expon)
      ENDWHERE
    ELSEIF ( nu .EQ. 0.0_dp ) THEN
      WHERE ( covariance .GE. EPSILON ( nu))
        covariance = 2.0_dp**( 1-nu) / GAMMA( nu + 1) / (4* PI) / &
          kappa**(2.0_dp * nu) * covariance**nu * BESSEL_K( covariance, nu)
      ELSEWHERE
        covariance = 2/gamma( 1.0_dp)/(4 * pi)* (-log(covariance/2) - 0.577215664901484)
      ENDWHERE
    ELSE
      WHERE ( covariance .GE. EPSILON ( nu))
        covariance = 2.0_dp**( 1-nu) / GAMMA( nu + 1) / (4* PI) / &
          kappa**(2.0_dp * nu) * covariance**nu * BESSEL_K( covariance, nu)
      ELSEWHERE
        covariance = ( 2**(1-nu)/GAMMA( nu + 1) / (4.0_dp* PI)/ &
          kappa**(2.0_dp * nu) * GAMMA( nu) * 2.0_dp**( nu - 1.0_dp)) * &
          (1.0_dp - (covariance / EPSILON( nu)) )+ 2**(1.0_dp - nu) / GAMMA(nu + dimen/2) /&
          (4.0_dp * PI) / kappa**(2.0_dp*nu) * EPSILON( nu)**nu * &
          BESSEL_K( EPSILON( nu), nu) * covariance/EPSILON( nu)
       ENDWHERE
    ENDIF
  ELSE
    w(1) = 0
    DO I = 2, freq_n
      w( I) = (I-1)* freq_max/ MAXVAL( y) / freq_n
    ENDDO
    dw = w(2) - w(1)
    spec = 1/(2* pi)**2/(kappa**4 + 2.0_dp * kappa**2 * cos(PI * theta) * w**2 + w**4)**((nu + 1)/2)
    covariance = 0
    DO K = 2, freq_n
      covariance = covariance + w(k) * BESSEL_J0( y * w(k)) * spec(k) * dw
    ENDDO
  ENDIF
  IF (norm) THEN
    covariance = covariance / covariance( 1,1)
  ENDIF

END SUBROUTINE MaternCov2d

