Covariance Subroutine

SUBROUTINE Covariance &
!
(dist, theta, cov, sill, range, varmodel)

IMPLICIT NONE

!Arguments with intent(in)
REAL (KIND = float), INTENT(in) ::  sill  !! partial sill (total sill - nugget) from automatic fitting
REAL (KIND = float), INTENT(in) :: range  !! range from automatic fitting
REAL (KIND = float), DIMENSION(:), INTENT(IN) :: dist !!separation distance vector
REAL (KIND = float), DIMENSION(:), INTENT(IN) :: theta !!anisotropy angle ??
INTEGER (KIND = short), INTENT(in) :: varmodel  !!semivariogram model: 1 = spherical, 2 = exponential, 3 = gaussian, 4 = power

!Arguments with intent(out):
REAL (KIND = float), DIMENSION(:), INTENT(OUT) :: cov !!covariance vector


!Local variable declarations
REAL (KIND = float), DIMENSION(SIZE(dist)) :: h_prime
INTEGER (KIND = short) :: i

!------------------------------end of declarations-----------------------------

!h_prime = dist ! da correggere con anisotropia quando sarà implementata
    
!Check for anisotropy
IF (anisotropyAngle == 0.) THEN
	h_prime = dist
ELSE
	!Transform distance vector by applying anisotropy correction factor:
	!The minor axis of variation reaches sill before major axis of variation;
	!this is achieved by adding extra distance to the point separations 
	!to 'pull back' the range towards the origin. Separation distances 
	!are reprojected into a transformed coordinate system h'. 
	!First theta is rotated to the major axis of semivariance lies on the horizontal.
	!The anisotropy correction factor is calculated by multiplying the SIN of the 
	!angle between the major axis of variation and the point separation angle by the 
	!eccentricity of the anisotropy ellipsoid (major semiaxis/minor semiaxis). 
	!The transformed coordinate system h' is then the product of the separation 
	!and the correction factor.
    
	!h_prime = dist + dist * ABS( SIN( theta - rotation(i) ) * (theta3(i)/theta2(i)) )
    h_prime = dist + dist * ABS( SIN( theta - anisotropyAngle ) * anisotropyRatio )
END IF
        
!covariance = nugget + partialsill - semivariogram(h)
!since semivariogram(h) = nugget + partialsill * f(h,range) => covariance = partialsill * ( 1- f(h,range) )
SELECT CASE(varmodel)
    CASE (1) !Spherical
		WHERE(h_prime < range)
        cov = sill * (1. - ( 1.5 * (h_prime / range) - 0.5 * (h_prime / range)**3. ) )
		ELSEWHERE
        cov = 0.
        END WHERE
             
    CASE (2) !Exponential
        WHERE(h_prime < range)
        !cov = sill * (1. - ( 1. - EXP (- h_prime / range) ) )
        cov = sill * (1. - ( 1. - EXP (- 3 * h_prime / range) ) )
		ELSEWHERE
        cov = 0.
        END WHERE
             
    CASE (3) !Gaussian: MUST have a nugget effect to function without instability
        WHERE(h_prime < range)
        !cov = sill * (1. - ( 1. - EXP (- h_prime * h_prime / (range * range) ) ) ) 
        cov = sill * (1. - ( 1. - EXP (- 10 * h_prime * h_prime / (range * range) ) ) ) 
		ELSEWHERE
        cov = 0.
        END WHERE
             
    CASE (4) !Power
        WHERE(h_prime < range)
        cov = sill * (1. - (h_prime / range) ) 
		ELSEWHERE
        cov = 0.
        END WHERE

END SELECT 

RETURN

END SUBROUTINE Covariance