view release on metacpan or  search on metacpan

XS.xs  view on Meta::CPAN

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923132169163975144
#endif

const double DEG_RADS = M_PI / 180.;

/* From Geo::Distance */
const double KILOMETER_RHO = 6371.64;

/* WGS 84 Ellipsoid */
const double A = 6378137.;
const double B = 6356752.314245;
const double F = 1. / 298.257223563;

static void
my_croak (char* pat, ...) {
    va_list args;
    SV *error_sv;

    dTHX;
    dSP;

    error_sv = newSV(0);

    va_start(args, pat);
    sv_vsetpvf(error_sv, pat, &args);
    va_end(args);

    ENTER;
    SAVETMPS;
    PUSHMARK(SP);
    XPUSHs(sv_2mortal(error_sv));
    PUTBACK;
    call_pv("Carp::croak", G_VOID | G_DISCARD);
    FREETMPS;
    LEAVE;
}

double
haversine (double lat1, double lon1, double lat2, double lon2) {
    lat1 *= DEG_RADS; lon1 *= DEG_RADS;
    lat2 *= DEG_RADS; lon2 *= DEG_RADS;
    double a = sin(0.5 * (lat2 - lat1));
    double b = sin(0.5 * (lon2 - lon1));
    double c = a * a + cos(lat1) * cos(lat2) * b * b;
    double d = 2. * atan2(sqrt(c), sqrt(fabs(1. - c)));
    return d;
}

double
cosines (double lat1, double lon1, double lat2, double lon2) {
    lat1 *= DEG_RADS; lon1 *= DEG_RADS;
    lat2 *= DEG_RADS; lon2 *= DEG_RADS;
    double a = sin(lat1) * sin(lat2);
    double b = cos(lat1) * cos(lat2) * cos(lon2 - lon1);
    double d = acos(a + b);
    /* Antipodal coordinates result in NaN */
    if (isnan(d))
        return haversine(lat1, lon1, lat2, lon2);
    return d;
}

double
polar (double lat1, double lon1, double lat2, double lon2) {
    double a = M_PI_2 - lat1 * DEG_RADS;
    double b = M_PI_2 - lat2 * DEG_RADS;
    double dlon = (lon2 - lon1) * DEG_RADS;
    double d = sqrt(a * a + b * b - 2. * a * b * cos(dlon));
    return d;
}

double
great_circle (double lat1, double lon1, double lat2 , double lon2) {
    lat1 *= DEG_RADS; lon1 *= DEG_RADS;
    lat2 *= DEG_RADS; lon2 *= DEG_RADS;
    double a = sin(0.5 * (lat2 - lat1));
    double b = sin(0.5 * (lon2 - lon1));
    double c = a * a + cos(lat1) * cos(lat2) * b * b;
    double d = 2. * asin(sqrt(c));
    return d;
}

double
vincenty (double lat1, double lon1, double lat2 , double lon2) {
    double dlon = (lon2 - lon1) * DEG_RADS;
    double u1 = atan((1. - F) * tan(lat1 * DEG_RADS));
    double u2 = atan((1. - F) * tan(lat2 * DEG_RADS));
    double sin_u1 = sin(u1), cos_u1 = cos(u1);
    double sin_u2 = sin(u2), cos_u2 = cos(u2);

    double lambda = dlon, lambda_p = 2. * M_PI;
    int iter_limit = 100;

    double sin_sigma, cos_sigma;
    double sigma;
    double cos_sq_alpha, cos_sigma_m;
    double u_sq, a, b, delta_sigma, d;

    while (fabs(lambda - lambda_p) > 1e-12 && iter_limit-- > 0) {
        double alpha, c;
        double sin_lambda = sin(lambda);
        double cos_lambda = cos(lambda);
        sin_sigma = sqrt((cos_u2 * sin_lambda) * (cos_u2 * sin_lambda) +
                         (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda) *
                         (cos_u1 * sin_u2-sin_u1 * cos_u2 * cos_lambda));
        if (sin_sigma == 0.) {
            return 0.;
        }
        cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
        sigma = atan2(sin_sigma, cos_sigma);
        alpha = asin(cos_u1 * cos_u2 * sin_lambda / sin_sigma);
        cos_sq_alpha = cos(alpha) * cos(alpha);
        cos_sigma_m = cos_sigma - 2. * sin_u1 * sin_u2 / cos_sq_alpha;
        if (isnan(cos_sigma_m)) {
            cos_sigma_m = 0.;
        }
        c = 0.0625 * F * cos_sq_alpha *