#include "geostat.h"

/*builds the sampled covariance function*/
/*dimensions are even*/
void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3])

{
    int i, j, k, maille, n2[3], symmetric;
    double di, dj, dk;

    for (i = 0; i < 3; i++)
        n2[i] = n[i] / 2;

    for (i = 0; i <= n2[0]; i++) {
        for (j = 0; j <= n2[1]; j++) {
            for (k = 0; k <= n2[2]; k++) {

                /*area 1*/
                maille = 1 + i + n[0] * (j + n[1] * k);
                di = (double)i * mesh.DX;
                dj = (double)j * mesh.DY;
                dk = (double)k * mesh.DZ;
                covar[maille] = (double)cov_value(variogram, di, dj, dk);

                if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) {
                    /*area 2*/
                    symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k));
                    covar[symmetric] = covar[maille];
                }

                if (i > 0 && i < n2[0]) {
                    /*area 4*/
                    di = -(double)i * mesh.DX;
                    dj = (double)j * mesh.DY;
                    dk = (double)k * mesh.DZ;
                    maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k);
                    covar[maille] = (double)cov_value(variogram, di, dj, dk);
                }

                if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) {
                    /*area 8*/
                    symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k));
                    covar[symmetric] = covar[maille];
                }

                if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) {
                    /*area 5*/
                    di = -(double)i * mesh.DX;
                    dj = -(double)j * mesh.DY;
                    dk = (double)k * mesh.DZ;
                    maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k);
                    covar[maille] = (double)cov_value(variogram, di, dj, dk);
                }

                if (k > 0 && k < n2[2]) {
                    /*area 6*/
                    symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k));
                    covar[symmetric] = covar[maille];
                }

                if (j > 0 && j < n2[1]) {
                    /*area 3*/
                    di = (double)i * mesh.DX;
                    dj = -(double)j * mesh.DY;
                    dk = (double)k * mesh.DZ;
                    maille = 1 + i + n[0] * (n[1] - j + n[1] * k);
                    covar[maille] = (double)cov_value(variogram, di, dj, dk);
                }

                if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) {
                    /*area 7*/
                    symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k));
                    covar[symmetric] = covar[maille];
                }
            }
        }
    }

    return;
}