haversine distance from R matrices

haversineDist.cpp

#include <Rcpp.h>
#include <math.h>
#include <iostream>
using namespace Rcpp;

/* Convert degrees to radians */
inline double deg2rad(const double& deg){return deg*M_PI/180;};

/* Earth mean radius [km] */
const static double R_earth = 6371;

/* Calculates the geodesic distance between two points specified by 
 * radian latitude/longitude using the Haversine formula 
 * Ouputs distance between sites 1 and 2 as meters 
 */
inline double gcd_hf(const double& long1, const double& lat1, const double& long2, const double& lat2){
  double deltaLong = (long2 - long1);
  double deltaLat = (lat2 - lat1);
  double a = pow(sin(deltaLat/2),2) + cos(lat1) * cos(lat2) * pow(sin(deltaLong/2),2);
  double sqrtA = sqrtf(a);
  double c = 2 * asin(std::min(1.0,sqrtA));
  double d = (R_earth * c)*1000;
  return d;
};

/* Fxn to calculate matrix of distances between each two sites
 * INPUT: a matrix in which longs are in first column and lats in second column
 * OUTPUT: a distance matrix between all pairwise sites
 * Output distances are in meters
 */
// [[Rcpp::export]]
Rcpp::NumericMatrix calcDists(const Rcpp::NumericMatrix& longlats){
  int n = longlats.nrow();
  Rcpp::NumericMatrix zz = Rcpp::NumericMatrix(n,n);
  for(int i=0; i<n; i++){
    for(int j=0; j<n; j++){
      zz(i,j) = gcd_hf(deg2rad(longlats(i,1)), deg2rad(longlats(i,0)),deg2rad(longlats(j,1)), deg2rad(longlats(j,0)));
    }
  }
  return zz;
};