function y = variogrammodel(type, r, model_para)
%   y = variogrammodel ( type, r, model_para)
%   |type| = model index for semivariogram/correlogram
%   type > 0: semi-variogram
%        < 0: correlogram  
%      r = vector lag distances
%  model_para = [ p1 p2 p3 p4 p5];
%      
%      p1: Nugt = nugget effect
%      p2: Sill = sill
%      p3: L = length scale for the main lobe
%      p4: b = length scale for hole effect     
%      p5: p = power for the expenential 
%
%      model type:
%          01 = spherical
%          02 = exponential
%          03 = gaussian
%          05 = linear

%      models including hole effects
%          20 = C * [ 1 - (sin b*r) / r ] + Nugt
%          21 = C * [ 1 - (exp(-r/L))     * cos(br) ] + Nugt
%          22 = C * [ 1 + (exp(-r/L))     * cos(br) ] + Nugt
%          23 = C * [ 1 - (exp(-(r)/L)    * cos(br) ] + Nugt
%          24 = C * [ 1 - (exp(-(r/L)^2)) * cos(br) ] + Nugt
%          25 = C * [ 1 -                   Jo (br) ] + Nugt
%          26 = C * [ 1 - (exp-(r/L))     * Jo (br) ] + Nugt
%          27 = C * [ 1 - (exp-(r/L)^2)   * Jo (br) ] + Nugt
%	        28 = C * [ 1 - (exp-(r/L)^2)   * (1 - br)] + Nugt
%          30 = C * [ 1 - (exp-(r/L)^p)   * Jo (br) ] + Nugt	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Author: D. Marcotte
% Version 2.1  97/aug/18
% Revised by Dezhang Chu,   10-29-98

Nugt=model_para(1);
Sill=model_para(2);
L=model_para(3);
b=model_para(4);
p=model_para(5);	
C=Sill-Nugt;
n = length(r);
rL = r ./ L;

switch abs(type) 
case 1
     indx1=find(rL < 1);
     indx2=find(rL >= 1);
     rL1=rL(indx1);
     rL2=rL(indx2);
     y1 = C * ( 1.5 .* rL1  -  0.5 .* rL1.^3)+Nugt;
     y2 = Sill * ones(size(rL2));
     y(indx1)=y1;
     y(indx2)=y2;
     y=reshape(y,size(r,1),size(r,2));
 case 2
     y = C * ( 1 - exp(-(r/L)))+Nugt;
 case 3
   y = C * ( 1 - exp(-(r/L).^2))+Nugt;
 case 5
   y = C .* r+Nugt;
 case 20
    y = C .* ( 1 - sin(b.*(r+eps))./(r+eps) )+Nugt;
 case 21
   y = C * ( 1 - exp(-r/L) .* cos(b*r) )+Nugt;
 case 22
   y = C * ( 1 + exp(-r/L) .* cos(b*r) )+Nugt;
 case 24
   y = C * (1 - exp(-(r/L).^2) .* cos(b*r) )+Nugt;
 case 25
   y = C * ( 1 - bessel(0, b*r) )+Nugt;
 case 26
   y = C * ( 1 - bessel(0, b*r) .* exp(-r/L) )+Nugt;
 case 27
   y = C * ( 1 - exp(-(r/L).^2) .* bessel(0, b*r) )+Nugt;
 case 28
    y = C * ( 1 - exp(-(r/L).^2) .* (1 - b*r.^2) )+Nugt;
 case 30							% generalized exponetial-Bessel 
    y= C* (1 - exp(- (r/L).^p).* bessel(0, b*r))+Nugt;		
end
