* efficient (non-iterative) inverse for this program! (i.e. a     * program to compute lat & long from V & H).     */    /** lat and lon are in degrees, positive north and east. */    public void toVH(double lat, double lon) {        lat = radians(lat);        lon = radians(lon);        /* Translate east by 52 degrees */        double lon1 = lon + radians(52.0);        /* Convert latitude to geocentric latitude using Horner's rule */        double latsq = lat * lat;        double lat1 = lat                * (K1 + (K2 + (K3 + (K4 + K5 * latsq) * latsq) * latsq) * latsq);        /*         * x, y, and z are the spherical coordinates corresponding to         * lat, lon.         */        double cos_lat1 = Math.cos(lat1);        double x = cos_lat1 * Math.sin(-lon1);        double y = cos_lat1 * Math.cos(-lon1);        double z = Math.sin(lat1);        /*         * e and w are the cosine of the angular distance (radians)         * between our point and the east and west centers.         */        double e = EX * x + EY * y + EZ * z;        double w = WX * x + WY * y + WZ * z;        e = e > 1.0 ? 1.0 : e;        w = w > 1.0 ? 1.0 : w;        e = M_PI_2 - Math.atan(e / Math.sqrt(1 - e * e));        w = M_PI_2 - Math.atan(w / Math.sqrt(1 - w * w));        /* e and w are now in radians. */        double ht = (e * e - w * w + .16) / .8;        double vt = Math.sqrt(Math.abs(e * e - ht * ht));        vt = (PX * x + PY * y + PZ * z) < 0 ? -vt : vt;        /* rotate and translate to get final v and h. */        double v = TRANSV + K9 * ht - K10 * vt;        double h = TRANSH + K10 * ht + K9 * vt;        this.resultV = v;        this.resultH = h;    }    /*     * Stu Jeffery Internet: stu@shell.portal.com 1072 Seena Ave.     * voice: 415-966-8199 Los Altos, CA. 94024 fax: 415-966-8456     * ////// Subject: V & H to Latitude and Longitude From:     * sicherman@lucent.com (Col. G.L. Sicherman) Date: 1997/03/05     * Message-ID: <telecom17.57.10@massis.lcs.mit.edu> Newsgroups:     * comp.dcom.telecom [More Headers] [Subscribe to     * comp.dcom.telecom]     *      *      * Recently I wanted to convert some Bell Labs "V&H" coordinates     * to latitude and longitude. A careful search through the     * Telecomm- unications Archives turned up a C program for     * converting in the other direction, and many pleas for what I     * was looking for. One poster even offered money!     *      * Since I work for Bell Labs, I had no trouble getting a copy of     * Erik Grimmelmann's legendary memorandum. (Don't get your hopes     * up - Bell Labs has no intention of releasing it to the public!)     * Thus armed, I hacked up the following C program, which ought to     * compile on any C platform. Its input and output agree with the     * output and input of ll_to_vh (as hacked by Tom Libert), and the     * comments summarize the math as explained by Grimmelmann. Enjoy!     */    /**     * 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.     * <P>     *      * 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.     * <P>     *      * 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.     * <P>     *      * 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     * <P>     * h = (square of distance to E - square of distance to W + square     * of distance between E and W) / twice distance between E and W;     * <BR>     * v = square root of absolute value of (square of distance to E -     * square of h).     * <P>     * Reduce by three-tenths of a percent, rotate by 76.597497     * degrees, and add 6363.235 to V and 2250.7 to H.     * <P>     *      * 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     * <P>     * A = (ExWz-EzWx)^2 + (EyWzx-EzWy)^2 + (ExWy-EyWx)^2; <BR>     * 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)]; <BR>     * 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.     * <P>     * Solve with the quadratic formula. The latitude is simply the     * arc sine of Pz. Px and Py satisfy     * <P>     * ExPx + EyPy + EzPz = cos(arc to E); <BR>     * WxPx + WyPy + WzPz = cos(arc to W).     * <P>     * 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.     */    public void toLatLon(double v0, double h0) {        /* GX = ExWz - EzWx; GY = EyWz - EzWy */        final double GX = 0.216507961908834992;        final double GY = -0.134633014879368199;        /* A = (ExWz-EzWx)^2 + (EyWz-EzWy)^2 + (ExWy-EyWx)^2 */        final double A = 0.151646645621077297;        /* Q = ExWy-EyWx; Q2 = Q*Q */        final double Q = -0.294355056616412800;        final double Q2 = 0.0866448993556515751;        final double EPSILON = .0000001;        double v = (double) v0;        double h = (double) h0;        double t1 = (v - TRANSV) / RADIUS;        double t2 = (h - TRANSH) / RADIUS;        double vhat = ROTC * t2 - ROTS * t1;        double hhat = ROTS * t2 + ROTC * t1;        double e = Math.cos(Math.sqrt(vhat * vhat + hhat * hhat));        double w = Math.cos(Math.sqrt(vhat * vhat + (hhat - 0.4) * (hhat - 0.4)));        double fx = EY * w - WY * e;        double fy = EX * w - WX * e;        double b = fx * GX + fy * GY;        double c = fx * fx + fy * fy - Q2;        double disc = b * b - A * c; /* discriminant */        double x, y, z, delta;        if (Math.abs(disc) < EPSILON) {            //      if (disc==0.0) { /* It's right on the E-W axis */            z = b / A;            x = (GX * z - fx) / Q;            y = (fy - GY * z) / Q;        } else {            delta = Math.sqrt(disc);            z = (b + delta) / A;            x = (GX * z - fx) / Q;            y = (fy - GY * z) / Q;            if (vhat * (PX * x + PY * y + PZ * z) < 0) { /*                                                          * wrong                                                          * direction                                                          */                z = (b - delta) / A;                x = (GX * z - fx) / Q;                y = (fy - GY * z) / Q;            }        }        double lat = Math.asin(z);        /*         * Use polynomial approximation for inverse mapping (sphere to         * spheroid):         */        final double[] bi = { 1.00567724920722457, -0.00344230425560210245,                0.000713971534527667990, -0.0000777240053499279217,                0.00000673180367053244284, -0.000000742595338885741395,                0.0000000905058919926194134 };        double lat2 = lat * lat;        /*         * KRA: Didn't seem to work at first, so i unrolled it. double         * earthlat = 0.0; for (int i=6; i>=0; i--) { earthlat =         * (earthlat + bi[i]) * (i > 0? lat2 : lat); }         */        double earthlat = lat                * (bi[0] + lat2                        * (bi[1] + lat2                                * (bi[2] + lat2                                        * (bi[3] + lat2                                                * (bi[4] + lat2                                                        * (bi[5] + lat2                                                                * (bi[6])))))));        earthlat = degrees(earthlat);        /*         * Adjust longitude by 52 degrees:         */        double lon = degrees(Math.atan2(x, y));        double earthlon = lon + 52.0;        this.resultLat = earthlat;        this.resultLon = -earthlon;        // Col. G. L. Sicherman.    }    public void toLatLon(int v0, int h0) {        toLatLon((double) v0, (double) h0);    }}