#include <math.h>
#include "geostat.h"
#include "genlib.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);
}