core_dcmg.c

Go to the documentation of this file.

#include "../include/exageostatcore.h"

// This function converts decimal degrees to radians
static double deg2rad(double deg) {
  return (deg * PI / 180);
}
//  This function converts radians to decimal degrees
static double rad2deg(double rad) {
  return (rad * 180 / PI);
}
static double distanceEarth(double lat1d, double lon1d, double lat2d, double lon2d) {
  double lat1r, lon1r, lat2r, lon2r, u, v;
  lat1r = deg2rad(lat1d);
  lon1r = deg2rad(lon1d);
  lat2r = deg2rad(lat2d);
  lon2r = deg2rad(lon2d);
  u = sin((lat2r - lat1r)/2);
  v = sin((lon2r - lon1r)/2);
  return 2.0 * earthRadiusKm * asin(sqrt(u * u + cos(lat1r) * cos(lat2r) * v * v));
}

static double calculateDistance(double x1, double y1, double x2, double y2, int distance_metric) {

    if(distance_metric == 1)
        return distanceEarth(x1, y1, x2, y2);
        return  sqrt(pow((x2 - x1), 2) + pow((y2 - y1), 2));
}

/***************************************************************************/
void core_dcmg (double *A, int m, int n, int m0, int n0, location  *l1, location *l2, double *localtheta, int distance_metric) {
    
    int i, j;
    double l1x, l1y, l2x, l2y;
    double expr = 0.0;
    double con;
    double sigma_square = localtheta[0];// * localtheta[0]; 

    con = pow(2,(localtheta[2]-1)) * tgamma(localtheta[2]);
    con = 1.0/con;
    con = sigma_square * con;

    for (j = 0; j < n; j++) {
        l1x = l1->x[j+n0];      
        l1y = l1->y[j+n0];  
        for (i = 0; i < m; i++) {
            l2x = l2->x[i+m0];
            l2y = l2->y[i+m0];
                        expr = calculateDistance(l1x, l1y, l2x, l2y, distance_metric)/localtheta[1];
            if(expr == 0)
                                A[i + j * m] = sigma_square; /* + 1e-4*/
                        else
                                A[i + j * m] = con*pow(expr, localtheta[2])*gsl_sf_bessel_Knu(localtheta[2],expr); // Matern Function

        }
    }

}