view release on metacpan or  search on metacpan

XS.c  view on Meta::CPAN

    /* e and w are now in radians. */
    ht = (e*e - w*w + .16)/.8;
    vt = sqrt(fabs(e*e - ht*ht));
    vt = (PX*x + PY*y + PZ*z) < 0 ? -vt : vt;

    /* rotate and translate to get final v and h. */
    v = TRANSV +  RADIUS*ROTC*ht - RADIUS*ROTS*vt;
    h = TRANSH + RADIUS*ROTS*ht +  RADIUS*ROTC*vt;

    XPUSHs(sv_2mortal(newSVnv(v)));
    XPUSHs(sv_2mortal(newSVnv(h)));
    PUTBACK;

}

/**
 *      V&H is a system of coordinates (V and H) for describing
 *      locations of rate centers in the United States.  The
 *      projection, devised by J. K. Donald, is an "elliptical,"
 *      or "doubly equidistant" projection, scaled down by a factor
 *      of 0.003 to balance errors.
 *
 *      The foci of the projection, from which distances are
 *      measured accurately (except for the scale correction),
 *      are at 37d 42m 14.69s N, 82d 39m 15.27s W (in Floyd Co.,
 *      Ky.) and 41d 02m 55.53s N, 112d 03m 39.35 W (in Webster
 *      Co., Utah).  They are just 0.4 radians apart.
 *
 *      Here is the transformation from latitude and longitude to V&H:
 *      First project the earth from its ellipsoidal surface
 *      to a sphere.  This alters the latitude; the coefficients
 *      bi in the program are the coefficients of the polynomial
 *      approximation for the inverse transformation.  (The
 *      function is odd, so the coefficients are for the linear
 *      term, the cubic term, and so on.)  Also subtract 52 degrees
 *      from the longitude.
 *
 *      For the rest, compute the arc distances of the given point
 *      to the reference points, and transform them to the coordinate
 *      system in which the line through the reference points is the
 *      X-axis and the origin is the eastern reference point.
 *      The solution is
 *              h = (square of distance to E - square of distance to W
 *                      + square of distance between E and W) /
 *                      twice distance between E and W;
 *              v = square root of absolute value of (square of
 *                      distance to E - square of h).
 *      Reduce by three-tenths of a percent, rotate by 76.597497
 *      degrees, and add 6363.235 to V and 2250.7 to H.
 *
 *      To go the other way, as this program does, undo the final translation,
 *      rotation, and scaling.  The z-value Pz of the point on the x-y-z sphere
 *      satisfies the quadratic Azz+Bz+c=0, where
 *              A = (ExWz-EzWx)^2 + (EyWzx-EzWy)^2 + (ExWy-EyWx)^2;
 *              B = -2[(Ex cos(arc to W) - Wx cos(arc to E))(ExWz-EzWx) -
 *                      (Ey cos(arc to W) -Wy cos(arc to E))(EyWz-EzWy)];
 *              C = (Ex cos(arc to W) - Wx cos(arc to E))^2 +
 *                      (Ey cos(arc to W) - Wy cos(arc to E))^2 -
 *                      (ExWy - EyWx)^2.
 *      Solve with the quadratic formula.  The latitude is simply the
 *      arc sine of Pz.  Px and Py satisfy
 *              ExPx + EyPy + EzPz = cos(arc to E);
 *              WxPx + WyPy + WzPz = cos(arc to W).
 *      Substitute Pz's value, and solve linearly to get Px and Py.
 *      The longitude is the arc tangent of Px/Py.
 *      Finally, this latitude and longitude are spherical; use the
 *      inverse polynomial approximation on the latitude to get the
 *      ellipsoidal earth latitude, and add 52 degrees to the longitude.
 */

void toLatLon (double v, double h) {

    double GX,GY,A,Q,Q2,EPSILON,t1,t2,vhat,hhat,e,w,fx,fy,b,c,disc,x,y,z;
    double delta,lat,lat2,earthlat,lon,earthlon;
    /*
     *  Use polynomial approximation for inverse mapping (sphere to spheroid)
     */
    double bi[] = {
        1.00567724920722457,
        -0.00344230425560210245,
        0.000713971534527667990,
        -0.0000777240053499279217,
        0.00000673180367053244284,
        -0.000000742595338885741395,
        0.0000000905058919926194134
    };
    dXSARGS;
    sp = mark;

    /* GX = ExWz - EzWx; GY = EyWz - EzWy */
    GX =  0.216507961908834992;
    GY = -0.134633014879368199;
    /* A = (ExWz-EzWx)^2 + (EyWz-EzWy)^2 + (ExWy-EyWx)^2 */
    A =   0.151646645621077297;
    /* Q = ExWy-EyWx; Q2 = Q*Q */
    Q =  -0.294355056616412800;
    Q2=   0.0866448993556515751;
    EPSILON = .0000001;

    t1 = (v - TRANSV) / RADIUS;
    t2 = (h - TRANSH) / RADIUS;
    vhat = ROTC*t2 - ROTS*t1;
    hhat = ROTS*t2 + ROTC*t1;
    e = cos(sqrt(vhat*vhat + hhat*hhat));
    w = cos(sqrt(vhat*vhat + (hhat-0.4)*(hhat-0.4)));
    fx = EY*w - WY*e;
    fy = EX*w - WX*e;
    b = fx*GX + fy*GY;
    c = fx*fx + fy*fy - Q2;
    disc = b*b - A*c;               /* discriminant */
    x, y, z, delta;
    if (fabs(disc) < EPSILON) {
        z = b/A;
        x = (GX*z - fx)/Q;
        y = (fy - GY*z)/Q;
    } else {
        delta = sqrt(disc);
        z = (b + delta)/A;
        x = (GX*z - fx)/Q;
        y = (fy - GY*z)/Q;