source: issm/branches/trunk-jpl-damage/src/c/classes/kriging/GaussianVariogram.cpp@ 12878

/*!\file GaussianVariogram.c
 * \brief: implementation of the GaussianVariogram object
 */

#ifdef HAVE_CONFIG_H
        #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <stdio.h>
#include <string.h>
#include "../objects.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../shared/shared.h"
#include "../../include/include.h"

/*GaussianVariogram constructors and destructor*/
/*FUNCTION GaussianVariogram::GaussianVariogram(){{{*/
GaussianVariogram::GaussianVariogram(){
        this->nugget = 0.2;
        this->sill   = 1;
        this->range  = SQRT3;
        return;
}
/*}}}*/
/*FUNCTION GaussianVariogram::GaussianVariogram(Options* options){{{*/
GaussianVariogram::GaussianVariogram(Options* options){

        /*Defaults*/
        this->nugget = 0.2;
        this->sill   = 1;
        this->range  = SQRT3;

        /*Overwrite from options*/
        if(options->GetOption("nugget")) options->Get(&this->nugget,"nugget");
        if(options->GetOption("sill"))   options->Get(&this->sill,"sill");
        if(options->GetOption("range"))  options->Get(&this->range,"range");

        /*Checks*/
        if(nugget==sill) _error2_("nugget and sill cannot be equal (constant semivariogram not allowed)");
}
/*}}}*/
/*FUNCTION GaussianVariogram::~GaussianVariogram(){{{*/
GaussianVariogram::~GaussianVariogram(){
        return;
}
/*}}}*/

/*Object virtual functions definitions:*/
/*FUNCTION GaussianVariogram::Echo {{{*/
void GaussianVariogram::Echo(void){
        _printLine_("GaussianVariogram");
        _printLine_("   nugget: " << this->nugget);
        _printLine_("   sill  : " << this->sill);
        _printLine_("   range : " << this->range);
}
/*}}}*/

/*Variogram function*/
/*FUNCTION GaussianVariogram::Covariance{{{*/
double GaussianVariogram::Covariance(double deltax,double deltay){
        /*The covariance can be deduced from the variogram from the following
         * relationship:
         *    2 gamma = C(x,x) + C(y,y) -2 C(x,y)
         * so
         *    C(h) = sill - gamma                                            */
        double h2,a,cova;

        /*Calculate length square*/
        h2=pow(deltax,2.)+pow(deltay,2.);

        /*If h is too small, return sill*/
        if(h2<0.0000001) return sill;

        /*compute covariance*/
        a     = 1./3.;
        cova = (sill-nugget)*exp(-h2/(a*range*range));

        return cova;
}
/*}}}*/
/*FUNCTION GaussianVariogram::SemiVariogram{{{*/
double GaussianVariogram::SemiVariogram(double deltax,double deltay){
        /*http://en.wikipedia.org/wiki/Variogram*/
        double h2,a,gamma;

        /*Calculate length square*/
        h2=pow(deltax,2.)+pow(deltay,2.);

        /*return semi-variogram*/
        a     = 1./3.;
        gamma = (sill-nugget)*(1.-exp(-h2/(a*range*range))) + nugget;

        //if(h2>1000*1000) _printLine_("gamma = " << gamma << " h= " << sqrt(h2));
        _printLine_("h = " << sqrt(h2) << " gamma = " << gamma);
        return gamma;
}
/*}}}*/