#include "genlib.h"
#include "geostat.h"
#include <math.h>
#include <time.h>

/*selection of model covariance*/
double cov_value(struct vario_mod variogram, double di, double dj, double dk) {
    double hx, hy, hz, h;
    double cov;
    int k;

    cov = 0.;

    for (k = 0; k < variogram.Nvario; k++) {

        hx = di * variogram.ap[9 * k] + dj * variogram.ap[9 * k + 1] + dk * variogram.ap[9 * k + 2];
        hy = di * variogram.ap[9 * k + 3] + dj * variogram.ap[9 * k + 4] + dk * variogram.ap[9 * k + 5];
        hz = di * variogram.ap[9 * k + 6] + dj * variogram.ap[9 * k + 7] + dk * variogram.ap[9 * k + 8];
        h = sqrt(hx * hx + hy * hy + hz * hz);

        switch (variogram.vario[k]) {
        case 1:
            cov += variogram.var[k] * exponential(h);
            break;
        case 2:
            cov += variogram.var[k] * gaussian(h);
            break;
        case 3:
            cov += variogram.var[k] * spherical(h);
            break;
        case 4:
            cov += variogram.var[k] * cardsin(h);
            break;
        case 5:
            cov += variogram.var[k] * stable(h, variogram.alpha[k]);
            break;
        case 6:
            cov += variogram.var[k] * gammf(h, variogram.alpha[k]);
            break;
        case 7:
            cov += variogram.var[k] * cubic(h);
            break;
        case 8:
            cov += variogram.var[k] * nugget(h);
            break;
        case 9:
            cov += variogram.var[k] * power(h, variogram.alpha[k]);
            break;
        }
    }
    return cov;
}