view release on metacpan or  search on metacpan

XS.xs  view on Meta::CPAN

    long double r_major    = 6378137;           // Equatorial Radius, WGS84
    long double r_minor    = 6356752.314245179; // defined as constant
    long double f          = 1/298.257223563;   // 1/f = ( $r_major - $r_minor ) / $r_major

    long double alpha1     = DEG2RAD( bearing );
    long double sinAlpha1  = sin( alpha1 );
    long double cosAlpha1  = cos( alpha1 );

    long double tanU1      = ( 1 - f ) * tan( DEG2RAD( self->latitude ) );
    long double cosU1      = 1 / sqrt( ( 1 + ( tanU1 * tanU1 ) ) );
    long double sinU1      = tanU1 * cosU1;
    long double sigma1     = atan2( tanU1, cosAlpha1 );
    long double sinAlpha   = cosU1 * sinAlpha1;
    long double cosSqAlpha = 1 - ( sinAlpha * sinAlpha );

    long double uSq        = cosSqAlpha * ( ( r_major * r_major ) - ( r_minor * r_minor ) ) / ( r_minor * r_minor );
    long double A          = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
    long double B          = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));

    long double sigma      = s / ( r_minor * A );
    long double sigmaP     = PI * 2;

    long double cos2SigmaM = cos(2*sigma1 + sigma);
    long double sinSigma   = sin(sigma);
    long double cosSigma   = cos(sigma);

    while ( abs(sigma - sigmaP) > 1 * pow( 10, -12 ) ) {
        cos2SigmaM = cos(2*sigma1 + sigma);
        sinSigma   = sin(sigma);
        cosSigma   = cos(sigma);

        long double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
                  B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
        sigmaP                 = sigma;
        sigma                  = s / (r_minor*A) + deltaSigma;
    }

    long double tmp    = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1;
    long double lat2   = atan2( sinU1*cosSigma + cosU1*sinSigma*cosAlpha1, (1-f)*sqrt(sinAlpha*sinAlpha + tmp*tmp) );

    long double lambda = atan2(sinSigma*sinAlpha1, cosU1*cosSigma - sinU1*sinSigma*cosAlpha1);
    long double C      = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
    long double L      = lambda - (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));

    long double lon2   = fmod( DEG2RAD( self->longitude ) + L + ( PI * 3 ), PI * 2 ) - PI;
    long double revAz  = atan2( sinAlpha, -tmp ); // final bearing, if required

    dest->lat = _precision_ld( RAD2DEG( lat2 ), precision );
    dest->lon = _precision_ld( RAD2DEG( lon2 ), precision );
    dest->final_bearing = _precision_ld( RAD2DEG( revAz ), precision );
}


MODULE = Geo::Calc::XS         PACKAGE = Geo::Calc::XS

PROTOTYPES: DISABLE

void new ( char *klass, ... )
    PREINIT:
        int add_count = items - 1;
    PPCODE:
    {
        SV *pv = NEWSV ( 0, sizeof( GCX ) );
        HV *options = newHV();
        int i = 0;
        SvPOK_only( pv );
        if( add_count % 2 != 0 )
            croak( "Please check your parameters while initiating the module\n" );

        for( i = 0; i < add_count; i = i + 2 ) {
            char *key = SvPV_nolen(ST(i + 1));
            if (! (0 == strcmp(key, "lat") || 0 == strcmp(key, "lon") || 0 == strcmp(key, "units")
                || 0 == strcmp(key, "radius"))) {
                croak("Unexpected key '%s'", key);
            }
            hv_store_ent( options, ST( i + 1 ), newSVsv(ST( i + 2)), 0);
        }

        geocalc_init( (GCX *)SvPVX( pv ), options );

        SvREFCNT_dec((SV *) options);

        XPUSHs( sv_2mortal( sv_bless( newRV_noinc( pv ), gv_stashpv( klass, 1 )) ) );
    }

long double
get_lat ( GCX *self )
    CODE:
        RETVAL = self->latitude;
    OUTPUT:
        RETVAL

long double
get_lon ( GCX *self )
    CODE:
        RETVAL = self->longitude;
    OUTPUT:
        RETVAL

char *
get_units ( GCX *self )
    CODE:
        switch ( self->unit_conv)
        {
            case UNIT_METERS:
                RETVAL = "m";
                break;
            case UNIT_KM:
                RETVAL = "k-m";
                break;
            case UNIT_YARDS:
                RETVAL = "yd";
                break;
            case UNIT_FEET:
                RETVAL = "ft";
                break;
            default:
                RETVAL = "mi";
        }
    OUTPUT:
        RETVAL