vhtransform.java

openmap java写的开源数字地图程序. 用applet实现,可以像google map 那样放大缩小地图.

     * 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: [e-mail] 1072 Seena Ave.
     * voice: 415-966-8199 Los Altos, CA. 94024 fax: [phone]
     * ////// Subject: V & H to Latitude and Longitude From:
     * [e-mail] (Col. G.L. Sicherman) Date: 1997/03/05
     * Message-ID: <[e-mail]> 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);
    }
}