update for EOFM-117

?? ?

2008-06-09 bd210ee7438fd203c19d3e8080ea12b79fe56159

 xdgnjobs/ximple-spatialjob/src/main/java/com/ximple/eofms/util/TWDDatumConverter.java

@@ -15,7 +15,7 @@
     *   Definition of math related value
     */
    private static final double COS67_5 = 0.3826834323650897717284599840304e0;
    private static final double PI      = 3.14159265358979323e0;
    private static final double PI = 3.14159265358979323e0;
    private static final double HALF_PI = 1.570796326794896615e0;
    private static final double DEG_RAD = 0.01745329251994329572e0;
    private static final double RAD_DEG = 57.295779513082321031e0;
@@ -23,10 +23,10 @@
    /*
     * Definition of datum related value
     */
    private static final double AD_C       = 1.0026000e0;
    private static final double TWD67_A    = 6378160.0e0;
    private static final double TWD67_B    = 6356774.7192e0;
    private static final double TWD67_ECC  = 0.00669454185458e0;
    private static final double AD_C = 1.0026000e0;
    private static final double TWD67_A = 6378160.0e0;
    private static final double TWD67_B = 6356774.7192e0;
    private static final double TWD67_ECC = 0.00669454185458e0;
    private static final double TWD67_ECC2 = 0.00673966079586e0;

    // different from garmin and already knowned value, but those all value only
@@ -36,17 +36,17 @@
    private static final double TWD67_DY = -361.32e0;

    // release to public. if can got more enough twd67/twd97 control point coordinare,
    private static final double TWD67_DZ   = -180.51e0;
    private static final double TWD67_RX   = -0.00000117e0;    // then we can calculate a better value than now.
    private static final double TWD67_RY   = 0.00000184e0;     //
    private static final double TWD67_RZ   = 0.00000098e0;     // and, also lack twd67/twd97 altitude convertion value...
    private static final double TWD67_S    = 0.00002329e0;     //
    private static final double TWD97_A    = 6378137.0e0;
    private static final double TWD97_B    = 6356752.3141e0;
    private static final double TWD97_ECC  = 0.00669438002290e0;
    private static final double TWD67_DZ = -180.51e0;
    private static final double TWD67_RX = -0.00000117e0;    // then we can calculate a better value than now.
    private static final double TWD67_RY = 0.00000184e0;     //
    private static final double TWD67_RZ = 0.00000098e0;     // and, also lack twd67/twd97 altitude convertion value...
    private static final double TWD67_S = 0.00002329e0;     //
    private static final double TWD97_A = 6378137.0e0;
    private static final double TWD97_B = 6356752.3141e0;
    private static final double TWD97_ECC = 0.00669438002290e0;
    private static final double TWD97_ECC2 = 0.00673949677556e0;
    private static final double TWD67_TM2  = 0.9999e0;         // TWD67->TM2 scale
    private static final double TWD97_TM2  = 0.9999e0;         // TWD97->TM2 scale
    private static final double TWD67_TM2 = 0.9999e0;         // TWD67->TM2 scale
    private static final double TWD97_TM2 = 0.9999e0;         // TWD97->TM2 scale

    /*
     * datum convert function
@@ -59,8 +59,8 @@
        double x1, y1, z1, x2, y2, z2;
        double q, q2, t, t1, s, s1, sum, sin_b, cos_b, sin_p, cos_p;

        lon    = pt.x * DEG_RAD;
        lat    = pt.y * DEG_RAD;
        lon = pt.x * DEG_RAD;
        lat = pt.y * DEG_RAD;
        height = pt.z * DEG_RAD;

        if ((lat < -HALF_PI) && (lat > -1.001 * HALF_PI))
@@ -81,14 +81,14 @@

        sin_lat = Math.sin(lat);
        cos_lat = Math.cos(lat);
        r       = TWD67_A / (Math.sqrt(1.0 - TWD67_ECC * sin_lat * sin_lat));
        x1      = (r + height) * cos_lat * Math.cos(lon);
        y1      = (r + height) * cos_lat * Math.sin(lon);
        z1      = ((r * (1 - TWD67_ECC)) + height) * sin_lat;
        x2      = x1 + TWD67_DX + TWD67_S * (lon + TWD67_RZ * lat - TWD67_RY * height);
        y2      = y1 + TWD67_DY + TWD67_S * (-TWD67_RZ * lon + lat + TWD67_RX * height);
        z2      = z1 + TWD67_DZ + TWD67_S * (TWD67_RY * lon - TWD67_RX * lat + height);
        pole    = 0.0;
        r = TWD67_A / (Math.sqrt(1.0 - TWD67_ECC * sin_lat * sin_lat));
        x1 = (r + height) * cos_lat * Math.cos(lon);
        y1 = (r + height) * cos_lat * Math.sin(lon);
        z1 = ((r * (1 - TWD67_ECC)) + height) * sin_lat;
        x2 = x1 + TWD67_DX + TWD67_S * (lon + TWD67_RZ * lat - TWD67_RY * height);
        y2 = y1 + TWD67_DY + TWD67_S * (-TWD67_RZ * lon + lat + TWD67_RX * height);
        z2 = z1 + TWD67_DZ + TWD67_S * (TWD67_RY * lon - TWD67_RX * lat + height);
        pole = 0.0;

        if (x2 != 0.0)
        {
@@ -104,7 +104,7 @@
            } else
            {
                pole = 1;
                lon  = 0;
                lon = 0;

                if (z2 > 0)
                {
@@ -114,7 +114,7 @@
                    lat = -HALF_PI;
                } else
                {
                    lat  = HALF_PI;
                    lat = HALF_PI;
                    newX = lon * RAD_DEG;
                    newY = lat * RAD_DEG;
                    newZ = -TWD97_B;
@@ -124,18 +124,18 @@
            }
        }

        q2    = x2 * x2 + y2 * y2;
        q     = Math.sqrt(q2);
        t     = z2 * AD_C;
        s     = Math.sqrt(t * t + q2);
        q2 = x2 * x2 + y2 * y2;
        q = Math.sqrt(q2);
        t = z2 * AD_C;
        s = Math.sqrt(t * t + q2);
        sin_b = t / s;
        cos_b = q / s;
        t1    = z2 + TWD97_B * TWD97_ECC2 * sin_b * sin_b * sin_b;
        sum   = q - TWD97_A * TWD97_ECC * cos_b * cos_b * cos_b;
        s1    = Math.sqrt(t1 * t1 + sum * sum);
        t1 = z2 + TWD97_B * TWD97_ECC2 * sin_b * sin_b * sin_b;
        sum = q - TWD97_A * TWD97_ECC * cos_b * cos_b * cos_b;
        s1 = Math.sqrt(t1 * t1 + sum * sum);
        sin_p = t1 / s1;
        cos_p = sum / s1;
        r     = TWD97_A / Math.sqrt(1.0 - TWD97_ECC * sin_p * sin_p);
        r = TWD97_A / Math.sqrt(1.0 - TWD97_ECC * sin_p * sin_p);

        if (cos_p >= COS67_5)
        {
@@ -148,7 +148,7 @@
            height = z2 / sin_p + r * (TWD97_ECC - 1.0);
        }

        if (pole !=0.0)
        if (pole != 0.0)
        {
            lat = Math.atan(sin_p / cos_p);
        }
@@ -167,8 +167,8 @@
        double x1, y1, z1, x2, y2, z2;
        double q, q2, t, t1, s, s1, sum, sin_b, cos_b, sin_p, cos_p;

        lon    = pt.x * DEG_RAD;
        lat    = pt.y * DEG_RAD;
        lon = pt.x * DEG_RAD;
        lat = pt.y * DEG_RAD;
        height = pt.z * DEG_RAD;

        if ((lat < -HALF_PI) && (lat > -1.001 * HALF_PI))
@@ -189,14 +189,14 @@

        sin_lat = Math.sin(lat);
        cos_lat = Math.cos(lat);
        r       = TWD97_A / (Math.sqrt(1.0 - TWD97_ECC * sin_lat * sin_lat));
        x1      = (r + height) * cos_lat * Math.cos(lon);
        y1      = (r + height) * cos_lat * Math.sin(lon);
        z1      = ((r * (1 - TWD97_ECC)) + height) * sin_lat;
        x2      = x1 - TWD67_DX - TWD67_S * (lon + TWD67_RZ * lat - TWD67_RY * height);
        y2      = y1 - TWD67_DY - TWD67_S * (-TWD67_RZ * lon + lat + TWD67_RX * height);
        z2      = z1 - TWD67_DZ - TWD67_S * (TWD67_RY * lon - TWD67_RX * lat + height);
        pole    = 0;
        r = TWD97_A / (Math.sqrt(1.0 - TWD97_ECC * sin_lat * sin_lat));
        x1 = (r + height) * cos_lat * Math.cos(lon);
        y1 = (r + height) * cos_lat * Math.sin(lon);
        z1 = ((r * (1 - TWD97_ECC)) + height) * sin_lat;
        x2 = x1 - TWD67_DX - TWD67_S * (lon + TWD67_RZ * lat - TWD67_RY * height);
        y2 = y1 - TWD67_DY - TWD67_S * (-TWD67_RZ * lon + lat + TWD67_RX * height);
        z2 = z1 - TWD67_DZ - TWD67_S * (TWD67_RY * lon - TWD67_RX * lat + height);
        pole = 0;

        if (x2 != 0.0)
        {
@@ -212,7 +212,7 @@
            } else
            {
                pole = 1;
                lon  = 0;
                lon = 0;

                if (z2 > 0)
                {
@@ -222,7 +222,7 @@
                    lat = -HALF_PI;
                } else
                {
                    lat  = HALF_PI;
                    lat = HALF_PI;
                    newX = lon * RAD_DEG;
                    newY = lat * RAD_DEG;
                    newZ = -TWD67_B;
@@ -232,18 +232,18 @@
            }
        }

        q2    = x2 * x2 + y2 * y2;
        q     = Math.sqrt(q2);
        t     = z2 * AD_C;
        s     = Math.sqrt(t * t + q2);
        q2 = x2 * x2 + y2 * y2;
        q = Math.sqrt(q2);
        t = z2 * AD_C;
        s = Math.sqrt(t * t + q2);
        sin_b = t / s;
        cos_b = q / s;
        t1    = z2 + TWD67_B * TWD67_ECC2 * sin_b * sin_b * sin_b;
        sum   = q - TWD67_A * TWD67_ECC * cos_b * cos_b * cos_b;
        s1    = Math.sqrt(t1 * t1 + sum * sum);
        t1 = z2 + TWD67_B * TWD67_ECC2 * sin_b * sin_b * sin_b;
        sum = q - TWD67_A * TWD67_ECC * cos_b * cos_b * cos_b;
        s1 = Math.sqrt(t1 * t1 + sum * sum);
        sin_p = t1 / s1;
        cos_p = sum / s1;
        r     = TWD67_A / Math.sqrt(1.0 - TWD67_ECC * sin_p * sin_p);
        r = TWD67_A / Math.sqrt(1.0 - TWD67_ECC * sin_p * sin_p);

        if (cos_p >= COS67_5)
        {
@@ -273,24 +273,24 @@
        double n, t, c, A;
        double newX, newY;

        x0   = x * DEG_RAD;
        y0   = y * DEG_RAD;
        x1   = lon * DEG_RAD;
        y1   = lat * DEG_RAD;
        m0   = mercator(y1, a, ecc);
        m1   = mercator(y0, a, ecc);
        n    = a / Math.sqrt(1 - ecc * Math.pow(Math.sin(y0), 2.0));
        t    = Math.pow(Math.tan(y0), 2.0);
        c    = ecc2 * Math.pow(Math.cos(y0), 2.0);
        A    = (x0 - x1) * Math.cos(y0);
        x0 = x * DEG_RAD;
        y0 = y * DEG_RAD;
        x1 = lon * DEG_RAD;
        y1 = lat * DEG_RAD;
        m0 = mercator(y1, a, ecc);
        m1 = mercator(y0, a, ecc);
        n = a / Math.sqrt(1 - ecc * Math.pow(Math.sin(y0), 2.0));
        t = Math.pow(Math.tan(y0), 2.0);
        c = ecc2 * Math.pow(Math.cos(y0), 2.0);
        A = (x0 - x1) * Math.cos(y0);
        newX = scale * n
               * (A + (1.0 - t + c) * A * A * A / 6.0
                  + (5.0 - 18.0 * t + t * t + 72.0 * c - 58.0 * ecc2) * Math.pow(A, 5.0) / 120.0);
                * (A + (1.0 - t + c) * A * A * A / 6.0
                + (5.0 - 18.0 * t + t * t + 72.0 * c - 58.0 * ecc2) * Math.pow(A, 5.0) / 120.0);
        newY = scale
               * (m1 - m0
                  + n * Math.tan(y0)
                    * (A * A / 2.0 + (5.0 - t + 9.0 * c + 4 * c * c) * Math.pow(A, 4.0) / 24.0
                       + (61.0 - 58.0 * t + t * t + 600.0 * c - 330.0 * ecc2) * Math.pow(A, 6.0) / 720.0));
                * (m1 - m0
                + n * Math.tan(y0)
                * (A * A / 2.0 + (5.0 - t + 9.0 * c + 4 * c * c) * Math.pow(A, 4.0) / 24.0
                + (61.0 - 58.0 * t + t * t + 600.0 * c - 330.0 * ecc2) * Math.pow(A, 6.0) / 720.0));
        return new Coordinate(newX, newY);
    }

@@ -300,30 +300,30 @@
        double x0, y0, x1, y1, phi, m, m0, mu, e1;
        double c1, t1, n1, r1, d;

        x0  = x;
        y0  = y;
        x1  = lon * DEG_RAD;
        y1  = lat * DEG_RAD;
        m0  = mercator(y1, a, ecc);
        m   = m0 + y0 / scale;
        e1  = (1.0 - Math.sqrt(1.0 - ecc)) / (1.0 + Math.sqrt(1.0 - ecc));
        mu  = m / (a * (1.0 - ecc / 4.0 - 3.0 * ecc * ecc / 64.0 - 5.0 * ecc * ecc * ecc / 256.0));
        x0 = x;
        y0 = y;
        x1 = lon * DEG_RAD;
        y1 = lat * DEG_RAD;
        m0 = mercator(y1, a, ecc);
        m = m0 + y0 / scale;
        e1 = (1.0 - Math.sqrt(1.0 - ecc)) / (1.0 + Math.sqrt(1.0 - ecc));
        mu = m / (a * (1.0 - ecc / 4.0 - 3.0 * ecc * ecc / 64.0 - 5.0 * ecc * ecc * ecc / 256.0));
        phi = mu + (3.0 * e1 / 2.0 - 27.0 * Math.pow(e1, 3.0) / 32.0) * Math.sin(2.0 * mu)
              + (21.0 * e1 * e1 / 16.0 - 55.0 * Math.pow(e1, 4.0) / 32.0) * Math.sin(4.0 * mu)
              + 151.0 * Math.pow(e1, 3.0) / 96.0 * Math.sin(6.0 * mu) + 1097.0 * Math.pow(e1, 4.0) / 512.0 * Math.sin(8.0 * mu);
        c1   = ecc2 * Math.pow(Math.cos(phi), 2.0);
        t1   = Math.pow(Math.tan(phi), 2.0);
        n1   = a / Math.sqrt(1 - ecc * Math.pow(Math.sin(phi), 2.0));
        r1   = a * (1.0 - ecc) / Math.pow(1.0 - ecc * Math.pow(Math.sin(phi), 2.0), 1.5);
        d    = x0 / (n1 * scale);
                + (21.0 * e1 * e1 / 16.0 - 55.0 * Math.pow(e1, 4.0) / 32.0) * Math.sin(4.0 * mu)
                + 151.0 * Math.pow(e1, 3.0) / 96.0 * Math.sin(6.0 * mu) + 1097.0 * Math.pow(e1, 4.0) / 512.0 * Math.sin(8.0 * mu);
        c1 = ecc2 * Math.pow(Math.cos(phi), 2.0);
        t1 = Math.pow(Math.tan(phi), 2.0);
        n1 = a / Math.sqrt(1 - ecc * Math.pow(Math.sin(phi), 2.0));
        r1 = a * (1.0 - ecc) / Math.pow(1.0 - ecc * Math.pow(Math.sin(phi), 2.0), 1.5);
        d = x0 / (n1 * scale);
        newX = (x1 + (d - (1.0 + 2.0 * t1 + c1) * Math.pow(d, 3.0) / 6.0
                      + (5.0 - 2.0 * c1 + 28.0 * t1 - 3.0 * c1 * c1 + 8.0 * ecc2 + 24.0 * t1 * t1) * Math.pow(d, 5.0)
                        / 120.0) / Math.cos(phi)) * RAD_DEG;
                + (5.0 - 2.0 * c1 + 28.0 * t1 - 3.0 * c1 * c1 + 8.0 * ecc2 + 24.0 * t1 * t1) * Math.pow(d, 5.0)
                / 120.0) / Math.cos(phi)) * RAD_DEG;
        newY = (phi
                - n1 * Math.tan(phi) / r1
                  * (d * d / 2.0 - (5.0 + 3.0 * t1 + 10.0 * c1 - 4.0 * c1 * c1 - 9.0 * ecc2) * Math.pow(d, 4.0) / 24.0
                     + (61.0 + 90.0 * t1 + 298.0 * c1 + 45.0 * t1 * t1 - 252.0 * ecc2 - 3.0 * c1 * c1) * Math.pow(d, 6.0)
                       / 72.0)) * RAD_DEG;
                * (d * d / 2.0 - (5.0 + 3.0 * t1 + 10.0 * c1 - 4.0 * c1 * c1 - 9.0 * ecc2) * Math.pow(d, 4.0) / 24.0
                + (61.0 + 90.0 * t1 + 298.0 * c1 + 45.0 * t1 * t1 - 252.0 * ecc2 - 3.0 * c1 * c1) * Math.pow(d, 6.0)
                / 72.0)) * RAD_DEG;
        return new Coordinate(newX, newY);
    }

@@ -335,9 +335,9 @@
        } else
        {
            return a * ((1.0 - ecc / 4.0 - 3.0 * ecc * ecc / 64.0 - 5.0 * ecc * ecc * ecc / 256.0) * y
                        - (3.0 * ecc / 8.0 + 3.0 * ecc * ecc / 32.0 + 45.0 * ecc * ecc * ecc / 1024.0) * Math.sin(2.0 * y)
                        + (15.0 * ecc * ecc / 256.0 + 45.0 * ecc * ecc * ecc / 1024.0) * Math.sin(4.0 * y)
                        - (35.0 * ecc * ecc * ecc / 3072.0) * Math.sin(6.0 * y));
                    - (3.0 * ecc / 8.0 + 3.0 * ecc * ecc / 32.0 + 45.0 * ecc * ecc * ecc / 1024.0) * Math.sin(2.0 * y)
                    + (15.0 * ecc * ecc / 256.0 + 45.0 * ecc * ecc * ecc / 1024.0) * Math.sin(4.0 * y)
                    - (35.0 * ecc * ecc * ecc / 3072.0) * Math.sin(6.0 * y));
        }
    }

@@ -492,6 +492,7 @@

    /**
     * ¥ÑTM2®y¼ÐÂà´«¦ÜTWD97®y¼Ð
     *
     * @param pt TM2¦ì¸m
     * @return ·sªºTWD97®y¼Ð
     */