polarcoordconverter.java

来自「world wind java sdk 源码」· Java 代码 · 共 467 行 · 第 1/2 页

        Polar_Delta_Northing = Northing * 2; // Increased range for accepted easting and northing values
        Polar_Delta_Northing = Math.abs(Polar_Delta_Northing) + epsilon;
        Polar_Delta_Easting = Polar_Delta_Northing;

        return (Error_Code);
    }

    /**
     * The function Convert_Geodetic_To_Polar_Stereographic converts geodetic coordinates (latitude and longitude) to
     * Polar Stereographic coordinates (easting and northing), according to the current ellipsoid and Polar
     * Stereographic projection parameters. If any errors occur, error code(s) are returned by the function, otherwise
     * POLAR_NO_ERROR is returned.
     *
     * @param Latitude  latitude, in radians
     * @param Longitude Longitude, in radians
     * @return error code
     */
    public long convertGeodeticToPolarStereographic(double Latitude, double Longitude)
    {
        double dlam;
        double slat;
        double essin;
        double t;
        double rho;
        double pow_es;
        long Error_Code = POLAR_NO_ERROR;

        if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
        {   /* Latitude out of range */
            Error_Code |= POLAR_LAT_ERROR;
        }
        if ((Latitude < 0) && (Southern_Hemisphere == 0))
        {   /* Latitude and Origin Latitude in different hemispheres */
            Error_Code |= POLAR_LAT_ERROR;
        }
        if ((Latitude > 0) && (Southern_Hemisphere == 1))
        {   /* Latitude and Origin Latitude in different hemispheres */
            Error_Code |= POLAR_LAT_ERROR;
        }
        if ((Longitude < -PI) || (Longitude > TWO_PI))
        {  /* Longitude out of range */
            Error_Code |= POLAR_LON_ERROR;
        }

        if (Error_Code == POLAR_NO_ERROR)
        {  /* no errors */

            if (Math.abs(Math.abs(Latitude) - PI_OVER_2) < 1.0e-10)
            {
                Easting = 0.0;
                Northing = 0.0;
            } else
            {
                if (Southern_Hemisphere != 0)
                {
                    Longitude *= -1.0;
                    Latitude *= -1.0;
                }
                dlam = Longitude - Polar_Origin_Long;
                if (dlam > PI)
                {
                    dlam -= TWO_PI;
                }
                if (dlam < -PI)
                {
                    dlam += TWO_PI;
                }
                slat = Math.sin(Latitude);
                essin = es * slat;
                pow_es = Math.pow((1.0 - essin) / (1.0 + essin), es_OVER_2);
                t = Math.tan(PI_Over_4 - Latitude / 2.0) / pow_es;

                if (Math.abs(Math.abs(Polar_Origin_Lat) - PI_OVER_2) > 1.0e-10)
                    rho = Polar_a_mc * t / tc;
                else
                    rho = two_Polar_a * t / e4;


                if (Southern_Hemisphere != 0)
                {
                    Easting = -(rho * Math.sin(dlam) - Polar_False_Easting);
                    //Easting *= -1.0;
                    Northing = rho * Math.cos(dlam) + Polar_False_Northing;
                } else
                    Easting = rho * Math.sin(dlam) + Polar_False_Easting;
                    Northing = -rho * Math.cos(dlam) + Polar_False_Northing;
            }
        }
        return (Error_Code);
    }

    public double getEasting()
    {
        return Easting;
    }

    public double getNorthing()
    {
        return Northing;
    }

    /**
     *  The function Convert_Polar_Stereographic_To_Geodetic converts Polar
     *  Stereographic coordinates (easting and northing) to geodetic
     *  coordinates (latitude and longitude) according to the current ellipsoid
     *  and Polar Stereographic projection Parameters. If any errors occur, the
     *  code(s) are returned by the function, otherwise POLAR_NO_ERROR
     *  is returned.
     *
     *  @param Easting Easting (X), in meters
     *  @param Northing Northing (Y), in meters
     *  @return error code
     */
    public long convertPolarStereographicToGeodetic (double Easting, double Northing)
    {
        double dy = 0, dx = 0;
        double rho = 0;
        double t;
        double PHI, sin_PHI;
        double tempPHI = 0.0;
        double essin;
        double pow_es;
        double delta_radius;
        long Error_Code = POLAR_NO_ERROR;
        double min_easting = Polar_False_Easting - Polar_Delta_Easting;
        double max_easting = Polar_False_Easting + Polar_Delta_Easting;
        double min_northing = Polar_False_Northing - Polar_Delta_Northing;
        double max_northing = Polar_False_Northing + Polar_Delta_Northing;

        if (Easting > max_easting || Easting < min_easting)
        { /* Easting out of range */
            Error_Code |= POLAR_EASTING_ERROR;
        }
        if (Northing > max_northing || Northing < min_northing)
        { /* Northing out of range */
            Error_Code |= POLAR_NORTHING_ERROR;
        }

        if (Error_Code == POLAR_NO_ERROR)
        {
            dy = Northing - Polar_False_Northing;
            dx = Easting - Polar_False_Easting;

            /* Radius of point with origin of false easting, false northing */
            rho = Math.sqrt(dx * dx + dy * dy);

            delta_radius = Math.sqrt(Polar_Delta_Easting * Polar_Delta_Easting + Polar_Delta_Northing * Polar_Delta_Northing);

            if(rho > delta_radius)
            { /* Point is outside of projection area */
                Error_Code |= POLAR_RADIUS_ERROR;
            }
        }

        if (Error_Code == POLAR_NO_ERROR)
        { /* no errors */
            if ((dy == 0.0) && (dx == 0.0))
            {
                Latitude = PI_OVER_2;
                Longitude = Polar_Origin_Long;

            }
            else
            {
                if (Southern_Hemisphere != 0)
                {
                    dy *= -1.0;
                    dx *= -1.0;
                }

                if (Math.abs(Math.abs(Polar_Origin_Lat) - PI_OVER_2) > 1.0e-10)
                    t = rho * tc / (Polar_a_mc);
                else
                    t = rho * e4 / (two_Polar_a);
                PHI = PI_OVER_2 - 2.0 * Math.atan(t);
                while (Math.abs(PHI - tempPHI) > 1.0e-10)
                {
                    tempPHI = PHI;
                    sin_PHI = Math.sin(PHI);
                    essin =  es * sin_PHI;
                    pow_es = Math.pow((1.0 - essin) / (1.0 + essin), es_OVER_2);
                    PHI = PI_OVER_2 - 2.0 * Math.atan(t * pow_es);
                }
                Latitude = PHI;
                Longitude = Polar_Origin_Long + Math.atan2(dx, -dy);

                if (Longitude > PI)
                    Longitude -= TWO_PI;
                else if (Longitude < -PI)
                    Longitude += TWO_PI;


                if (Latitude > PI_OVER_2)  /* force distorted values to 90, -90 degrees */
                    Latitude = PI_OVER_2;
                else if (Latitude < -PI_OVER_2)
                    Latitude = -PI_OVER_2;

                if (Longitude > PI)  /* force distorted values to 180, -180 degrees */
                    Longitude = PI;
                else if (Longitude < -PI)
                    Longitude = -PI;

            }
            if (Southern_Hemisphere != 0)
            {
                Latitude *= -1.0;
                Longitude *= -1.0;
            }

        }
        return (Error_Code);
    } 

    /**
     * @return Latitude in radians.
     */
    public double getLatitude()
    {
        return Latitude;
    }

    /**
     * @return Longitude in radians.
     */
    public double getLongitude()
    {
        return Longitude;
    }

}