#include "geostat.h"
#include "chunk_array.h"
#include <time.h>

/*builds the sampled covariance function*/
/*dimensions are even*/
void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) {
    int i, j, k, maille, n2[3], symmetric;
    double di, dj, dk;

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

    //chunk_array_read(covar);

    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;
                chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));

                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));
                    double value;
                    chunk_array_get(covar, maille, &value);
                    chunk_array_save(covar, symmetric, value);
                }

                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);
                    chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
                }

                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));
                    double value;
                    chunk_array_get(covar, maille, &value);
                    chunk_array_save(covar, symmetric, value);               
                }

                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);
                    chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
                }

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

                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);
                    chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
                }

                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));
                    double value;
                    chunk_array_get(covar, maille, &value);
                    chunk_array_save(covar, symmetric, value);                      
                }
            }
        }
    }
}