view release on metacpan or  search on metacpan

XS.xs  view on Meta::CPAN

    M = radius * ( ( 1 - 0.25 * eccentricity_2 - (3./64.) * eccentricity_4 - (5./256.) * eccentricity_6 ) * latitude
                   - ( (3./8.) * eccentricity_2 + (3./32.) * eccentricity_4 + (45./1024.) * eccentricity_6) * sin_2_latitude
                   + ( ( 15./256.) * eccentricity_4 + (45./1024.) * eccentricity_6) * sin_4_latitude
                   - (35./3072.) * eccentricity_6 * sin_6_latitude);

    A_2 = A * A;
    A_3 = A_2 * A;
    A_4 = A_3 * A;
    A_5 = A_4 * A;
    A_6 = A_5 * A;

    tan_latitude_4 = tan_latitude_2 * tan_latitude_2;
        
    *easting =   500000.0 + k0 * N * (A + (1 - tan_latitude_2 + C) * A_3 * (1./6.) + (5 - 18 * tan_latitude_2 + tan_latitude_4 + 72 * C - 58 * eccentricity_prime_2) * A_5 * (1./120.));

    *northing= k0 * ( M + N * tan_latitude * ( A * A / 2 + (5 - tan_latitude_2 + 9 * C + 4 * C * C) * A_4 * (1./ 24.)
                                             + (61 - 58 * tan_latitude_2 + tan_latitude_4 + 600 * C - 330 * eccentricity_prime_2) * A_6 * (1./720.)));

    if ((*zone_letter) < 'N')
        *northing += 10000000.0;
}
    

MODULE = Geo::Coordinates::UTM::XS		PACKAGE = Geo::Coordinates::UTM::XS		
PROTOTYPES: ENABLE

BOOT:
memset(ellipsoids, 0, sizeof(ellipsoids));
ellipsoid_hv = GvHV(gv_fetchpv("Geo::Coordinates::UTM::XS::_ellipsoid", TRUE, SVt_PVHV));

void
_set_ellipsoid_info(index, radius, eccentricity_2)
    int index
    double radius
    double eccentricity_2
PROTOTYPE: @
CODE:
    {
        if (index < 1 || index >= MAX_ELLIPSOIDS || ellipsoids[index].defined) {
            Perl_croak(aTHX_ "invalid ellipsoid index %i", index);
        }
        ellipsoids[index].radius = radius;
        ellipsoids[index].inv_radius = 1.0/radius;
        ellipsoids[index].eccentricity_2 = eccentricity_2;
        ellipsoids[index].eccentricity_4 = eccentricity_2 * eccentricity_2;
        ellipsoids[index].eccentricity_6 = eccentricity_2 * eccentricity_2 * eccentricity_2;
        ellipsoids[index].eccentricity_prime_2 = eccentricity_2/(1-eccentricity_2);
        ellipsoids[index].defined = 1;
    }

void
_latlon_to_utm(ename, latitude_deg, longitude_deg)
    SV *ename
    double latitude_deg
    double longitude_deg
PROTOTYPE: $$$
PREINIT:
    int zone = 0;
    char zone_letter = 0;
    double easting, northing;
PPCODE:
    _latlon_to_utm(aTHX_ ename, latitude_deg, longitude_deg,
                   &zone, &zone_letter, &easting, &northing);
    XPUSHs(sv_2mortal(newSVpvf("%d%c", zone, zone_letter)));
    XPUSHs(sv_2mortal(newSVnv(easting)));
    XPUSHs(sv_2mortal(newSVnv(northing)));
    XSRETURN(3);

void
_latlon_to_utm_force_zone(ename, zone, latitude_deg, longitude_deg)
    SV *ename
    SV *zone
    double latitude_deg
    double longitude_deg
PROTOTYPE: $$$$
PREINIT:
    int zone_number;
    char zone_letter;
    double easting, northing;
PPCODE:
    zone_letter = 0;
    _zonesv_to_number_letter(aTHX_ zone, &zone_number, &zone_letter);
    if (zone_number < 0 || zone_number > 60)
        Perl_croak(aTHX_ "Zone value (%d) invalid.", zone_number);
    _latlon_to_utm(aTHX_ ename, latitude_deg, longitude_deg,
                   &zone_number, &zone_letter, &easting, &northing);
    XPUSHs(sv_2mortal(newSVpvf("%d%c", zone_number, zone_letter)));
    XPUSHs(sv_2mortal(newSVnv(easting)));
    XPUSHs(sv_2mortal(newSVnv(northing)));
    XSRETURN(3);



void
_utm_to_latlon(ename, zone, easting, northing)
    SV *ename
    SV *zone
    double easting
    double northing
PROTOTYPE: $$$$
PPCODE:
    {
        int index;
        double radius, inv_radius, eccentricity_2, eccentricity_4, eccentricity_6, eccentricity_prime_2;
        double x, y;
        double longitude_origin_deg;
        double M, mu, e1, e1_2, e1_3, e1_4, phi1, N1, N2, N3, tan_phi1_2, C1, C1_2, R1, D, D_2, D_3, D_4, D_5, D_6;
        double sin_phi1, cos_phi1, inv_cos_phi1, tan_phi1, sin_phi1_2, tan_phi1_4;
        double latitude_deg, longitude_deg;
        int zone_number;
        char zone_letter;
        int i, zone_ok;

        index = ellipsoid_index(aTHX_ ename);
        
        if (index < 1 || index >= MAX_ELLIPSOIDS || !ellipsoids[index].defined) {
            Perl_croak(aTHX_ "invalid ellipsoid index %i", index);
        }

        radius = ellipsoids[index].radius;
        inv_radius = ellipsoids[index].inv_radius;

        eccentricity_2 = ellipsoids[index].eccentricity_2;
        eccentricity_4 = ellipsoids[index].eccentricity_4;
        eccentricity_6 = ellipsoids[index].eccentricity_6;
        eccentricity_prime_2 = ellipsoids[index].eccentricity_prime_2;

        zone_letter = 'N';
        _zonesv_to_number_letter(aTHX_ zone, &zone_number, &zone_letter);

        x = easting - 500000.0;
        y = northing;
        if (zone_letter < 'N')
            y -= 10000000.0;
        
        M = inv_k0 * y;
        mu = M / (radius * (1 - eccentricity_2 * (1./4.) - eccentricity_4 * (3./64.) - eccentricity_6 * (5./256.)));
        e1 = -1 + 2 * (1 - sqrt(1 - eccentricity_2)) / eccentricity_2;

        e1_2 = e1 * e1;
        e1_3 = e1_2 * e1;
        e1_4 = e1_3 * e1;
        
        phi1 = mu + (1.5 * e1 - e1_3 * (27./32.)) * sin(2 * mu) + (e1_2 * (21./16.) - e1_4 * (55./32.)) * sin(4 * mu) + (e1_3 * (151./96.)) * sin(6 * mu);

        sin_phi1 = sin(phi1);
        cos_phi1 = cos(phi1);
        inv_cos_phi1 = 1.0 / cos_phi1;
        tan_phi1 = sin_phi1 * inv_cos_phi1;

        sin_phi1_2 = sin_phi1 * sin_phi1;

        N3 = sqrt(1 - eccentricity_2 * sin_phi1_2);
        N2 = 1.0 / N3;
        N1 = radius * N2;
        
        tan_phi1_2 = tan_phi1 * tan_phi1;
        tan_phi1_4 = tan_phi1_2 * tan_phi1_2;
        
        C1 = eccentricity_2 * cos_phi1 * cos_phi1;
        C1_2 = C1 * C1;