function model = variogrammodels(x, p, type)
% model = variogramModel(x, p, type),
% Calculates a variogram model from the paramters p, for distances x
% The following variogram model types are supported:
% - 'linear' 
% - 'linear_with_sill'
% - 'spherical'
% - 'exponential' 
% - 'gausian'
% - 'polynomial' experimental model, not an idealized description of
%    semivariance, but a close fit to the experimental variogram, 
%    often produces nice results, but can cause strange extrapolation 
%    effects when used for kriging.
% 
% The parameter p should be a 3x1 vector, except for polynomial, where it can contain more elements.
% 
% Yannick Kremer 2007

% C0 : nugget
% C : sill break
% a : slope 0..range

switch lower(type)
    case 'linear' 
        C0 =  p(1); 
        if C0<0, C0=0;, end;
        
        a = p(2);
        model = C0 + a*x;
    case 'linear_with_sill' 
        C0 = p(1); a =  p(2); C = p(3);
        if C0<0, C0=0;, end;
        
        model = zeros(size(x));        
        model( x<C  ) = C0 + a*x( x<C );
        model( x>=C ) = C0 + a*C; 
        model( x==0 ) = 0;
    case 'spherical'
        C0 = p(1); a = p(2); C = p(3);
        if C0<0, C0=0;, end;
        
        model = zeros(size(x));
        model( x<=a ) = C0 + C*( 1.5*( x(x<=a)/a) - 0.5*( x(x<=a)/a ).^3);
        model( x>a ) = C0 + C; 
        model( x==0 ) = 0;
    case 'exponential'
        C0 = p(1); a = p(2); C = p(3);  
        if C0<0, C0=0;, end;
        
        model = C0 + C*(1 - exp( -x/a ));
        model( x==0) = 0;
    case 'gaussian'
        C0 = p(1); a = p(2); C = p(3);
        if C0<0, C0=0;, end;
        
        model = C0 + C*(1 - exp(-(x/a).^2));
        model( x==0) = 0;
    case 'polynomial'
        model = polyval(p, x);
        model( x==0 ) = 0;
        model( model<0) = 0;
    otherwise
        error('Error -  invalid model specified, please choose linear, linear_with_sill, spherical, exponential, gaussian or polynomial');
end   