<|start_filename|>library/src/main/java/org/opensextant/geodesy/UPS.java<|end_filename|>
/****************************************************************************************
 *  UPS.java
 *
 *  Created: Jan 20, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.io.Serializable;
import java.text.DecimalFormat;

/**
 * UPS coordinates are based on a family of two Polar Stereographic map projections,
 * one for each pole. The origin of the UPS coordinate system is the pole (north or
 * south), where X=2,000,000m and Y=2,000,000m.
 * <p/>
 * The X-axis lies along the meridians 90&deg;E and 90&deg;W.
 * Moving from the pole (north or south), X-values (Eastings) increase along
 * the 90&deg;E meridian.
 * The Y-axis lies along the meridians 0&deg; and 180&deg;.
 * Moving from the north pole, Y-values (Northings) increase along the 180&deg;
 * meridian.  Moving from the south pole, Y-values (Northings) increase along
 * the 0&deg; meridian.
 */
public class UPS implements Serializable {
	private static final long serialVersionUID = 1L;
	
	// Static Class Constants
    public static final double MIN_NORTH_LATDEG = +83.5;
    public static final double MAX_SOUTH_LATDEG = -79.5;
    public static final int FALSE_EASTING = 2000000;
    public static final int FALSE_NORTHING = 2000000;

    private static final double MIN_EASTING = 0.0;
    private static final double MAX_EASTING = 4000000.0;
    private static final double MIN_NORTHING = 0.0;
    private static final double MAX_NORTHING = 4000000.0;
    private static final Latitude NORTH_POLAR = new Latitude(+81, 6, 52.3);
    private static final Latitude SOUTH_POLAR = new Latitude(-81, 6, 52.3);

    // Static Class Methods

    /**
     * Validate UPS polar zone to see if it is valid (by itself). Used in place of UTM
     * lat band char.
     *
     * @param polarZone UPS Polar Zone ('A' or 'B' for Southern, 'Y' or 'Z' for Northern)
     * @throws IllegalArgumentException error if UPS Polar Zone is invalid
     */
    public static void validatePolarZone(char polarZone) {
        if ((polarZone != 'A') && (polarZone != 'B') && (polarZone != 'Y') && (polarZone != 'Z'))
            throw new IllegalArgumentException("UPS polar zone '" + polarZone + "' is not valid");
    }

    /*
     * Validate hemisphere character to see if it is valid (by itself).
     *
     * @param hemisphere UPS Hemisphere character ('N' for North, 'S' for South)
     * @throws IllegalArgumentException error if hemisphere character is invalid
     */
    private static void validateHemisphere(char hemisphere) {
        if ((hemisphere != 'N') && (hemisphere != 'S'))
            throw new IllegalArgumentException("Invalid hemisphere '" +
                    hemisphere + "', should be 'N' or 'S'");
    }

    /*
     * Validate easting value to see if it is within allowed positive numerical range.
     *
     * @param easting UPS Easting value in meters
     * @throws IllegalArgumentException error if easting value is out of range
     */
    private static void validateEasting(double easting) {
        if ((easting < MIN_EASTING) || (MAX_EASTING < easting))
            throw new IllegalArgumentException("Easting value '" + easting +
                    "' is outside of valid range (0 to 4,000,000 meters)");
    }

    /*
     * Validate northing value to see if it is within allowed positive numerical range.
     *
     * @param northing UPS Northing value in meters
     * @throws IllegalArgumentException error if northing value is out of range
     */
    private static void validateNorthing(double northing) {
        if ((northing < MIN_NORTHING) || (MAX_NORTHING < northing))
            throw new IllegalArgumentException("Northing value '" + northing +
                    "' is outside of valid range (0 to 4,000,000 meters)");
    }

    /*
     * Validate latitude value to see if it is within allowed numerical range.
     *
     * @param lat - latitude
     * @throws IllegalArgumentException error if latitude value is out of range
     */
    private static void validateLatitude(Latitude lat) {
        double latDeg = lat.inDegrees();
        if ((0.0 <= latDeg) && (latDeg < MIN_NORTH_LATDEG)) {
            throw new IllegalArgumentException("Latitude value '" + latDeg +
                    "' is out of legal range (+83.5 deg to +90 deg) for UPS Northern Hemisphere");
        } else if ((MAX_SOUTH_LATDEG < latDeg) && (latDeg < 0.0)) {
            throw new IllegalArgumentException("Latitude value '" + latDeg +
                    "' is out of legal range (-79.5 deg to -90 deg) for UPS Southern Hemisphere");
        }
    }

    // *************************** End of Static Definitions ********************************

    // Instance Variables
    @NonNull private PolarStereographic ps;  // Polar Stereographic projection object (keeps ellipsoid)
    private char hemisphere;        // Hemisphere char ('N' for Northern, 'S' for Southern)
    private char polarZone;         // East-West zone identifier character (like UTM lat band)
    private double easting;         // positive meters east of false adjusted pole
    private double northing;        // positive meters north of false adjusted pole

    @NonNull private Geodetic2DPoint lonLat; // Longitude & Latitude coordinates for this UPS position

    // private initializer method for geodetic constructors
    private void initGeodetic(Ellipsoid ellip, Longitude lon, Latitude lat)
            throws IllegalArgumentException {
        hemisphere = lat.getHemisphere();
        Latitude latOfTrueScale = (hemisphere == 'N') ? NORTH_POLAR : SOUTH_POLAR;
        ps = new PolarStereographic(latOfTrueScale);
        if (ellip != null) ps.setEllipsoid(ellip);
        lonLat = new Geodetic2DPoint(lon, lat);
        fromGeodetic();
    }

    // private initializer method for UPS parameter constructors
    private void initUPS(Ellipsoid ellip, char hemisphere, double easting, double northing)
            throws IllegalArgumentException {
        Latitude latOfTrueScale = (hemisphere == 'N') ? NORTH_POLAR : SOUTH_POLAR;
        ps = new PolarStereographic(latOfTrueScale);
        if (ellip != null) ps.setEllipsoid(ellip);
        this.hemisphere = hemisphere;
        this.easting = easting;
        this.northing = northing;
        // Determine polar zone and cache value in instance variable
        if (hemisphere == 'S') polarZone = (easting < 0.0) ? 'A' : 'B';
        else polarZone = (easting < 0.0) ? 'Y' : 'Z';
        toGeodetic();
    }

    /**
     * This constructor takes an ellipsoid and Geodetic2DPoint object.  If WGS 84 is desired,
     * the caller should instead use one of the constructors without the ellipsoid parameter,
     * since WGS 84 is the default. The specified Geodetic2DPoint object (lon-lat point) is
     * converted to its UPS equivalent.
     *
     * @param ellip  Ellipsoid data model for earth
     * @param lonLat Geodetic2DPoint coordinates (lon-lat) to be converted to UPS
     * @throws IllegalArgumentException error if arguments are invalid or out of legal range
     */
    public UPS(Ellipsoid ellip, Geodetic2DPoint lonLat) {
        initGeodetic(ellip, lonLat.getLongitude(), lonLat.getLatitude());
    }

    /**
     * This constructor takes an ellipsoid, Longitude and Latitude Objects.
     * If WGS 84 is desired, the caller should instead use one of the constructors without
     * the ellipsoid parameter, since WGS 84 is the default. The specified longitude and
     * latitude object values are converted to their UPS equivalent.
     *
     * @param ellip Ellipsoid data model for earth
     * @param lon   Longitude of point to convert to UPS
     * @param lat   Latitude of point to convert to UPS
     * @throws IllegalArgumentException error if arguments are invalid or out of legal range
     */
    public UPS(Ellipsoid ellip, Longitude lon, Latitude lat) {
        initGeodetic(ellip, lon, lat);
    }

    /**
     * This constructor takes an ellipsoid and the UPS parameters for hemisphere, easting, and
     * northing.  If WGS 84 is desired, the caller should instead use one of the constructors
     * without the ellipsoid parameter, since WGS 84 is the default.  The UPS parameters are
     * converted to their equivalent geodetic (lon-lat) position, along with the correct polar
     * zone.
     *
     * @param ellip      Ellipsoid data model for earth
     * @param hemisphere character 'N' for Northern or 'S' for Southern hemisphere
     * @param easting    positive meters east of false adjusted pole
     * @param northing   positive meters north of false adjusted pole
     * @throws IllegalArgumentException error if arguments are invalid or out of legal range
     */
    public UPS(Ellipsoid ellip, char hemisphere, double easting, double northing)
            throws IllegalArgumentException {
        initUPS(ellip, hemisphere, easting, northing);
    }

    /**
     * This constructor assumes the WGS 84 Ellipsoid (default for PolarStereographic).
     * The specified Geodetic2DPoint object (lon-lat point) is converted to its UPS equivalent.
     *
     * @param lonLat Geodetic2DPoint coordinates (lon-lat) to be converted to UPS
     * @throws IllegalArgumentException error if arguments are invalid or out of legal range
     */
    public UPS(Geodetic2DPoint lonLat) {
        initGeodetic(null, lonLat.getLongitude(), lonLat.getLatitude());
    }

    /**
     * This constructor assumes the WGS 84 Ellipsoid (default for PolarStereographic).
     * The specified longitude and latitude object values are converted to their UPS equivalent.
     *
     * @param lon Longitude of point to convert to UPS
     * @param lat Latitude of point to convert to UPS
     * @throws IllegalArgumentException error if latitude is not close enough to a pole
     */
    public UPS(Longitude lon, Latitude lat) {
        initGeodetic(null, lon, lat);
    }

    /**
     * This constructor assumes the WGS 84 Ellipsoid (default for PolarStereographic).
     * The UPS parameters are converted to their equivalent geodetic (lon-lat) position,
     * along with the correct polar zone.
     *
     * @param hemisphere character 'N' for Northern or 'S' for Southern hemisphere
     * @param easting    positive meters east of false adjusted central meridian for lonZone
     * @param northing   positive meters north of false adjusted origin (equator for 'N' hemisphere)
     * @throws IllegalArgumentException error if arguments are invalid or out of legal range
     */
    public UPS(char hemisphere, double easting, double northing) {
        initUPS(null, hemisphere, easting, northing);
    }

    /*
     * This method converts UPS projection (hemisphere, easting and northing)
     * coordinates to geodetic (longitude and latitude) coordinates, based on the
     * current ellipsoid model.
     *
     * @throws IllegalArgumentException input parameter error(s)
     */
    private void toGeodetic() {
        // Validate input parameters
        validateHemisphere(hemisphere);
        validateEasting(easting);
        validateNorthing(northing);

        double e = easting - FALSE_EASTING;
        double n = northing - FALSE_NORTHING;

        // Un-project to geodetic coordinates
        lonLat = ps.toGeodetic(e, n);
    }

    /*
     * This method converts from Geodetic2DPoint coordinates (lon-lat) to UPS parameters.
     */
    private void fromGeodetic() {
        Longitude lon = lonLat.getLongitude();
        Latitude lat = lonLat.getLatitude();

        // Validate latitude
        validateLatitude(lat);
        hemisphere = lat.getHemisphere();

        Topocentric2DPoint en = ps.toPolarStereographic(lon, lat);
        easting = en.getEasting() + UPS.FALSE_EASTING;
        northing = en.getNorthing() + UPS.FALSE_NORTHING;

        // Determine polar zone and cache value in instance variable
        if (hemisphere == 'S') polarZone = (easting < UPS.FALSE_EASTING) ? 'A' : 'B';
        else polarZone = (easting < UPS.FALSE_EASTING) ? 'Y' : 'Z';
    }

    /**
     * This accessor method is used to get the currently set Ellipsoid earth model for
     * making projections during conversion to and from geodetic coordinates.
     *
     * @return the currently set Ellipsoid earth model being used by this UPS coordinate
     */
    public Ellipsoid getEllipsoid() {
        return ps.getEllipsoid();
    }

    /**
     * This accessor method returns the geodetic (lon-lat) point for this UPS coordinate.
     *
     * @return Geodetic2DPoint object (lon-lat) coordinates
     */
	@NonNull
    public Geodetic2DPoint getGeodetic() {
        return lonLat;
    }

    /**
     * This accessor method returns the Longitude part of the geodetic coordinates for
     * this UPS coordinate.
     *
     * @return Longitude coordinate of the geodetic (lon-lat) point for this UPS position
     */
	@NonNull
    public Longitude getLongitude() {
        return lonLat.getLongitude();
    }

    /**
     * This accessor method returns the Latitude part of the geodetic coordinates for
     * this UPS coordinate.
     *
     * @return Latitude coordinate of the geodetic (lon-lat) point for this UPS position
     */
	@NonNull
    public Latitude getLatitude() {
        return lonLat.getLatitude();
    }

    /**
     * This accessor method returns the easting part of this UPS coordinate.
     *
     * @return number of meters east of polar false easting
     */
    public double getEasting() {
        return easting;
    }

    /**
     * This accessor method returns the northing part of this UPS coordinate.
     *
     * @return number of meters north of polar false easting
     */
    public double getNorthing() {
        return northing;
    }

    /**
     * This accessor method returns the hemisphere character for this UPS coordinate ('N' or 'S').
     *
     * @return hemisphere character ('N' or 'S')
     */
    public char getHemisphere() {
        return hemisphere;
    }

    /**
     * This accessor method returns the polar zone character for this UPS coordinate
     * ('A' or 'B' for South Pole, 'Y' or 'Z' for North Pole).
     *
     * @return polar zone character ('A' or 'B' for South Pole, 'Y' or 'Z' for North Pole)
     */
    public char getPolarZone() {
        return polarZone;
    }

	@Override
    public boolean equals(Object o) {
        if (this == o) return true;
        if (o == null || getClass() != o.getClass()) return false;

        UPS ups = (UPS) o;
		if (hemisphere != ups.hemisphere) return false;
        if (Double.compare(ups.easting, easting) != 0) return false;
        if (Double.compare(ups.northing, northing) != 0) return false;
        if (!lonLat.equals(ups.lonLat)) return false;
        return true;
    }

    @Override
    public int hashCode() {
        // int result;
        long temp;
		int result = (int) hemisphere;
        temp = easting != +0.0d ? Double.doubleToLongBits(easting) : 0L;
        result = 31 * result + (int) (temp ^ (temp >>> 32));
        temp = northing != +0.0d ? Double.doubleToLongBits(northing) : 0L;
        result = 31 * result + (int) (temp ^ (temp >>> 32));
        result = 31 * result + lonLat.hashCode();
        return result;
    }

    /**
     * This method converts the UPS object to a printable string, using the specified number of
     * fractional digits for the easting and northing offsets.
     *
     * @param fractDig number of digits past the decimal point for easting and northing
     * @return String formatted UPS coordinate
     */
    public String toString(int fractDig) {
        String fract = (fractDig > 0) ? "." : "";
        while (fractDig-- > 0) fract += "0";
        DecimalFormat offFmt = new DecimalFormat("0" + fract);
        return (getEllipsoid().getName() + " UPS " + polarZone + " " +
                offFmt.format(easting) + "m E, " + offFmt.format(northing) + "m N");
    }

    /**
     * This method converts the UPS object to a printable string, with easting and northing offsets
     * given in integer meters.
     *
     * @return String formatted UPS coordinate
     */
    public String toString() {
        return toString(0);
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/ISO8601DateTimeInterval.java<|end_filename|>
/************************************************************************************
 * ISO8601DateTimeInterval.java 12/22/10 10:40 AM psilvey
 *
 *  (C) Copyright MITRE Corporation 2010
 *
 * The program is provided "as is" without any warranty express or implied, including
 * the warranty of non-infringement and the implied warranties of merchantability and
 * fitness for a particular purpose.  The Copyright owner will not be liable for any
 * damages suffered by you as a result of using the Program. In no event will the
 * Copyright owner be liable for any special, indirect or consequential damages or
 * lost profits even if the Copyright owner has been advised of the possibility of
 * their occurrence.
 *
 ***********************************************************************************/
package org.opensextant.geodesy;

import java.util.Calendar;
import java.util.Date;
import java.util.GregorianCalendar;

/**
 * The ISO8601DateTimeInterval class is a simple wrapper for a time interval
 * with a start and end time.
 */
public class ISO8601DateTimeInterval extends ISO8601DateTimePoint {
    private static final String INVALID_ORDER = "invalid time order";

    private long endTime;

    /**
     * Constructor using the current local time for both start and end times
     */
    public ISO8601DateTimeInterval() {
        super();
        this.endTime = this.startTime;
    }

    /**
     * Constructor using the specified start and end times such that
	 * the start time is earlier or equal to the end time.
     *
     * @param startTime long start time for interval
     * @param endTime   long end time for interval
     * @throws IllegalArgumentException if points are not in correct time order
     */
    public ISO8601DateTimeInterval(long startTime, long endTime)
            throws IllegalArgumentException {
        if (endTime < startTime) throw new IllegalArgumentException(INVALID_ORDER);
        this.startTime = startTime;
        this.endTime = endTime;
    }

    /**
     * Constructor that takes an ISO Date Time String time point or interval.
     *
     * @param isoDateTimeStr String containing ISO 8601 Date Time point or interval
     * @throws IllegalArgumentException if a parsing error occurs
     */
    public ISO8601DateTimeInterval(String isoDateTimeStr) {
        try {
            int i = isoDateTimeStr.indexOf("--");
            if (i > 0) {
                ISO8601DateTimeInterval tSeg1 =
                        new ISO8601DateTimeInterval(isoDateTimeStr.substring(0, i));
                ISO8601DateTimeInterval tSeg2 =
                        new ISO8601DateTimeInterval(isoDateTimeStr.substring(i + 2));
                long startTime1 = tSeg1.getStartTimeInMillis();
                long startTime2 = tSeg2.getStartTimeInMillis();
                long endTime1 = tSeg1.getEndTimeInMillis();
                long endTime2 = tSeg2.getEndTimeInMillis();
                if (startTime1 > startTime2 || endTime1 > endTime2)
                    throw new IllegalArgumentException(INVALID_ORDER  + " for " + isoDateTimeStr);
                this.startTime = startTime1;
                this.endTime = endTime2;
            } else {
                String toParse = isoDateTimeStr;
                int eoy = toParse.indexOf("-");
                if (eoy < 0) eoy = toParse.length();
                int n = toParse.length();
                toParse = toParse.replace(" ", "T");
                toParse += dtSuffix.substring(n - eoy);
                Integer calField = null;
                if (n == eoy) calField = Calendar.YEAR;
                else if (n == eoy + 3) calField = Calendar.MONTH;
                else if (n == eoy + 6) calField = Calendar.DAY_OF_MONTH;
                else if (n == eoy + 9) calField = Calendar.HOUR_OF_DAY;
                else if (n == eoy + 12) calField = Calendar.MINUTE;
                else if (n == eoy + 15) calField = Calendar.SECOND;
                GregorianCalendar cal = new GregorianCalendar(UTC_TIMEZONE);
                Date d1 = DF.parse(toParse);
                cal.setTime(d1);
                if (calField == null) {
                    this.startTime = d1.getTime();
                    this.endTime = this.startTime;
                } else {
                    cal.add(calField, 1);
                    Date d2 = cal.getTime();
                    this.startTime = d1.getTime();
                    this.endTime = d2.getTime() - 1;
                }
                // Final validation that yyyy-MM-dd is a valid day (round trip test)
                String ymdInput = toParse.substring(0, toParse.indexOf("T"));
                String ymdOutput = this.toString();
                ymdOutput = ymdOutput.substring(0, ymdOutput.indexOf("T"));
                if (!ymdInput.equals(ymdOutput))
					throw new IllegalArgumentException("Invalid ISO 8601 date and time, " +
                    	isoDateTimeStr);
			}
		} catch (IllegalArgumentException ex) {
			throw ex;
        } catch (Exception ex) {
            throw new IllegalArgumentException("Invalid ISO 8601 date and time, " +
                    isoDateTimeStr, ex);
        }
    }

    /**
     * Getter method for this intervals's end time in milliseconds
     *
     * @return long end time
     */
    public long getEndTimeInMillis() {
        return this.endTime;
    }

    /**
     * Setter method for this interval's end time in milliseconds
     *
     * @param endTime long end time
     */
    public void setEndTimeInMillis(long endTime) {
        this.endTime = endTime;
    }

    @Override
    public boolean equals(Object o) {
        if (this == o) return true;
        if (o == null || getClass() != o.getClass()) return false;
        if (!super.equals(o)) return false;

        ISO8601DateTimeInterval that = (ISO8601DateTimeInterval) o;
        return this.endTime == that.endTime;
    }

    @Override
    public int hashCode() {
        int result = super.hashCode();
        result = 31 * result + (int) (this.endTime ^ (this.endTime >>> 32));
        return result;
    }

    /**
     * Return int indicating whether this ISO8601DateTimeInterval is equal to (0),
     * starts before or starts equal but ends before (-1), or starts after or starts
     * equal and ends after (+1) the specified ISO8601DateTimePoint argument.
     *
     * @param anotherTime ISO8601DateTimePoint to compare to this ISO8601DateTimeInterval
     * @return -1, 0, or +1 depending if this interval is <, ==, or > <tt>anotherTime</tt>
	 * @throws NullPointerException if <tt>anotherTime</tt> is <tt>null</tt>
     */
    public int compareTo(ISO8601DateTimePoint anotherTime) {
        int result = (this.startTime == anotherTime.startTime) ?
                0 : (this.startTime < anotherTime.startTime) ? -1 : +1;
        if (result == 0) {
            long et = (anotherTime instanceof ISO8601DateTimeInterval) ?
                    ((ISO8601DateTimeInterval) anotherTime).endTime : anotherTime.startTime;
            result = (this.endTime == et) ?
                    0 : (this.endTime < et) ? -1 : +1;
        }
        return result;
    }

    /**
     * Format this ISO8601DateTimeInterval object using the ISO 8601 interval syntax,
     * which can be parsed using the String argument constructor method of this
     * class.
     *
     * @return String IOS 8601 Date Time interval syntax for this time segment
     */
    @Override
    public String toString() {
        return DF.format(new Date(this.startTime)) +
                "--" + DF.format(new Date(this.endTime));
    }
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestGeoPoint.java<|end_filename|>
/****************************************************************************************
 *  TestGeoPoint.java
 *
 *  Created: Mar 27, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantibility and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy.test;

import junit.framework.TestCase;
import junit.framework.TestSuite;
import junit.textui.TestRunner;

import org.opensextant.geodesy.*;

import java.util.Random;

public class TestGeoPoint extends TestCase {

    private final Random r = new Random();

    /**
     * This method is used to generate a random Geodetic2DPoint.
     *
     * @param r a Random object to use for selecting the coordinate components
     * @return a random Geodetic2DPoint suitable for testing conversions with
     */
    public static Geodetic2DPoint randomGeodetic2DPoint(Random r) {
        // Define a random Geodetic point
        double lonDeg = (r.nextDouble() * 360.0) - 180.0;  // lonDeg
        Longitude lon = new Longitude(lonDeg, Angle.DEGREES);
        double latDeg = (r.nextDouble() * 180.0) - 90.0;  // latDeg
        Latitude lat = new Latitude(latDeg, Angle.DEGREES);
        return new Geodetic2DPoint(lon, lat);
    }

    public static Geodetic3DPoint randomGeodetic3DPoint(Random r) {
        Geodetic2DPoint pt = randomGeodetic2DPoint(r);
        double elev = r.nextInt(1000000);  // stress test at 1,000 km.
        return new Geodetic3DPoint(pt.getLongitude(), pt.getLatitude(), elev);
    }

    /**
     * This method performs pair-wise format conversions to test for consistency
     * among the GeoPoint representations, using the default FrameOfReference.
     */
    public void testConversions() {
        FrameOfReference f = new FrameOfReference();
        for (int i = 0; i < 1000; i++) {
            GeoPoint a1 = randomGeodetic2DPoint(r);
            GeoPoint a2 = f.toTopocentric(a1);
            GeoPoint a3 = f.toGeocentric(a1);
            GeoPoint a4 = f.toGeodetic(a2);
            assertTrue(f.proximallyEquals(a1, a2));
            assertTrue(f.proximallyEquals(a1, a3));
            assertTrue(f.proximallyEquals(a1, a4));
        }
    }

    public void testReference() {
        Geodetic3DPoint p = randomGeodetic3DPoint(r);
        FrameOfReference f = new FrameOfReference(p);
        assertEquals(p, f.getTopographicOrigin());

        final Ellipsoid ellip = Ellipsoid.getInstance("WGS 84");
        f = new FrameOfReference(ellip);
        assertEquals(ellip, f.getEllipsoid());
    }

	public void testBadReference() {
		Ellipsoid ellip = null;
		try {
			new FrameOfReference(ellip);
			fail("Expected NullPointerException");
		} catch(NullPointerException npe) {
			// expected
		}
		try {
			new FrameOfReference(ellip, randomGeodetic3DPoint(r));
			fail("Expected NullPointerException");
		} catch(NullPointerException npe) {
			// expected
		}
	}

    /**
     * This method performs format conversions with string constructors in degrees and dms
     */
    public void testRandomEquals() {
        for (int i = 0; i < 1000; i++) {
            Geodetic2DPoint a1 = randomGeodetic2DPoint(r);
            Geodetic2DPoint a2 = new Geodetic2DPoint(a1.getLongitude().inDegrees() + "," + a1.getLatitude().inDegrees());
            // note by making equals() work in round off errors (such as in this case) using Angle.equals() vs phi1=phi2 && lamb1==lamb2
            // but break contract in hashCode such that a.equals(b) -> true but hashCode(a) may not equal hashCode(b)
            assertEquals(a1, a2);
            assertEquals(a1.hashCode(), a2.hashCode());

            String aStr = a1.toString(7);
            Geodetic2DPoint a3 = new Geodetic2DPoint(aStr);

            /*
            System.out.println(a1.hashCode() + " " + a2.hashCode() + " " + a1.equals(a2)
                    + "\n\t" + aStr
                    + "\n\t" + a3.toString(6)
                    + "\n\t" + a1.getLongitude().inRadians() + " " + a1.getLatitude().inRadians()
                    + "\n\t" + a2.getLongitude().inRadians() + " " + a2.getLatitude().inRadians()
                    + "\n\t" + a3.getLongitude().inRadians() + " " + a3.getLatitude().inRadians());
            */

            /*
               Following test uncovered a bug near poles with decimal angle = 0 in
               south or west resuling in angle deg = -0 being used in floating point
               calculations and sign being incorrectly assigned.

                   (0 deg 46' 41.23530" W, 6 deg 18' 6.01751" N) parsed with wrong sign on longitude:
                   (0 deg 46' 41.23530" E, 6 deg 18' 6.01751" N)
                   a1 lon-lat (rads) = -0.013580771984344651 0.10998491659692418
                   a3 lon-lat (rads) =  0.013580771974469753 0.1099849165873849

                   dealing with -0 had problems with sign test: fixed in Angle(String valStr) constructor.

                   Note: 6 decimal places in seconds should be enough decimals such that
                   round off error adjustment in equals() is true in nearly all cases.

               Here's cases that fails if only using 5 decimal places for seconds where Angle.EPISON = 1e-7:
                   (60 deg 41' 16.04803" E, 0 deg 0' 12.67506" N)
                   (60 deg 41' 16.04803" E, 0 deg 0' 12.67506" N)
                   a1 lon-lat (rads) = 1.0592017708190238 6.145040119951688E-5
                   a3 lon-lat (rads) = 1.0592017707968808 6.145042496884236E-5
                   a1 != a3

                   If change EPSILON in Angle equals(double, double) from current value of 1e-7 then
                   will need to adjust # fractional digits of precision used above.
                   */
            assertEquals(a1, a3);
        }
    }

    private final Geodetic2DPoint a = new Geodetic2DPoint(new Longitude(-1), new Latitude(1));
    private final Geodetic2DPoint b = new Geodetic2DPoint(new Longitude(-1), new Latitude(1));
    private final Geodetic2DPoint c = new Geodetic2DPoint(new Longitude(-1.001), new Latitude(1));

	public void testLatLonAngle() {
		Latitude lat = new Latitude(new Angle());
		Longitude lon = new Longitude(new Angle());
		Geodetic2DPoint a = new Geodetic2DPoint(lon, lat);
		Geodetic2DPoint b = new Geodetic2DPoint( new Longitude(new Angle(0, Angle.DEGREES)),
				new Latitude(new Angle(0, Angle.DEGREES)));
		assertEquals(a, b);
	}

    public void testEquals() {
        // test equality with known geo-points
        assertEquals(a, b);
        assertFalse(a.equals(c));

        Geodetic2DPoint a2 = new Geodetic2DPoint(new Longitude(a.getLongitude().inRadians()),
                new Latitude(a.getLatitude().inRadians()));
        assertEquals(a, a2);

        // test equality with any possible round off errors
        Geodetic2DPoint a3 = new Geodetic2DPoint(new Longitude(Math.toRadians(a.getLongitude().inDegrees())),
                new Latitude(Math.toRadians(a.getLatitude().inDegrees())));

        assertEquals(a, a3);
    }

    public void testNullCompare() {
        Geodetic2DPoint other = null;
        assertFalse(a.equals(other));
    }

    public void testHashCode() {
        assertEquals(a.hashCode(), b.hashCode());
        assertTrue(a.hashCode() != c.hashCode()); // this test isn't required by equals-hashCode contract but by how the hashCode is computed
    }

    public void testStringCreation() {
        String[] latLons = {
                // "37d25'19.07\"N, 122d05'06.24\"W", // this format does not work
                "(122\u00B0 5' 6\" W, 37\u00B0 25' 19\" N)",
                "37 25 19.07 N, 122 05 06.24 W", // latitude first
                "37 25'19.07\"N, 122 05'06.24\"W", // no whitespace between min and sec fields
                "(42\u00B0 22' 11.77\" N, 71\u00B0 1' 40.30\" W)", // latitude first
                "(51 deg 28' 15.19\" N, 45 deg 27' 33.41\" W)",
                "51 deg 28' 15.19\" N, 45 deg 27' 33.41\" W",
                "(12 34 56E, 45 34 23N)",
                "12 34 56E, 45 34 23N"
        };
        for (String latLon : latLons) {
            // System.out.println(latLon);
            try {
                new Geodetic2DPoint(latLon);
            } catch (IllegalArgumentException e) {
                System.out.println("Failed on " + latLon);
                throw e;
            }
            //System.out.println("\t\t\t\t" + pt);
        }
    }

    public void testToString() {
        Geodetic2DPoint p = randomGeodetic2DPoint(r);
        int prevLen = 0;

        String base = p.toString();
        String prefix = base.substring(0, base.indexOf(' '));   // e.g. "75 deg"
        String suffix = base.substring(base.length() - 2);      // e.g. "N)"

        /*
        length of string should increase two digits for each fractional # of digits in lat & lon
        (75 deg 10' 31" E, 89 deg 17' 7" N)
        (75 deg 10' 31.5" E, 89 deg 17' 7.2" N)
        (75 deg 10' 31.46" E, 89 deg 17' 7.18" N)
        (75 deg 10' 31.455" E, 89 deg 17' 7.177" N)
        (75 deg 10' 31.4551" E, 89 deg 17' 7.1768" N)
        (75 deg 10' 31.45508" E, 89 deg 17' 7.17678" N)
         ...
         */
        for (int i = 0; i < 6; i++) {
            String s = p.toString(i);
            System.out.println(s);
            int len = s.length();
            assertTrue(len >= prevLen + 2);
            // simple tests: 1) each should start with same degrees and 2) end with same hemisphere letter
            assertTrue(s.startsWith(prefix));
            assertTrue(s.endsWith(suffix));
            prevLen = len;
        }
    }

    /**
     * Main method for running class tests.
     *
     * @param args standard command line arguments - ignored.
     */
    public static void main(String[] args) {
        TestSuite suite = new TestSuite(TestGeoPoint.class);
        new TestRunner().doRun(suite);
    }
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestAngle.java<|end_filename|>
/***************************************************************************
 * $Id$
 *
 * (C) Copyright MITRE Corporation 2006-2008
 *
 * The program is provided "as is" without any warranty express or implied,
 * including the warranty of non-infringement and the implied warranties of
 * merchantability and fitness for a particular purpose.  The Copyright
 * owner will not be liable for any damages suffered by you as a result of
 * using the Program.  In no event will the Copyright owner be liable for
 * any special, indirect or consequential damages or lost profits even if
 * the Copyright owner has been advised of the possibility of their
 * occurrence.
 *
 ***************************************************************************/
package org.opensextant.geodesy.test;

import java.io.*;
import java.util.Random;

import junit.framework.JUnit4TestAdapter;

import org.junit.Assert;
import org.junit.Test;
import org.opensextant.geodesy.Angle;

/**
 * 
 * @author <NAME>
 */
public class TestAngle {

    /**
     * This is the unit test case for functionality provided by the Angle class
     */
	@Test
    public void testAngle() {
        Random r = new Random();
        int fractDig = 2;
        for (int i = 0; i < 1000; i++) {
            // Test radians constructor and conversion to degrees
            double rad = (r.nextDouble() * (2.0 * Math.PI)) - Math.PI;
            Angle r1 = new Angle(rad, Angle.RADIANS);
            Angle r2 = new Angle(r1.inDegrees(), Angle.DEGREES);
            Assert.assertEquals(r1.toString(fractDig), r2.toString(fractDig));

            // Test degrees, minutes, seconds constructor and conversion to radians
            int degrees = r.nextInt(360) - 180;
            int minutes = r.nextInt(60);
            double seconds = r.nextDouble() * 60.0;
            Angle d1 = new Angle(degrees, minutes, seconds);
            Angle d2 = new Angle(d1.inRadians());
			Assert.assertEquals(d1.toString(fractDig), d2.toString(fractDig));

            // Test String constructor
            Angle s1 = new Angle(r1.toString(5));
            Angle s2 = new Angle(s1.inDegrees(), Angle.DEGREES);
			Assert.assertEquals(s1.toString(fractDig), s2.toString(fractDig));
        }

        Angle z000 = new Angle();   // default constructor makes 0 deg angle
        // Some common positive angles (testing add method)
        Angle p045 = new Angle(45.0, Angle.DEGREES);
        Angle p090 = p045.add(p045);
        Angle p135 = p090.add(p045);
        Angle p180 = p090.add(p090);
        // Some common negative angles (testing difference method)
        Angle m045 = p045.difference(z000);
        Angle m090 = p090.difference(z000);
        Angle m135 = p090.difference(m045);
        Angle m180 = p090.difference(m090);

        Assert.assertTrue(p135.add(p090).equals(m135));
        Assert.assertEquals(m180.toString(fractDig), p180.toString(fractDig));
        Assert.assertEquals(m180, p180);

//        System.out.println("\rString constructor parse tests:");
//        String[] parseTests = new String[]{
//                "77.5 degrees", "77 deg:15 min", "1.4r", "-23 deg, 15 min", "3.141592654 radians",
//                "-180deg, 15:23.2", "1:2:5", "4 23 15", "60 61", "25:25.3:56"
//        };
//        for (String s : parseTests) {
//            System.out.print("Angle string \"" + s + "\" interpreted as ");
//            try {
//                System.out.println(new Angle(s).toString(fractDig));
//            } catch (Exception ex) {
//                System.out.println(ex);
//            }
//        }
    }

	@Test
	public void testEquals() {
		Random r = new Random();
        for (int i = 0; i < 1000; i++) {
            // Test radians constructor and conversion to degrees with equals() and hashCode() test
            double rad = (r.nextDouble() * (2.0 * Math.PI)) - Math.PI;
            Angle r1 = new Angle(rad, Angle.RADIANS);
			Angle r2 = new Angle(r1.inDegrees(), Angle.DEGREES);
			Assert.assertEquals(r1, r2);
			Assert.assertEquals(r1.hashCode(), r2.hashCode());
			// now test from degrees back to radians
			Angle r3 = new Angle(r2.inRadians(), Angle.RADIANS);
			Assert.assertEquals(r2, r3);
			Assert.assertEquals(r2.hashCode(), r3.hashCode());
		}
		// small EPSILON delta but equal
		Angle r1 = new Angle(1);
		Angle r2 = new Angle(r1.inRadians() + 1e-8);
		Assert.assertEquals(r1, r2);
		Assert.assertEquals(r1.hashCode(), r2.hashCode());
		// "large" enough EPSILON delta so not equal
		Angle r3 = new Angle(r1.inRadians() + 1e-7);
		Assert.assertFalse(r1.equals(r3));
	}

	@Test
	public void testCompare() {
		Angle a = new Angle(0);
		Angle b = new Angle(0, Angle.DEGREES);
		Assert.assertEquals(0, a.compareTo(b));
		double delta = 0.1;
		for(int i=0; i < 8; i++) {
			Angle c = new Angle(a.inRadians() + delta);
			Assert.assertEquals(-1, a.compareTo(c));
			delta /= 10.0;
		}
		Assert.assertEquals(0.0, a.getValue(Angle.RADIANS), 1e-6);
		Assert.assertEquals(0.0, a.getValue(Angle.DEGREES), 1e-6);
	}

	@Test
    public void testNullAngleCompare() {
		Angle r1 = new Angle(Math.PI, Angle.RADIANS);
		Angle r2 = null;
		Assert.assertFalse(r1.equals(r2));
	}

	@Test
	public void testDegreeSymbol() throws IOException {
		InputStreamReader reader = null;
		final InputStream stream = getStream("DEGREE_SYMBOL.txt");
		Assert.assertNotNull(stream);
		try {
			reader = new InputStreamReader(stream, "ISO-8859-1");
			final StringBuilder buff = new StringBuilder();
			for (int c = reader.read(); c != -1; c = reader.read()) {
				buff.append((char) c);
			}
			final String file_degree = buff.toString().trim();
			Assert.assertEquals(Angle.DEGSYM, file_degree);
		} finally {
			if (reader != null)
				try {
					reader.close();
				} catch (IOException e) { }
				try {
					stream.close();
				} catch (IOException e) { }
		}
	}

    private InputStream getStream(String filename) throws IOException {
        File file = new File("src/test/resources/org/opensextant/geodesy/test/" + filename);
        if (file.exists()) return new java.io.FileInputStream(file);
        return getClass().getResourceAsStream(filename);
    }

	public static junit.framework.Test suite() {
		return new JUnit4TestAdapter(TestAngle.class);
	}

	/**
     * Main method for running class tests.
     *
     * @param args standard command line arguments - ignored.
     */
    public static void main(String[] args) {
		new junit.textui.TestRunner().doRun(suite());
		System.exit(0);
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Topocentric3DPoint.java<|end_filename|>
/****************************************************************************************
 *  Topocentric3DPoint.java
 *
 *  Created: Mar 14, 2007
 *
 *  @author <NAME>, based on work by <NAME>, <NAME>, and <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import java.text.DecimalFormat;

/**
 * The Topocentric3DPoint class extends the simple Topocentric2DPoint class by adding an
 * explicit elevation value, in meters. To be correctly interpreted, one needs a particular
 * reference Topocentric origin point and Ellipsoid model of the earth (see FrameOfReference
 * class for ways to define and use this in coordinate conversions). The coordinate system
 * used here is the one defined by Nato Standards for Airborne Ground Survelliance (AGS)
 * data. In particular, the x-axis is positive in the East direction, the y-axis is
 * positive in the North direction, and the z-axis is positive above the tangent plane that
 * touches the surface of the Ellipsoid at the topocentric origin point. All units of length
 * are specified in meters. Note that there are other systems of topocentric reference that
 * make different assumptions about the directional positioning of the x, y, and z axis.
 * This class is lightweight and simple in the interest of efficient storage and processing.
 */
public class Topocentric3DPoint extends Topocentric2DPoint {
	
	private static final long serialVersionUID = 1L;
	
    private double elevation;

    /**
     * The constructor takes double precision floating point values for easting and
     * northing distances in meters from the topographic origin, along with an elevation
     * value in meters.  The origin itself is a contextual variable set in the
     * FrameOfReference object and used there by the coordinate transform methods that
     * need it.
     *
     * @param easting   number of meters along x-axis from topographic origin
     * @param northing  number of meters along y-axis from topographic origin
     * @param elevation number of meters along z-axis coordinate from topographic origin
     */
    public Topocentric3DPoint(double easting, double northing, double elevation) {
        super(easting, northing);
        this.elevation = elevation;
    }

    /**
     * This accessor method is used to get the elevation (assumed in meters).
     *
     * @return elevation of this Topocentric3DPoint point (in meters).
     */
    public double getElevation() {
        return elevation;
    }

    /**
     * This settor method is used to update the elevation value for this Topocentric3DPoint point.
     *
     * @param elevation height of this Topocentric3DPoint point above reference plane (in meters).
     */
    public void setElevation(double elevation) {
        this.elevation = elevation;
    }

    /**
     * This method returns a hash code for this Topocentric3DPoint object. The
     * result is the exclusive OR of the component values to maintain the general
     * contract for the hashCode method, which states that equal objects must have
     * equal hash codes.
     *
     * @return a hash code value for this object.
     */
    public int hashCode() {
        return 31 * super.hashCode() +
                Double.valueOf(elevation).hashCode();
    }

    /**
     * This method is used to test whether two points are equal in the sense that
     * have the same coordinate value.
     *
     * @param that Topocentric3DPoint point to compare against this one.
     * @return true if specified Topocentric3DPoint point is equal in value to this
     *         Topocentric3DPoint point.
     */
    public boolean equals(Topocentric3DPoint that) {
        return that != null && this.easting == that.easting &&
				this.northing == that.northing &&
				this.elevation == that.elevation;
    }

    /**
     * This method is used to test whether two points are equal in the sense that
     * have the same coordinate value.
     *
     * @param that Topocentric3DPoint point to compare against this one.
     * @return true if specified Topocentric3DPoint point is equal in value to this
     *         Topocentric3DPoint point.
     */
    public boolean equals(Object that) {
        return that instanceof Topocentric3DPoint && equals((Topocentric3DPoint) that);
    }

    /**
     * This method abstracts this Topocentric3DPoint object as a general purpose GeoPoint.
     *
     * @param fRef the FrameOfReference in which to interpret this coordinate.
     * @return the equivalent Geodetic3DPoint
     */
    public Geodetic3DPoint toGeodetic3D(FrameOfReference fRef) {
        return fRef.toGeodetic(this);
    }

    /**
     * This toString method returns a string formatted version of this Topocentric3DPoint,
     * printed with the specified number of decimal digits for each component.
     *
     * @param fractDig number of digits to the right of the decimal point to be formatted.
     * @return string version of this Topocentric3DPoint, suitable for debugging purposes.
     */
    public String toString(int fractDig) {
        String result = super.toString(fractDig);
        String fract = (fractDig > 0) ? "." : "";
        while (fractDig-- > 0) fract += "0";
        DecimalFormat fmt = new DecimalFormat("0" + fract);
        return (result + " @ " + fmt.format(elevation) + "m");
    }

    /**
     * This toString method returns a string formatted version of this Topocentric3DPoint,
     * printed with integer values for each component.
     *
     * @return string version of this Topocentric3DPoint, suitable for debugging purposes.
     */
    public String toString() {
        return toString(0);
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/GeocentricPoint.java<|end_filename|>
/****************************************************************************************
 *  GeocentricPoint.java
 *
 *  Created: Mar 13, 2007
 *
 *  @author <NAME>, based on work by <NAME>, <NAME>, and <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import java.io.Serializable;
import java.text.DecimalFormat;

/**
 * The GeocentricPoint class is used to represent a Earth-Centered, Earth-Fixed (ECEF)
 * point in 3D Cartesian coordinate space. The origin of the 3-space is the center of
 * the earth, which coincides with Ellipsoid model center points like WGS-84. The x-axis
 * lies on the Equitorial Plane and passes through the Prime Meridian in the positive
 * direction.  The y-axis is likewise in the Equitorial Plane, 90&deg East of the x-axis.
 * The z-axis passes through the North Pole in the positive direction, making this
 * reference system a right-hand-system.
 */
public class GeocentricPoint implements GeoPoint, Serializable {
	private static final long serialVersionUID = 1L;
	
	private double x;
    private double y;
    private double z;

    /**
     * The constructor takes the 3 coordinate points, measured in meters from the center
     * of the earth, ordered x, y, and z.
     *
     * @param x double distance in meters along the geocentric x-axis
     * @param y double distance in meters along the geocentric y-axis
     * @param z double distance in meters along the geocentric z-axis
     */
    public GeocentricPoint(double x, double y, double z) {
        this.x = x;
        this.y = y;
        this.z = z;
    }

    /**
     * This accessor method returns the x coordinate (in meters) for this GeocentricPoint.
     *
     * @return double distance in meters along the geocentric x-axis
     */
    public double getX() {
        return x;
    }

    /**
     * This settor method updates the x coordinate (in meters) for this GeocentricPoint.
     *
     * @param x double distance in meters along the geocentric x-axis
     */
    public void setX(double x) {
        this.x = x;
    }

    /**
     * This accessor method returns the y coordinate (in meters) for this GeocentricPoint.
     *
     * @return double distance in meters along the geocentric y-axis
     */
    public double getY() {
        return y;
    }

    /**
     * This settor method updates the y coordinate (in meters) for this GeocentricPoint.
     *
     * @param y double distance in meters along the geocentric y-axis
     */
    public void setY(double y) {
        this.y = y;
    }

    /**
     * This accessor method returns the z coordinate (in meters) for this GeocentricPoint.
     *
     * @return double distance in meters along the geocentric z-axis
     */
    public double getZ() {
        return z;
    }

    /**
     * This settor method updates the z coordinate (in meters) for this GeocentricPoint.
     *
     * @param z double distance in meters along the geocentric z-axis
     */
    public void setZ(double z) {
        this.z = z;
    }

    /**
     * This method returns a hash code for this GeocentricPoint object. The
     * result is the exclusive OR of the component values to maintain the general
     * contract for the hashCode method, which states that equal objects must have
     * equal hash codes.
     *
     * @return a hash code value for this object.
     */
    public int hashCode() {
        return 31 * Double.valueOf(x).hashCode() +
                Double.valueOf(y).hashCode() ^
                Double.valueOf(z).hashCode();
    }

    /**
     * This method is used to test whether two points are equal in the sense that
     * have the same coordinate value.
     *
     * @param that GeocentricPoint point to compare against this one.
     * @return true if specified GeocentricPoint point is equal in value to this
     *         GeocentricPoint point.
     */
    public boolean equals(GeocentricPoint that) {
        return that != null && ((this.x == that.x) && (this.y == that.y) && (this.z == that.z));
    }

    /**
     * This method is used to test whether two points are equal in the sense that
     * have the same coordinate value.
     *
     * @param that GeocentricPoint point to compare against this one.
     * @return true if specified GeocentricPoint point is equal in value to this
     *         GeocentricPoint point.
     */
    public boolean equals(Object that) {
        return that instanceof GeocentricPoint && equals((GeocentricPoint) that);
    }

    /**
     * This method abstracts this Geocentric object as a general purpose GeoPoint.
     *
     * @param fRef the FrameOfReference in which to interpret this coordinate.
     * @return the equivalent Geodetic3DPoint
     */
    public Geodetic3DPoint toGeodetic3D(FrameOfReference fRef) {
        return fRef.toGeodetic(this);
    }

    /**
     * This toString method returns a string formatted version of this GeocentricPoint,
     * as '(x, y, z)', where each component is printed with the specified number of decimal
     * digits.
     *
     * @param fractDig number of digits to the right of the decimal point to be formatted.
     * @return string version of this GeocentricPoint, suitable for debugging purposes.
     */
    public String toString(int fractDig) {
        String fract = (fractDig > 0) ? "." : "";
        while (fractDig-- > 0) fract += "0";
        DecimalFormat fmt = new DecimalFormat("0" + fract);
        return ("(" + fmt.format(x) + "m, " + fmt.format(y) + "m, " + fmt.format(z) + "m)");
    }

    /**
     * The toString method returns a string formatted version of this GeocentricPoint,
     * as '(x, y, z)', where each component is rounded to integer meters.
     *
     * @return string version of this GeocentricPoint, suitable for debugging purposes.
     */
    public String toString() {
        return toString(0);
    }

}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Geodetic3DBounds.java<|end_filename|>
/****************************************************************************************
 *  Geodetic3DBounds.java
 *
 *  Created: Apr 3, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

/**
 * The Geodetic3DBounds class extends the simple Geodetic2DBounds class by adding explicit
 * min and max elevation values, in meters.
 */
public class Geodetic3DBounds extends Geodetic2DBounds {
	
	private static final long serialVersionUID = 1L;

    // using 1e-3 meters for elevation equality gives a precision up to 1 millimeter
    private static final double DELTA = 1e-3;

    public double minElev;
    public double maxElev;

    /**
     * The default constructor returns a bounding box containing only the 0-0-0 point.
     */
    public Geodetic3DBounds() {
        //minElev = 0.0;
        //maxElev = 0.0;
    }

    /**
     * This constructor will accept west and east Geodetic3DPoint coordinates and set the
     * bounding box limits accordingly.  Note that the source coordinate components are
     * used by reference (i.e. they are not cloned).  The first coordinate is used to set
     * the west Longitude even if its value is greater than the east, so boxes that span
     * the international date line can be accommodated by this convention of listing lons
     * in reverse order.  On the other hand, latitude values are normalized so that south
     * is always less than north. Min and Max elevation are set accordingly.
     *
     * @param westCoordinate western most point
     * @param eastCoordinate eastern most point
     */
    public Geodetic3DBounds(Geodetic3DPoint westCoordinate, Geodetic3DPoint eastCoordinate) {
        super(westCoordinate, eastCoordinate);
        // Normalize elevation values so min < max
        if (westCoordinate.getElevation() < eastCoordinate.getElevation()) {
            minElev = westCoordinate.getElevation();
            maxElev = eastCoordinate.getElevation();
        } else {
            minElev = eastCoordinate.getElevation();
            maxElev = westCoordinate.getElevation();
        }
    }

    /**
     * This constructor takes a single Geodetic3DPoint seed point and constructs a degenerate
     * bounding box that contains it.  It is primarily used with the include method to grow a
     * bounding box around a set of points.
     *
     * @param seedPoint Initial point to use as limits of bounding box.
     */
    public Geodetic3DBounds(Geodetic3DPoint seedPoint) {
        super(seedPoint);
        minElev = seedPoint.getElevation();
        maxElev = minElev;
    }

    /**
     * This constructor takes a single Geodetic3DBounds as a seed and constructs a new
     * bounding box by cloning its values. It is primarily used with the include methods
     * to grow a bounding box around a set of points or other bounding boxes.
     *
     * @param seedBBox Initial bounding box to use as limits of bounding box.
     */
    public Geodetic3DBounds(Geodetic3DBounds seedBBox) {
        super(seedBBox);
        minElev = seedBBox.minElev;
        maxElev = seedBBox.maxElev;
    }

    /**
     * This constructor takes a {@code Geodetic2DBounds} as a seed and constructs a new
     * bounding box by cloning its values in addition to the min/max elevation bounds.
     *
     * @param seedBBox Initial bounding box to use as limits of bounding box
     * @param minElev min elevation (in meters)
     * @param maxElev max elevation (in meters)
     */
    public Geodetic3DBounds(Geodetic2DBounds seedBBox, double minElev, double maxElev) {
        super(seedBBox);
        if (minElev > maxElev) {
            // min/max elevations in reverse order. swap values.
            this.maxElev = minElev;
            this.minElev = maxElev;
        } else {
            this.minElev = minElev;
            this.maxElev = maxElev;
        }
    }

    /**
     * This method is used to extend this bounding box to include a new point. If the new
     * point is already inside the existing bounding box, no change will result. The bounding
     * box is always extended by the smallest amount necessary to include the new point.
     * Sometimes this will cause the bounding box to wrap around the international date line.
     *
     * @param newPoint new Geodetic3DPoint point to include in this bounding box.
     */
    public void include(Geodetic3DPoint newPoint) {
        super.include(newPoint);
        double elev = newPoint.getElevation();
        if (elev < minElev) minElev = elev;
        else if (elev > maxElev) maxElev = elev;
    }

    /**
     * This method is used to extend this bounding box to include a new point. If the new
     * point is already inside the existing bounding box, no change will result. The bounding
     * box is always extended by the smallest amount necessary to include the new point.
     * Sometimes this will cause the bounding box to wrap around the international date line.
     *
     * @param newPoint new Geodetic2DPoint point to include in this bounding box.
     */
    public void include(Geodetic2DPoint newPoint) {
        if (newPoint instanceof Geodetic3DPoint)
            include((Geodetic3DPoint)newPoint);
        else
            super.include(newPoint);
    }

    /**
     * This method is used to extend this bounding box to include another bounding box. If the
     * newly specified bounding box is inside the existing bounding box, no change will result.
     * This bounding box is always extended by the smallest amount necessary to include the new
     * bounding box. Sometimes this will cause this bounding box to wrap around the international
     * date line.
     *
     * @param bbox additional Geodetic3DBounds bounding box to include in this bounding box.
     */
    public void include(Geodetic3DBounds bbox) {
        this.include(new Geodetic3DPoint(bbox.getWestLon(), bbox.getSouthLat(), bbox.minElev));
        this.include(new Geodetic3DPoint(bbox.getEastLon(), bbox.getNorthLat(), bbox.maxElev));
    }

    /**
     * This method is used to extend this bounding box to include another bounding box. If the
     * newly specified bounding box is inside the existing bounding box, no change will result.
     * This bounding box is always extended by the smallest amount necessary to include the new
     * bounding box. Sometimes this will cause this bounding box to wrap around the international
     * date line.
     *
     * @param bbox additional Geodetic2DBounds bounding box to include in this bounding box.
     */
	public void include(Geodetic2DBounds bbox) {
		if (bbox instanceof Geodetic3DBounds)
            include((Geodetic3DBounds)bbox);
        else
            super.include(bbox);
	}

    /**
     * This predicate method determines whether the specified geodetic point is
     * contained within this bounding box.
     *
     * @param testPoint Geodetic3DPoint to test for containment within this bounding box
     * @return true if specified point is within this Geodetic3DBounds, false otherwise.
     */
    public boolean contains(Geodetic3DPoint testPoint) {
        double testElev = testPoint.getElevation();
        return (super.contains(testPoint) && (minElev <= testElev) && (testElev <= maxElev));
    }

    /**
     * This predicate method determines whether the specified geodetic point is
     * contained within this bounding box.
     *
     * @param testPoint Geodetic2DPoint to test for containment within this bounding box
     * @return true if specified point is within this Geodetic3DBounds, false otherwise.
     */
    public boolean contains(Geodetic2DPoint testPoint) {
        if (testPoint instanceof Geodetic3DPoint)
            return contains((Geodetic3DPoint)testPoint);
        else
            return super.contains(testPoint);
    }

    /**
     * This predicate method determines whether the specified geodetic bounding box is
     * contained within this bounding box.
     *
     * @param testBox Geodetic3DBounds to test for containment within this bounding box
     * @return true if the specified box is contained within this bounding box
     */
    public boolean contains(Geodetic3DBounds testBox) {
        double testMinElev = testBox.minElev;
        double testMaxElev = testBox.maxElev;
        return (super.contains(testBox) &&
                (minElev <= testMinElev) && (testMinElev <= maxElev) &&
                (minElev <= testMaxElev) && (testMaxElev <= maxElev));
    }

    /**
     * This predicate method determines whether the specified geodetic bounding box is
     * contained within this bounding box.
     *
     * @param testBox Geodetic2DBounds to test for containment within this bounding box
     * @return true if the specified box is contained within this bounding box
     */
    public boolean contains(Geodetic2DBounds testBox) {
        if (testBox instanceof Geodetic3DBounds)
            return contains((Geodetic3DBounds)testBox);
        else
            return super.contains(testBox);
    }

    /**
     * This method is used to determine the geodetic point that lies at the center of
     * this bounding box.
     *
     * @return Geodetic3DPoint that lies at the center of this Geodetic2DBounds box
     */
    @NonNull
    public Geodetic3DPoint getCenter() {
        double westLonRad = this.getWestLon().inRadians;
        double eastLonRad = this.getEastLon().inRadians;
        // If longitudes wrap, adjust east to be greater than west before calc
        if (westLonRad > eastLonRad) eastLonRad += Angle.TWO_PI;
        double centLonRad = westLonRad + ((eastLonRad - westLonRad) / 2.0);
        double southLatRad = this.getSouthLat().inRadians;
        double northLatRad = this.getNorthLat().inRadians;
        double centLatRad = southLatRad + ((northLatRad - southLatRad) / 2.0);
        double centElev = minElev + ((maxElev - minElev) / 2.0);
        // Note that Longitude constructor will re-normalize angle if > 360 deg
        return new Geodetic3DPoint(new Longitude(centLonRad), new Latitude(centLatRad), centElev);
    }

    /**
     * Returns a hash code for this <code>Geodetic3DBounds</code> object. The
     * result is the exclusive OR of the latitude, longitude, and elevation values to maintain
     * the general contract for the hashCode method, which states
     * that equal objects must have equal hash codes.
     *
     * @return a <code>hash code</code> value for this object.
     */
    public int hashCode() {
        return 31 * super.hashCode() + Double.valueOf(maxElev).hashCode() ^ Double.valueOf(minElev).hashCode();
    }

    /**
     * Compares this object against the specified object.  The result
     * is <code>true</code> if and only if the argument is not
     * <code>null</code> and is a <code>Geodetic3DBounds</code> object that
     * has the same latitude, longitude, and elevation values.
     *
     * @param that the object to compare with.
     * @return <code>true</code> if the objects are the same;
     *         <code>false</code> otherwise.
     */
    public boolean equals(Object that) {
		if (that instanceof Geodetic3DBounds)
			return equals((Geodetic3DBounds) that);
		return that instanceof Geodetic2DBounds && equals((Geodetic2DBounds) that);
	}

    /**
     * The equals method tests for bounding box coordinate numeric equality
     *
     * @param that Geodetic3DBounds object to compare to this one
     * @return true if specified Geodetic3DBounds is spatially equivalent to this one
     */
    public boolean equals(Geodetic3DBounds that) {
        return this == that || that != null && equals(that, that.minElev, that.maxElev);
    }

	/**
     * The equals method tests for bounding box coordinate numeric equality
     *
     * @param that Geodetic3DBounds object to compare to this one
     * @return true if specified Geodetic3DBounds is spatially equivalent to this one
     */
    public boolean equals(Geodetic2DBounds that) {
	if (that instanceof Geodetic3DBounds)
		return equals((Geodetic3DBounds) that);
	return equals(that, 0, 0);
    }

    private boolean equals(Geodetic2DBounds that, double minElev, double maxElev) {
	return super.equals(that) && Math.abs(this.maxElev - maxElev) < DELTA &&
			Math.abs(this.minElev - minElev) < DELTA;
    }

    /**
     * The toString method formats the bounding box coordinates for printing,
     * as SW .. NE points, followed by min and max elevation in meters.
     *
     * @return String representation of corner points of bounding box plus elevation range.
     */
    public String toString() {
        String result = super.toString();
        return (result + " .. (" + minElev + "m, " + maxElev + "m)");
    }
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestEllipsoid.java<|end_filename|>
/****************************************************************************************
 *  TestEllipsoid.java
 *
 *  Created: Mar 27, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantibility and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy.test;

import junit.framework.TestSuite;
import junit.framework.TestCase;
import junit.textui.TestRunner;

import org.opensextant.geodesy.Ellipsoid;
import org.opensextant.geodesy.Geodetic2DPoint;

import java.util.Set;

public class TestEllipsoid extends TestCase {

    /**
     * This method exercises some of the Ellipsoid API, including performing
     * some surface distance calculation tests between airports of major cities.
     */
    public void testDistances() {
        Ellipsoid ellip = null;
        Set<String> names = Ellipsoid.getEllipsoidNames();
        for (String name : names) {
            ellip = Ellipsoid.getInstance(name);
            if (name.equals("WGS 84")) break;
        }
        assertNotNull(ellip);
        assertEquals(ellip, Ellipsoid.getInstance("WGS 84"));

        double distInM;
        long distInKm;

        // Do some distance calculations
        // Note: for Angle constructor parsing, degree symbol can be encoded in either
        // ISO-8859-1 (Latin-1) or unicode, but toString formatting always uses unicode
        Geodetic2DPoint London =
                new Geodetic2DPoint("51 deg 28' 15.19\" N, 0 deg 27' 33.41\" W");
        Geodetic2DPoint Boston =
                new Geodetic2DPoint("42\u00B0 22' 11.77\" N, 71\u00B0 1' 40.30\" W");
        Geodetic2DPoint Honolulu =
                new Geodetic2DPoint("21 deg 19' 55.36\" N, 157 deg 55' 11.08\" W");
        Geodetic2DPoint Sydney =
                new Geodetic2DPoint("33\u00B0 56' 43.70\" S, 151\u00B0 10' 51.43\" E");

        distInM = ellip.orthodromicDistance(London, Boston);
        assertEquals((long) distInM, (long) ellip.orthodromicDistance(Boston, London)); // round off to nearest meter
        distInKm = Math.round(distInM / 1000.0);
        assertEquals(5256L, distInKm);

        distInM = ellip.orthodromicDistance(Honolulu, Sydney);
        assertEquals(distInM, ellip.orthodromicDistance(Sydney, Honolulu));
        distInKm = Math.round(distInM / 1000.0);
        assertEquals(8155L, distInKm);

        distInM = ellip.orthodromicDistance(Honolulu, London);
        assertEquals(distInM, ellip.orthodromicDistance(London, Honolulu));
        distInKm = Math.round(distInM / 1000.0);
        assertEquals(11646L, distInKm);
    }

    /**
     * Main method for running class tests.
     *
     * @param args standard command line arguments - ignored.
     */
    public static void main(String[] args) {
        TestSuite suite = new TestSuite(TestEllipsoid.class);
        new TestRunner().doRun(suite);
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Geodetic3DPoint.java<|end_filename|>
/****************************************************************************************
 *  Geodetic3DPoint.java
 *
 *  Created: Mar 12, 2007
 *
 *  @author <NAME>, based on work by <NAME>, <NAME>, and <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.text.DecimalFormat;

/**
 * The Geodetic3DPoint class extends the simple Geodetic2DPoint class by adding an explicit
 * elevation value, in meters. To be correctly interpreted, one needs a particular Ellipsoid
 * model of the earth (see FrameOfReference class for ways to define and use this in coordinate
 * conversions). This class is simple in the interest of efficient storage and processing.
 * A Geodetic3DPoint can be converted to a GeocentricPoint and back with loss of precision
 * to within a tenth of an arc second, for any points on the surface of the earth and any
 * ocean depths, as well as up to altitudes that include those achievable by "air-breather"
 * aircraft (tested up to 26km).
 */
public class Geodetic3DPoint extends Geodetic2DPoint {

	private static final long serialVersionUID = 1L;

    private double elevation;

    /**
     * This constructor takes a Longitude, Latitude, and an elevation value in meters.
     * The default elevation reference point (vertical datum) is the tangent surface plane
     * touching the WGS-84 Ellipsoid at the given longitude and latitude.  The elevation
     * distance is measured along the normal vector through the surface plane. Note that for
     * Geodetic representations on an Ellipsoid earth model, this vector does not generally
     * pass through the center of the Ellipsoid.
     *
     * @param lon       Longitude of this Geodetic3DPoint.
     * @param lat       Latitude of this Geodetic3DPoint.
     * @param elevation elevation in meters from the assumed reference point
	 * @throws NullPointerException if lon or lat are null
     */
    public Geodetic3DPoint(Longitude lon, Latitude lat, double elevation) {
        super(lon, lat);
        this.elevation = elevation;
    }

    /**
     * Default constructor makes a geodetic point at the central meridian on the equator (0, 0, 0).
     */
    public Geodetic3DPoint() {
        super();
    }

    /**
     * This accessor method is used to get the elevation (assumed in meters).
     *
     * @return elevation of this Geodetic3DPoint point (in meters).
     */
    public double getElevation() {
        return elevation;
    }

    /**
     * This settor method is used to update the elevation value for this Geodetic3DPoint.
     *
     * @param elevation elevation of this Geodetic3DPoint point (in meters).
     */
    public void setElevation(double elevation) {
        this.elevation = elevation;
    }

    // using 1e-3 meters for elevation equality gives a precision up to 1 millimeter
    private static final double DELTA = 1e-3;

    /**
     * This method returns a hash code for this Geodetic3DPoint object. The
     * result is the exclusive OR of the component values to maintain the general
     * contract for the hashCode method, which states that equal objects must have
     * equal hash codes.
     *
     * @return a hash code value for this object.
     */
    public int hashCode() {
        // Note we're using approximate equals vs absolute equals on floating point numbers
        // so we ignore beyond ~3 decimal places in computing the hashCode, otherwise
        // we break the equals-hashCode contract. Changing EPSILON or equals(Geodetic3DPoint)
        // may require changing the logic used here also.
        return 31 * lon.hashCode() + lat.hashCode() ^ ((int) (elevation / DELTA));
    }

   /**
     * This method is used to test whether two points are equal in the sense that have
     * the same angular coordinate values and elevations, to within epsilon. See also
     * {@code proximallyEquals} method in the {@link FrameOfReference} class for
     * additional equality tests.
     *
     * @param that Geodetic3DPoint point to compare against this one.
     * @return <tt>true</tt> if specified Geodetic3DPoint point is equal in value to this
     *			Geodetic3DPoint. If point is a {@link Geodetic2DPoint} then 0 elevation
     * 			is used for comparison.
     */
    public boolean equals(Object that) {
		if (this == that) {
	    	return true;
		}
        if (that instanceof Geodetic3DPoint)
			return this.eq((Geodetic3DPoint) that);
		if (that instanceof Geodetic2DPoint) {
			Geodetic2DPoint pt = (Geodetic2DPoint)that;
			return eq(pt.lon, pt.lat, 0);  // test with 0 elevation
		}
		return false;
    }

	// Inherited Javadoc
    private boolean eq(Geodetic3DPoint that) {
        // angular diff on surface Earth quite different than Earth to nearest star
        return that != null && eq(that.lon, that.lat, that.elevation);
    }

	// Inherited Javadoc
	private boolean eq(Longitude otherLon, Latitude otherLat, double otherElevation) {
		return Angle.equals(this.lon.inRadians, otherLon.inRadians) &&
                Angle.equals(this.lat.inRadians, otherLat.inRadians) &&
                (Double.compare(this.elevation, otherElevation) == 0 ||
                        Math.abs(this.elevation - otherElevation) <= DELTA);
	}

    /**
     * This method abstracts this Geodetic3DPoint object as a general purpose GeoPoint.
     *
     * @param fRef the FrameOfReference in which to interpret this coordinate.
     * @return the equivalent Geodetic3DPoint
     */
	@NonNull
    public Geodetic3DPoint toGeodetic3D(FrameOfReference fRef) {
        return this;
    }

    /**
     * This toString method returns a string formatted version of this Geodetic3DPoint,
     * as '(lon, lat) @ h', lon and lat are printed with the specified number of decimal
     * digits of arc seconds, and h is printed as fractional meters with the number of
     * digits of precision specified.
     *
     * @param fractDig number of digits to the right of the decimal point to be formatted.
     * @return string version of this Geodetic3DPoint, suitable for debugging purposes.
     */
	@NonNull
    public String toString(int fractDig) {
        String result = super.toString(fractDig);
        String fract = (fractDig > 0) ? "." : "";
        while (fractDig-- > 0) fract += "0";
        DecimalFormat fmt = new DecimalFormat("0" + fract);
        return result + " @ " + fmt.format(elevation) + "m";
    }

    /**
     * The toString method returns a string formatted version of this Geodetic3DPoint,
     * as '(lon, lat) @ h', where lon and lat are printed in integer degrees, minutes,
     * and seconds, and elevation h is rounded to integer meters.
     *
     * @return string version of this Geodetic3DPoint, suitable for debugging purposes.
     */
	@NonNull
    public String toString() {
        return toString(0);
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/ISO8601DateTimePoint.java<|end_filename|>
/************************************************************************************
 * ISO8601DateTimePoint.java 12/22/10 10:16 AM psilvey
 *
 *  (C) Copyright MITRE Corporation 2010
 *
 * The program is provided "as is" without any warranty express or implied, including
 * the warranty of non-infringement and the implied warranties of merchantability and
 * fitness for a particular purpose.  The Copyright owner will not be liable for any
 * damages suffered by you as a result of using the Program. In no event will the
 * Copyright owner be liable for any special, indirect or consequential damages or
 * lost profits even if the Copyright owner has been advised of the possibility of
 * their occurrence.
 *
 ***********************************************************************************/
package org.opensextant.geodesy;

import java.util.Date;
import java.util.TimeZone;

/**
 * The ISO8601DateTimePoint class is a simple wrapper for a time instant standardized
 * as a number of milliseconds since midnight Jan 1, 1970 GMT. This class has a
 * String constructor and toString method that respectively parse and format dates
 * and times in the ISO 8601 format. Partial dates are parsed to the beginning of the
 * time interval they define.
 */
public class ISO8601DateTimePoint implements Comparable<ISO8601DateTimePoint> {
    protected static final String INVALID_POINT = "invalid time point";

    protected final static TimeZone UTC_TIMEZONE = TimeZone.getTimeZone("UTC");
    protected final static String dtSuffix = "-01-01T00:00:00.000Z";
    protected final static SafeDateFormat DF;

    static {
        DF = new SafeDateFormat("yyyy-MM-dd'T'HH:mm:ss.SSS'Z'");
        DF.setTimeZone(UTC_TIMEZONE);
    }

    protected long startTime;

    /**
     * Constructor using the current local time for start time
     */
    public ISO8601DateTimePoint() {
        this.startTime = System.currentTimeMillis();
    }

    /**
     * Constructor using the specified start time
     *
     * @param startTime long start time for interval
     */
    public ISO8601DateTimePoint(long startTime) {
        this.startTime = startTime;
    }

    /**
     * Constructor that takes an ISO Date Time String to be interpreted as an
     * instant in time (start of interval if partial date time specified).
     *
     * @param isoDateTimeStr String containing an ISO 8601 Date Time point or interval
     * @throws IllegalArgumentException if a parsing error occurs
     */
    public ISO8601DateTimePoint(String isoDateTimeStr) {
		String ymdInput, ymdOutput;
        try {
            String toParse = isoDateTimeStr;
            int eoy = toParse.indexOf("-");
            if (eoy < 0) eoy = toParse.length();
            int n = toParse.length();
            toParse = toParse.replace(" ", "T");
            toParse += dtSuffix.substring(n - eoy);
            this.startTime = DF.parse(toParse).getTime();
            // Final validation that yyyy-MM-dd is a valid day (round trip test)
            ymdInput = toParse.substring(0, toParse.indexOf("T"));
            ymdOutput = this.toString();
            ymdOutput = ymdOutput.substring(0, ymdOutput.indexOf("T"));
        } catch (Exception ex) {
            throw new IllegalArgumentException("Invalid ISO 8601 date and time, " +
                    isoDateTimeStr, ex);
        }
		if (!ymdInput.equals(ymdOutput))
				throw new IllegalArgumentException("Invalid ISO 8601 date and time, " +
                    isoDateTimeStr);
    }

    /**
     * Getter method for this point's start time in milliseconds
     *
     * @return long start time
     */
    public long getStartTimeInMillis() {
        return this.startTime;
    }

    /**
     * Setter method for this point's start time in milliseconds
     *
     * @param startTime long start time
     */
    public void setStartTimeInMillis(long startTime) {
        this.startTime = startTime;
    }

    @Override
    public boolean equals(Object o) {
        if (this == o) return true;
        if (o == null || getClass() != o.getClass()) return false;

        ISO8601DateTimePoint that = (ISO8601DateTimePoint) o;
        return this.startTime == that.startTime;
    }

    @Override
    public int hashCode() {
        return (int) (this.startTime ^ (this.startTime >>> 32));
    }

    /**
     * Return int indicating whether this ISO8601DateTimePoint is equal to (0),
     * is before (-1), or is after (+1) the specified ISO8601DateTimePoint argument.
     *
     * @param anotherTime ISO8601DateTimePoint to compare to this ISO8601DateTimePoint
     * @return -1, 0, or +1 depending if this point is <, ==, or > <tt>anotherTime</tt>
	 * @throws NullPointerException if <tt>anotherTime</tt> is <tt>null</tt>
     */
    public int compareTo(ISO8601DateTimePoint anotherTime) {
        return (this.startTime == anotherTime.startTime) ?
                0 : (this.startTime < anotherTime.startTime) ? -1 : + 1;
    }

    /**
     * Format this TmPoint object using the ISO 8601 point syntax, which can
     * be parsed using the String argument constructor method of this class.
     *
     * @return String ISO 8601 Date Time point syntax for this TmPoint
     */
    @Override
    public String toString() {
        return DF.format(new Date(this.startTime));
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/FrameOfReference.java<|end_filename|>
/****************************************************************************************
 *  FrameOfReference.java
 *
 *  Created: Mar 14, 2007
 *
 *  @author <NAME>, based on work by <NAME>, <NAME>, and <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.io.Serializable;

/**
 * The FrameOfReference class is a contextual class to allow our geo-spatial point
 * objects to remain lightweight and context free by themselves, yet allow them to
 * be accurately converted between formats or compared when necessary. This class
 * holds the global and local reference points and mathematical modeling assumptions
 * for converting geo-spatial points from one representational system to another.
 * In particular, it holds the Ellipsoid model of the earth and the Topocentric point
 * of origin, along with methods to convert points between Geodetic, Geocentric, and
 * Topocentric point formats, given this FrameOfReference. The default Ellipsoid is
 * WGS-84, and the default Topocentric point of origin is where the equator meets the
 * prime meridian (0 deg E, 0 deg N). All of the point types are objects, and there are
 * three different kinds of equivalence tests that may be applied to them. As objects
 * with unique handles, they will only be '=' if they map to the same address in
 * memory (i.e. are the exact same object in the Java Virtual Machine). If they have
 * the same values for their components but different addresses in memory, they are
 * considered 'equal' (each individual point class has an 'equals' method to test this).
 * This is the sense of equality that is used for their hashcode identity. Finally,
 * since coordinate conversion is subject to some loss of precision and round-off error,
 * the FrameOfReference class provides a method (called proximallyEquals) to compare
 * points for near-equality in space. As currently coded, points are considered equal
 * if they are within 1 meter of each other on the surface of the Ellipsoid, and also
 * within 1 meter of each other in elevation value from the surface.  As a result, this
 * equality test requires that non-Geodetic representations be first converted to their
 * Geodetic form, and surface distance is measured by the distance method of the
 * Ellipsoid class (which is a more accurate form of surface distance than Greater
 * Circle distances on a Sphere).
 */
public class FrameOfReference implements Serializable {

	private static final long serialVersionUID = 1L;
	
	@NonNull private Ellipsoid ellip;
	@NonNull private Geodetic3DPoint topoOrigLle;
	@NonNull private GeocentricPoint topoOrigEcf;

    private static final Ellipsoid DEFAULT_ELLIPSOID =
            Ellipsoid.getInstance("WGS 84");
    private static final Geodetic3DPoint DEFAULT_TOPOORIG =
            new Geodetic3DPoint(new Longitude(0.0), new Latitude(0.0), 0.0);

    /*
     * This method converts a Geodetic3DPoint to a GeocentricPoint (lle 2 ecf).
     *
     * @param lle the Geodetic3DPoint (Longitude, Latitude, Elevation or lle)
     * @return the GeocentricPoint (Earth Centered Fixed or ecf)
     * @throws NullPointerException if ellip or lle are null
     */
    private GeocentricPoint toGeocentric(Geodetic3DPoint lle) {
        double lambda = lle.getLongitude().inRadians();
        Latitude lat = lle.getLatitude();
        double phi = lat.inRadians();
        double h = lle.getElevation();

        double N = ellip.getRadius(lat);

        double cosPhi = Math.cos(phi);
        double X = (N + h) * cosPhi * Math.cos(lambda);
        double Y = (N + h) * cosPhi * Math.sin(lambda);
        double Z = ((N * (1.0 - ellip.getEccentricitySquared())) + h) * Math.sin(phi);

        return new GeocentricPoint(X, Y, Z);
    }

    /*
     * This method converts a GeocentricPoint to a TopocentricPoint (ecf 2 tcs).
     *
     * @param ecf the GeocentricPoint (Earth Centered Fixed or ecf)
     * @return the TopocentricPoint (Topo Centric System or tcs)
     * @throws NullPointerException if ecf is null
     */
    private Topocentric3DPoint toTopocentric(GeocentricPoint ecf) {
        double x0 = ecf.getX() - topoOrigEcf.getX();
        double y0 = ecf.getY() - topoOrigEcf.getY();
        double z0 = ecf.getZ() - topoOrigEcf.getZ();

        double lonRad = topoOrigLle.getLongitude().inRadians();
        double cosOrLon = Math.cos(lonRad);
        double sinOrLon = Math.sin(lonRad);

        double latRad = topoOrigLle.getLatitude().inRadians();
        double cosOrLat = Math.cos(latRad);
        double sinOrLat = Math.sin(latRad);

        double x = cosOrLon * x0 + sinOrLon * y0;
        double tcsX = -sinOrLon * x0 + cosOrLon * y0;
        double tcsY = -sinOrLat * x + cosOrLat * z0;
        double tcsZ = cosOrLat * x + sinOrLat * z0;
        return new Topocentric3DPoint(tcsX, tcsY, tcsZ);
    }

    /*
     * This method converts a Geodetic3DPoint to a TopocentricPoint (lle 2 ecf 2 tcs).
     * It does this via an intermediate conversion using a GeocentricPoint (ecf).
     *
     * @param lle the Geodetic3DPoint (Longitude, Latitude, Elevation or lle)
     * @return the TopocentricPoint (Topo Centric System or tcs)
     */
    private Topocentric3DPoint toTopocentric(Geodetic3DPoint lle) {
        GeocentricPoint ecf = this.toGeocentric(lle);
        return this.toTopocentric(ecf);
    }

    /*
     * This method converts a TopocentricPoint to a GeocentricPoint (tcs 2 ecf).
     *
     * @param tcs the TopocentricPoint (Topo Centric System or tcs)
     * @return the GeocentricPoint (Earth Centered Fixed or ecf)
     */
    private GeocentricPoint toGeocentric(Topocentric3DPoint tcs) {
        double lonRad = topoOrigLle.getLongitude().inRadians();
        double sinLambda = Math.sin(lonRad);
        double cosLambda = Math.cos(lonRad);

        double latRad = topoOrigLle.getLatitude().inRadians();
        double sinPhi = Math.sin(latRad);
        double cosPhi = Math.cos(latRad);

        double x = (-1.0 * sinPhi * tcs.getNorthing() + cosPhi * tcs.getElevation());
        double y = tcs.getEasting();

        double ecfX = cosLambda * x - sinLambda * y + topoOrigEcf.getX();
        double ecfY = sinLambda * x + cosLambda * y + topoOrigEcf.getY();
        double ecfZ = (cosPhi * tcs.getNorthing() +
                sinPhi * tcs.getElevation() + topoOrigEcf.getZ());
        return new GeocentricPoint(ecfX, ecfY, ecfZ);
    }

    /*
     * This method converts a GeocentricPoint to a Geodetic3DPoint (ecf 2 lle).
     *
     * @param ecf the GeocentricPoint (Earth Centered Fixed or ecf)
     * @return the Geodetic3DPoint (Longitude, Latitude, Elevation or lle)
     */
    private Geodetic3DPoint toGeodetic(GeocentricPoint ecf) {
        double X = ecf.getX();
        double Y = ecf.getY();
        double Z = ecf.getZ();
        double a = ellip.getEquatorialRadius();
        double b = ellip.getPolarRadius();

        double p = Math.sqrt((X * X) + (Y * Y));
        Longitude lon;
        Latitude lat;
        double h;
        if (p == 0.0) {
            // At a Pole all Longitude values are at the same place, so we might normalize
            // lon to be zero, but we're using the Topocentric Origin's lon instead
            lon = topoOrigLle.getLongitude();
            lat = new Latitude((Z > 0.0) ? +90.0 : -90.0, Angle.DEGREES);
            h = (Z > 0.0) ? Z - b : b - Z;
        } else {
            lon = new Longitude(Math.atan2(Y, X));
            double theta = Math.atan((Z * a) / (p * b));
            double sinTheta = Math.sin(theta);
            double cosTheta = Math.cos(theta);
            double phi = Math.atan(
                    (Z + (ellip.get2ndEccentricitySquared() * b * sinTheta * sinTheta * sinTheta)) /
                            (p - (ellip.getEccentricitySquared() * a * cosTheta * cosTheta * cosTheta)));
            lat = new Latitude(phi);
            h = (p / Math.cos(phi)) - ellip.getRadius(lat);
        }
        return new Geodetic3DPoint(lon, lat, h);
    }

    /*
     * This method converts a TopocentricPoint to a Geodetic3DPoint (tcs 2 lle).
     * It does this via an intermediate conversion using a GeocentricPoint (ecf).
     *
     * @param tcs the TopocentricPoint (Topo Centric System or tcs)
     * @return the Geodetic3DPoint (Longitude, Latitude, Elevation or lle)
     */
    private Geodetic3DPoint toGeodetic(Topocentric3DPoint tcs) {
        GeocentricPoint ecf = this.toGeocentric(tcs);
        return this.toGeodetic(ecf);
    }

    /**
     * This constructor takes an Ellipsoid object earth model and a Geodetic3DPoint
     * to use as the Topographic Origin point for coordinate conversions.
     *
     * @param ellip    the Ellipsoid model of the earth to use in conversions
     * @param topoOrig the topographic origin point as a Geodetic3DPoint object.
	 * @throws NullPointerException if ellip or topoOrig are null
     */
    public FrameOfReference(Ellipsoid ellip, Geodetic3DPoint topoOrig) {
        this.ellip = ellip;
        this.topoOrigLle = topoOrig;
        this.topoOrigEcf = this.toGeocentric(topoOrig);
    }

    /**
     * This constructor takes an Ellipsoid object earth model and defaults the
     * Topographic Origin to be (0&deg N, 0&deg E, 0m).
     *
     * @param ellip the Ellipsoid model of the earth to use in conversions.
	 * @throws NullPointerException if ellip is null
     */
    public FrameOfReference(Ellipsoid ellip) {
        this(ellip, DEFAULT_TOPOORIG);
    }

    /**
     * This constructor takes a Geodetic3DPoint to use as the Topographic Origin
     * point for coordinate conversions, and defaults the Ellipsoid to be WGS-84.
     *
     * @param topoOrig the topographic origin point as a Geodetic3DPoint object.
	 * @throws NullPointerException if topoOrig is null
	 */
    public FrameOfReference(Geodetic3DPoint topoOrig) {
        this(DEFAULT_ELLIPSOID, topoOrig);
    }

    /**
     * This constructor uses the default WGS-84 Ellipsoid and the default
     * Topographic origin of (0&deg N, 0&deg E, 0m) for conversions.
     */
    public FrameOfReference() {
        this(DEFAULT_ELLIPSOID, DEFAULT_TOPOORIG);
    }

    /**
     * This accessor method is used to get the current Ellipsoid model of the
     * earth for this FrameOfReference.
     *
     * @return Ellipsoid object being used.
     */
    @NonNull
	public Ellipsoid getEllipsoid() {
        return ellip;
    }

    /**
     * This settor method is used to update the current Ellipsoid model of the
     * earth for this FrameOfReference.
     *
     * @param ellip the new Ellipsoid object to be used.
	 * @throws NullPointerException if ellip is null
     */
    public void setEllipsoid(Ellipsoid ellip) {
		if (ellip == null) throw new NullPointerException();
        this.ellip = ellip;
    }

    /**
     * This accessor method is used to get the current Topographic Origin point
     * (as a Geodetic3DPoint object) for this FrameOfReference.
     *
     * @return the topographic origin point as a Geodetic3DPoint object.
     */
    public Geodetic3DPoint getTopographicOrigin() {
        return topoOrigLle;
    }

    /**
     * This settor method is used to update the current Topographic Origin point
     * (as a Geodetic3DPoint object) for this FrameOfReference.
     *
     * @param topoOrig the new topographic origin point as a Geodetic3DPoint object.
	 * @throws NullPointerException if topoOrig is null
	 */
    public void setTopographicOrigin(Geodetic3DPoint topoOrig) {
        this.topoOrigLle = topoOrig;
        this.topoOrigEcf = this.toGeocentric(topoOrig);
    }

    /**
     * This method calculates the orthodromic distance in meters between two GeoPoints, which
     * is the shortest distance between two geodetic points on the surface of the ellipsoid
     * used by this FrameOfReference. Note that any elevation values are ignored as the
     * points are projected onto the ellipsoid at their geodetic (lon-lat) coordinates.
     *
     * @param p1 the starting GeoPoint.
     * @param p2 the ending GeoPoint.
     * @return the approximate distance in surface meters between the two specified points.
     */
    public double orthodromicDistance(GeoPoint p1, GeoPoint p2) {
        return ((p1 instanceof Geodetic2DPoint) && (p2 instanceof Geodetic2DPoint)) ?
                ellip.orthodromicDistance((Geodetic2DPoint) p1, (Geodetic2DPoint) p2) :
                ellip.orthodromicDistance(this.toGeodetic(p1), this.toGeodetic(p2));
    }

    /**
     * This method computes the shortest straight line distance (Euclidean Distance)
     * through space between the specified GeoPoints. Note that it may pass through the
     * earth.
     *
     * @param p1 the starting GeoPoint.
     * @param p2 the ending GeoPoint.
     * @return the shortest straight line distance between the two specified points.
     */
    public double euclideanDistance(GeoPoint p1, GeoPoint p2) {
        GeocentricPoint g1 = this.toGeocentric(p1);
        GeocentricPoint g2 = this.toGeocentric(p2);
        double xDiff = g1.getX() - g2.getX();
        double yDiff = g1.getY() - g2.getY();
        double zDiff = g1.getZ() - g2.getZ();
        return (Math.sqrt((xDiff * xDiff) + (yDiff * yDiff) + (zDiff * zDiff)));
    }

    /**
     * This predicate method is used to determine if two GeoPoints are close enough
     * in this FrameOfReference to be considered 'equal'.  The criteria is that they
     * are within a meter of each other on the surface of the Ellipsoid, and also
     * within a meter of each other in elevation.
     *
     * @param p1 the first GeoPoint to be compared.
     * @param p2 the second GeoPoint to be compared.
     * @return true if the two points are within a meter of each other on the surface and
     *         also in elevation.
     */
    public boolean proximallyEquals(GeoPoint p1, GeoPoint p2) {
        boolean closeEnough;
        if ((p1 instanceof Geodetic2DPoint) && (p2 instanceof Geodetic2DPoint)) {
            double dist = ellip.orthodromicDistance((Geodetic2DPoint) p1, (Geodetic2DPoint) p2);
            closeEnough = (dist < 1.0);
        } else {
            Geodetic3DPoint g1 = this.toGeodetic(p1);
            Geodetic3DPoint g2 = this.toGeodetic(p2);
            closeEnough = ((this.orthodromicDistance(g1, g2) < 1.0) &&
                    (Math.abs(g1.getElevation() - g2.getElevation()) < 1.0));
        }
        return closeEnough;
    }

    /**
     * This method is the public conversion method for using this FrameOfReference
     * to convert any GeoPoint object into a Geodetic3DPoint.
     *
     * @param gp the GeoPoint to convert.
     * @return the equivalent Geodetic3DPoint.
     * @throws IllegalArgumentException error if the GeoPoint instance is not recognized.
     */
    public Geodetic3DPoint toGeodetic(GeoPoint gp) {
        if (gp instanceof Geodetic3DPoint)
            return (Geodetic3DPoint) gp;
        else if ((gp instanceof Geodetic2DPoint) || (gp instanceof MGRS))
            return gp.toGeodetic3D(this);
        else if (gp instanceof GeocentricPoint)
            return this.toGeodetic((GeocentricPoint) gp);
        else if (gp instanceof Topocentric3DPoint)
            return this.toGeodetic((Topocentric3DPoint) gp);
        else if (gp instanceof Topocentric2DPoint)
            return this.toGeodetic(((Topocentric2DPoint) gp).toTopocentric3D());
        else throw new IllegalArgumentException("unknown GeoPoint instance");
    }

    /**
     * This method is the public conversion method for using this FrameOfReference
     * to convert any GeoPoint object into a GeocentricPoint.
     *
     * @param gp the GeoPoint to convert.
     * @return the equivalent GeocentricPoint.
     * @throws IllegalArgumentException error if the GeoPoint instance is not recognized.
     */
    public GeocentricPoint toGeocentric(GeoPoint gp) {
        if (gp instanceof GeocentricPoint)
            return (GeocentricPoint) gp;
        else if (gp instanceof Geodetic3DPoint)
            return this.toGeocentric((Geodetic3DPoint) gp);
        else if ((gp instanceof Geodetic2DPoint) || (gp instanceof MGRS))
            return this.toGeocentric(gp.toGeodetic3D(this));
        else if (gp instanceof Topocentric3DPoint)
            return this.toGeocentric((Topocentric3DPoint) gp);
        else if (gp instanceof Topocentric2DPoint)
            return this.toGeocentric(((Topocentric2DPoint) gp).toTopocentric3D());
        else throw new IllegalArgumentException("unknown GeoPoint instance");
    }

    /**
     * This method is the public conversion method for using this FrameOfReference
     * to convert any GeoPoint object into a Topocentric3DPoint.
     *
     * @param gp the GeoPoint to convert.
     * @return the equivalent Topocentric3DPoint.
     * @throws IllegalArgumentException error if the GeoPoint instance is not recognized.
     */
    public Topocentric3DPoint toTopocentric(GeoPoint gp) {
        if (gp instanceof Topocentric3DPoint)
            return (Topocentric3DPoint) gp;
        else if (gp instanceof Topocentric2DPoint)
            return ((Topocentric2DPoint) gp).toTopocentric3D();
        else if (gp instanceof Geodetic3DPoint)
            return this.toTopocentric((Geodetic3DPoint) gp);
        else if ((gp instanceof Geodetic2DPoint) || (gp instanceof MGRS))
            return this.toTopocentric(gp.toGeodetic3D(this));
        else if (gp instanceof GeocentricPoint)
            return this.toTopocentric((GeocentricPoint) gp);
        else throw new IllegalArgumentException("unknown GeoPoint instance");
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Geodetic2DEllipse.java<|end_filename|>
/****************************************************************************************
 *  Geodetic2DEllipse.java
 *
 *  Created: Nov 2, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.io.Serializable;
import java.lang.ref.WeakReference;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

/**
 * The Geodetic2DEllipse class represents an ellipse on the surface of the earth
 * (itself modeled as an Ellipsoid). Similar to the Geodetic2DPoint class, instances
 * of this object are lightweight in the sense that they represent ellipses in terms
 * of simple parameters: center, semiMajorAxis, semiMinorAxis, and orientation.
 * The center point is specified as a Geodetic2DPoint, the axes are given in meters,
 * and the orientation is an Angle (clockwise positive from North to the semiMajor
 * axis). Just as with the Geodetic2DPoint class, the reference Ellipsoid earth model
 * is not explicitly carried with the objects, but may be needed for some uses of the
 * ellipse. Because an ellipse is symmetrical along either axis, two ellipses whose
 * orientations differ by 180 degrees are equal if all the other parameters are equal.
 */
public class Geodetic2DEllipse implements Serializable {
    private static final long serialVersionUID = 1L;

    @NonNull private Geodetic2DPoint center;
    private double semiMajorAxis;
    private double semiMinorAxis;
    @NonNull private Angle orientation;
    private transient WeakReference<List<Geodetic2DPoint>> boundary;
    
    private static final String NULL_ERROR = "Ellipse parameters can not be null";
    private static final String AXIS_ERROR = "Semi-major axis must be greater than semi-minor axis";

    /**
     * Default constructor makes a geodetic ellipse at the central meridian on the
     * equator (0, 0), with 0 semi major and minor axes, and 0 degrees orientation.
     */
    public Geodetic2DEllipse() {
        this.center = new Geodetic2DPoint();
        this.semiMajorAxis = 0.0;
        this.semiMinorAxis = 0.0;
        this.orientation = new Angle();
    }

    /**
     * The constructor takes a center point, semi major and minor axes in meters,
     * and angle or orientation from North to the semi major axis.
     *
     * @param center        Geodetic2DPoint at the center of the ellipse
     * @param semiMajorAxis length in meters of the ellipse's semi major axis
     * @param semiMinorAxis length in meters of the ellipse's semi minor axis
     * @param orientation   Angle from North to semi major axis for this ellipse
     * @throws IllegalArgumentException if null parameter or axis ordering error
     */
    public Geodetic2DEllipse(Geodetic2DPoint center, double semiMajorAxis,
                             double semiMinorAxis, Angle orientation)
            throws IllegalArgumentException {
        if ((center == null) || (orientation == null))
            throw new IllegalArgumentException(NULL_ERROR);
        if (semiMajorAxis < semiMinorAxis)
            throw new IllegalArgumentException(AXIS_ERROR);
        this.center = center;
        this.semiMajorAxis = semiMajorAxis;
        this.semiMinorAxis = semiMinorAxis;
        this.orientation = orientation;
    }

    
    /**
     * Get the boundary list for the ellipse
     * @param count the number of slices that the list should include, non zero
     * and should be a multiple of four for best results.
     * @return the boundary iterator
     */
	@NonNull
    public Iterable<Geodetic2DPoint> boundary(int count) {
    	List<Geodetic2DPoint> blist;
    	if (boundary != null) {
    		blist = boundary.get();
    		if (blist != null && blist.size() == count) {
    			return blist;
    		}
    	}
    	blist = Collections.unmodifiableList(makeBoundary(count));
    	boundary = new WeakReference<List<Geodetic2DPoint>>(blist);
    	return blist;
    }
    
    /**
     * Make a boundary list for a given count
     * @param count the number of slices in the boundary list
     * @return the boundary list, never <code>null</code> or empty
     */
	@NonNull
    private List<Geodetic2DPoint> makeBoundary(int count) {
    	// omega is the ellipse orientation angle (final rotation adjustment)
        Angle omega = new Angle(getOrientation().inRadians());
        // theta is the normalized circular sweep angle that
        // varys from -180 to +180 degrees
        Angle theta = new Angle(-Math.PI);
        // delta is the variable angular increment value
        Angle delta;

        // Use polar coordinate equation for an ellipse, where e is the
        // eccentricity.  Angle (theta) sweeps around the circle, and r
        // changes as a function of it
        double t, r;
        double a = getSemiMajorAxis();
        double b = getSemiMinorAxis();
        double e = Math.sqrt(1.0 - ((b * b) / (a * a)));

        // nominal delta in degrees (average angle for n steps around circle)
        double nDelta = 360.0 / count;
        Angle nominalDelta = new Angle(nDelta, Angle.DEGREES);

        // We compute r as a function of theta,
        // then rotate by ellipse orientation amount
        Geodetic2DArc arc = new Geodetic2DArc(getCenter(), 0.0, omega);
        Angle stepAngle = new Angle(0.0, Angle.DEGREES);
        List<Geodetic2DPoint> rval = new ArrayList<Geodetic2DPoint>();
        for (int step = 0; step < count; step++) {
            // y cycles from -1 to +1 to -1 from 0 deg to 180 deg, repeats
            double y = ((1.0 - Math.cos(stepAngle.inRadians() * 2.0)) - 1.0);
            // delta varies from 0 to 2 * nominal delta
            // to make small steps near major axis and larger ones near minor
            delta = new Angle(nDelta + (nDelta * y), Angle.DEGREES);
            theta = theta.add(delta);
            t = Math.cos(theta.inRadians());
            r = b / (Math.sqrt(1.0 - (e * e * t * t)));
            arc.setDistanceAndAzimuth(r, theta.add(omega));
            rval.add(arc.getPoint2());
            stepAngle = stepAngle.add(nominalDelta);
        }
        return rval;
    }
    
    /**
     * This method returns the Geodetic2DPoint at the center of this ellipse.
     *
     * @return Geodetic2DPoint center of ellipse
     */
	@NonNull
    public Geodetic2DPoint getCenter() {
        return center;
    }

    /**
     * This method is used to set the geodetic center point of this ellipse.
     *
     * @param center new Geodetic2DPoint center of this ellipse
     * @throws IllegalArgumentException error if center is null
     */
    public void setCenter(Geodetic2DPoint center) {
        if (center == null)
            throw new IllegalArgumentException(NULL_ERROR);
        this.center = center;
    }

    /**
     * This method returns the semi major axis of this ellipse, in meters.
     *
     * @return double distance in meters of this ellipse's semi major axis
     */
    public double getSemiMajorAxis() {
        return semiMajorAxis;
    }

    /**
     * This method is used to set the semi major axis (in meters) for this ellipse.
     *
     * @param semiMajorAxis double distance in meters of this ellipse's semi major axis
     * @throws IllegalArgumentException if axis ordering error
     */
    public void setSemiMajorAxis(double semiMajorAxis) {
        if (semiMajorAxis < this.semiMinorAxis)
            throw new IllegalArgumentException(AXIS_ERROR);
        this.semiMajorAxis = semiMajorAxis;
    }

    /**
     * This method returns the semi minor axis of this ellipse, in meters.
     *
     * @return double distance in meters of this ellipse's semi minor axis
     */
    public double getSemiMinorAxis() {
        return semiMinorAxis;
    }

    /**
     * This method is used to set the semi minor axis (in meters) for this ellipse.
     *
     * @param semiMinorAxis double distance in meters of this ellipse's semi minor axis
     * @throws IllegalArgumentException if axis ordering error
     */
    public void setSemiMinorAxis(double semiMinorAxis) {
        if (this.semiMajorAxis < semiMinorAxis)
            throw new IllegalArgumentException(AXIS_ERROR);
        this.semiMinorAxis = semiMinorAxis;
    }

    /**
     * This method is used to set the semi major and minor axes of this ellipse, in meters.
     *
     * @param semiMajorAxis length in meters of the ellipse's semi major axis
     * @param semiMinorAxis length in meters of the ellipse's semi minor axis
     * @throws IllegalArgumentException if axis ordering error
     */
    public void setSemiAxes(double semiMajorAxis, double semiMinorAxis)
            throws IllegalArgumentException {
        if (semiMajorAxis < semiMinorAxis)
            throw new IllegalArgumentException(AXIS_ERROR);
        this.semiMajorAxis = semiMajorAxis;
        this.semiMinorAxis = semiMinorAxis;
    }

    /**
     * This method returns the Angle of orientation from North to the semi major axis
     * of this ellipse.
     *
     * @return Angle of orientation of this ellipse
     */
	@NonNull
    public Angle getOrientation() {
        return orientation;
    }

    /**
     * This method is used to set the Angle of orientation from North to the semi major
     * axis for this ellipse.
     *
     * @param orientation Angle from North to the semi major axis for this ellipse
     * @throws IllegalArgumentException error if orientation is null
     */
    public void setOrientation(Angle orientation) {
        if (orientation == null)
            throw new IllegalArgumentException(NULL_ERROR);
        this.orientation = orientation;
    }

    /**
     * This method returns a hash code for this Geodetic2DEllipse object. Ellipses that
     * have equal parameters (or only differ in orientation by 180 degrees) have the
     * same hash code, since they are considered equal objects.
     *
     * @return a hash code value for this object.
     */
    @Override
    public int hashCode() {
        // Normalize orientation in radians to between -PI and 0.0 inclusive
        double normAng = orientation.inRadians;
        if (normAng > 0) normAng = normAng - Math.PI;
        return 31 * center.hashCode() +
                Double.valueOf(semiMajorAxis).hashCode() ^
                Double.valueOf(semiMinorAxis).hashCode() ^
                Double.valueOf(normAng).hashCode();
    }

    /**
     * This method is used to compare this ellipse to another Geodetic2DEllipse object.
     * Ellipses that have equal parameters (or only differ in orientation by 180 degrees)
     * are considered equal objects.
     *
     * @param ellipse Geodetic2DEllipse to compare to this one
     * @return boolean indicating whether this ellipse is equal to the specified ellipse
     */
    public boolean equals(Geodetic2DEllipse ellipse) {
        if (this == ellipse) {
            return true;
        }
        boolean result = false;
        if (ellipse != null && this.getCenter().equals(ellipse.getCenter())) {
            if ((this.semiMajorAxis == ellipse.semiMajorAxis) &&
                    (this.semiMinorAxis == ellipse.semiMinorAxis)) {
                if (this.orientation.equals(ellipse.orientation))
                    result = true;
                else {
                    Angle temp = new Angle(Math.PI).add(this.orientation);
                    if (temp.equals(ellipse.orientation)) result = true;
                }
            }
        }
        return result;
    }

    /**
     * This method is used to compare this ellipse to another Geodetic2DEllipse object.
     * Ellipses that have equal parameters (or only differ in orientation by 180 degrees)
     * are considered equal objects.
     *
     * @param that Geodetic2DEllipse to compare to this one
     * @return boolean indicating whether this ellipse is equal to the specified ellipse
     */
    @Override
    public boolean equals(Object that) {
        return that instanceof Geodetic2DEllipse && equals((Geodetic2DEllipse) that);
    }

    /**
     * The toString method formats the ellipse for printing.
     *
     * @return String representation of center point, axes, and orientation
     */
    @Override
    public String toString() {
        return "(" + center + "," + semiMajorAxis + "m by " +
                semiMinorAxis + "m at " + orientation + ")";
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/SafeDateFormat.java<|end_filename|>
/****************************************************************************************
 *  SafeDateFormat.java
 *
 *  Created: Apr 10, 2009
 *
 *  @author DRAND
 *
 *  (C) Copyright MITRE Corporation 2009
 *
 *  The program is provided "as is" without any warranty express or implied, including
 *  the warranty of non-infringement and the implied warranties of merchantibility and
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any
 *  damages suffered by you as a result of using the Program.  In no event will the
 *  Copyright owner be liable for any special, indirect or consequential damages or
 *  lost profits even if the Copyright owner has been advised of the possibility of
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.text.ParseException;
import java.text.SimpleDateFormat;
import java.util.Date;
import java.util.TimeZone;

/**
 * {@link SimpleDateFormat} is not a thread safe class and therefore needs to be
 * instantiated on a per-thread basis. This class is a simplistic wrapper that
 * wraps a per-thread formatter and creates them on the fly, exposing only 
 * needed functionality.
 * 
 * @author DRAND
 */
public class SafeDateFormat {

	private final ThreadLocal<SimpleDateFormat> ms_dateFormatter =
		new ThreadLocal<SimpleDateFormat>();
	
	@NonNull private String pattern;
	@NonNull private TimeZone timeZone;

	/**
	 * Constructs a <code>SafeDateFormat</code> using the given pattern
	 * and UTC (aka GMT) as the default time zone.
	 *  
	 * @param pattern the pattern describing the date and time format
	 * @throws IllegalArgumentException if pattern is empty string or null
	 */
	public SafeDateFormat(String pattern) {
		this(pattern, null); // default time zone to UTC
	}

	/**
	 * Constructs a <code>SafeDateFormat</code> using the given pattern and TimeZone.
	 *  
	 * @param pattern the pattern describing the date and time format
	 * @param tz the given new time zone, if null UTC is used by default.
	 * @throws IllegalArgumentException if pattern is empty string, blank or null
	 */
	public SafeDateFormat(String pattern, TimeZone tz) {
		if (pattern == null || pattern.trim().length() == 0) {
			throw new IllegalArgumentException(
					"pattern should never be null or empty");
		}
		this.pattern = pattern;
		this.timeZone = (tz == null) ? TimeZone.getTimeZone("UTC") : tz;
	}
	
	private SimpleDateFormat getInstance() {
		SimpleDateFormat df = ms_dateFormatter.get();
		if (df == null) {
			df = new SimpleDateFormat(pattern);
			df.setTimeZone(timeZone);
			ms_dateFormatter.set(df);
		}
		return df;
	}
	
	/**
	 * Format the value
	 * @param value the value, never <code>null</code>
	 * @return the formatted value, never <code>null</code> or empty
	 */
	@NonNull
	public String format(Date value) {
		return getInstance().format(value);
	}

	/**
	 * Format the value
	 * @param value    The value to format
	 * @return the formatted value, never <code>null</code> or empty
	 * @exception IllegalArgumentException if the Format cannot format the given
	 *            object
	 */
	@NonNull
	public String format(Number value) {
		return getInstance().format(value);
	}
	
	/**
	 * Parse the value
	 * @param value the value, never <code>null</code>
	 * @return the parsed value, never <code>null</code>
	 * @exception ParseException if the beginning of the specified string
	 *            cannot be parsed.
	 */
	@NonNull
	public Date parse(String value) throws ParseException {
		return getInstance().parse(value);
	}

	/**
	 * Sets the time zone for the calendar of this SafeDateFormat object.
	 * 
     * @param timeZone the given new time zone.
	 */
	public void setTimeZone(TimeZone timeZone) {
		if (timeZone == null) timeZone = TimeZone.getTimeZone("UTC");
		this.timeZone = timeZone;
		getInstance().setTimeZone(timeZone);
	}

	/**
     * Gets the time zone.
     * @return the time zone associated with the calendar of this SafeDateFormat.
     */
	@NonNull
	public TimeZone getTimeZone() {
		return timeZone;
	}

	public String toPattern() {
		return pattern;
	}
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Angle.java<|end_filename|>
/****************************************************************************************
 *  Angle.java
 *
 *  (C) Copyright MITRE Corporation 2007-2008
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.io.Serializable;
import java.text.DecimalFormat;
import java.util.StringTokenizer;

/**
 * The Angle class is used to store and manipulate angular quantities. By having
 * conventions for the position of zero on the circle and the direction of increasing
 * angular value, angles can be used to represent either positions on the circle, or
 * angular distances between positions on the circle such that those distances are
 * less than one full revolution. As such, this class maintains the value of an angle
 * internally in normalized radians between -PI inclusive and +PI exclusive, but has
 * constructors for decimal radians or degrees, numeric components of degrees-minutes-
 * seconds, and a String parser for degrees, minutes, and seconds. There are accessor
 * methods for getting the value of this Angle instance in decimal degrees or in radians.
 * Also provided is a method to add another Angle to this Angle instance, a method to
 * find the smallest angular difference and direction from this Angle instance. Longitude
 * and Latitude are objects that extend Angle with additional methods and constraints.
 * The toString methods will pretty print the angle in a degrees, minutes, and seconds
 * String format which the appropriate constructor can parse if needed. An Angle can
 * be constructed from a String, and the degree symbol '&deg;' can be encoded either
 * using the ISO-8859-1 (Latin-1 Western) character encoding or in Unicode. The toString
 * method currently writes the degree symbol using the Latin-1 Western encoding. This
 * usually gets escaped as an HTML entity for proper viewing across platforms.
 *
 * @author <NAME>
 */
public class Angle implements Serializable, Comparable<Angle> {
    private static final long serialVersionUID = 1L;

    // Instance variable holds angular value in normalized radians (-PI .. +PI)
    double inRadians;

    // Constant values
    static final double TWO_PI = 2.0 * Math.PI;  // 360 degrees

    /**
     * Unicode degree symbol
     * <br>
     * This is used by {@link #toString()}.
     */
    public static final String DEGSYM = Character.toString('\u00B0'); //

    /**
     * unitType parameter for constructor specifying angle is in decimal radians
     */
    public static final int RADIANS = 0;
    /**
     * unitType parameter for constructor specifying angle is in decimal degrees
     */
    public static final int DEGREES = 1;

    /*
     * Reduce angle in radians to normalized range form [-PI .. +PI)
     *
     * @param rad value in radians
     * @return normalized value in radians (-PI .. PI)
     * @throws IllegalArgumentException - error if angle is too large
     */
    static double normalize(double rad) {
        // This is somewhat arbitrary to limit angles to 4 revolutions
        if (Math.abs(rad) > (4.0 * TWO_PI))
            throw new IllegalArgumentException("Angle " + rad + " radians is too big");
        while (rad >= Math.PI) rad -= TWO_PI;
        while (rad < -Math.PI) rad += TWO_PI;
        return rad;
    }

    /*
     * Base initializer takes angle value and units designator
     */
    protected void init(double value, int unitType) {
        if (unitType == Angle.RADIANS)
            this.inRadians = normalize(value);
        else if (unitType == Angle.DEGREES)
            this.inRadians = normalize(Math.toRadians(value));
        else throw new IllegalArgumentException("Invalid Angle unit type");
    }

    /*
     * This initializer takes integer degrees and minutes with decimal seconds
     * Degrees may be any signed integer value, others are range checked.
     */
    protected void init(double deg, double min, double sec) {
        if ((min < 0.0) || (60.0 <= min))
            throw new IllegalArgumentException("Arc minutes value '" + min +
                    "' is out of legal range");
        else if ((sec < 0.0) || (60.0 <= sec))
            throw new IllegalArgumentException("Arc seconds value '" + sec +
                    "' is out of legal range");
        double ddeg = Math.abs(deg) + (min / 60.0) + (sec / 3600.0);
        double sign = (deg < 0.0) ? -1.0 : +1.0;
        this.init(sign * ddeg, Angle.DEGREES);
    }

    /**
     * This constructor takes angle value and unitType designator.
     *
     * @param value    angular value in specified units (Angle.DEGREES or Angle.RADIANS)
     * @param unitType one of Angle.DEGREES or Angle.RADIANS
     * @throws IllegalArgumentException error if units or value are out of range
     */
    public Angle(double value, int unitType) {
        this.init(value, unitType);
    }

    /**
     * This constructor assumes units are radians.
     *
     * @param valueInRadians angular value in decimal radians
     */
    public Angle(double valueInRadians) {
        this.init(valueInRadians, Angle.RADIANS);
    }

    /**
     * This constructor takes numeric degree, minutes, seconds components.
     *
     * @param degrees signed, integer degree component of angular value
     * @param minutes integer minutes component of angular value
     * @param seconds double seconds component of angular value
     * @throws IllegalArgumentException error if a component value is out of range
     */
    public Angle(int degrees, int minutes, double seconds) {
        this.init(degrees, minutes, seconds);
    }

    /**
     * This constructor takes a string and parses it as an angle.  It breaks the string into
     * one, two, or three numeric components with optional units designation.  Units are assumed
     * to be in degrees unless marked as radians, and components after the first are assumed to
     * be arc minutes and arc seconds.
     *
     * @param valStr the String to be parsed as an instance of Angle
     * @throws IllegalArgumentException error if string can not be successfully parsed.
     */
    public Angle(String valStr) {
        String token;
        String[] components = new String[]{"", "", ""};
        double ang, min, sec;

        // Normalize to upper case chars, allow several possible delimiters
        valStr = valStr.trim().toUpperCase();
        StringTokenizer st = new StringTokenizer(valStr, " \t,;:");

        // Combine unit tokens with their preceding number tokens using optimistic assumptions
        // This allows optional whitespace between numeric values and unit strings to be accepted
        int i = 0;
        while (st.hasMoreTokens() && (i < 3)) {
            token = st.nextToken();
            boolean containsDigits = false;
            for (int j = 0; j < token.length(); j++) {
                if (Character.isDigit(token.charAt(j))) {
                    containsDigits = true;
                    components[i++] = token;
                    break;
                }
            }
            if (!containsDigits && (i > 0)) components[i - 1] += token;
        }
        if (st.hasMoreTokens())
            throw new IllegalArgumentException("Unexpected token '" + st.nextToken() +
                    "' in Angle string \"" + valStr + "\"");

        // Note which components we actually have
        boolean hasAng = (components[0].length() > 0);
        boolean hasMin = (components[1].length() > 0);
        boolean hasSec = (components[2].length() > 0);
        if (hasMin && !hasSec && components[1].endsWith("\"")) {
            // special case: minutes and seconds without whitespace in-between
            // e.g., 37 25'19.07"N
            int ind = components[1].indexOf('\'');
            if (ind > 0) {
                components[2] = components[1].substring(ind + 1);
                components[1] = components[1].substring(0, ind + 1);
                hasSec = true;
            }
        }

        if (!hasAng) throw new IllegalArgumentException("No valid tokens found in Angle string \"" +
                valStr + "\"");

        // No units suffix assumes degrees
        int units = Angle.DEGREES;
        // Look for angle in degrees or radians
        token = components[0];
        // See if there is a units designator string suffix and handle accordingly
        int suffixLen = 0;
        if (token.endsWith(DEGSYM)) suffixLen = 1;
        else if (token.endsWith("D")) suffixLen = 1;
        else if (token.endsWith("DEG")) suffixLen = 3;
        else if (token.endsWith("DEGS")) suffixLen = 4;
        else if (token.endsWith("DEGREES")) suffixLen = 7;
        else if (token.endsWith("R")) {
            units = Angle.RADIANS;
            suffixLen = 1;
        } else if (token.endsWith("RAD")) {
            units = Angle.RADIANS;
            suffixLen = 3;
        } else if (token.endsWith("RADS")) {
            units = Angle.RADIANS;
            suffixLen = 4;
        } else if (token.endsWith("RADIANS")) {
            units = Angle.RADIANS;
            suffixLen = 7;
        }
        // Remove the units suffix string (if any) and parse the numeric part for angle
        token = token.substring(0, token.length() - suffixLen);
        boolean negativeZeroAngle = false;
        if ("-0".equals(token)) {
            //System.out.println("XXX: negative zero");
            // don't want to deal with negative zero in double calculations: flip sign later
            negativeZeroAngle = true;
            ang = 0;
        } else
            ang = Double.parseDouble(token);

        if (hasMin) {
            // Do some additional error checking if there is more than one component
            if (units == Angle.RADIANS)
                throw new IllegalArgumentException("Unexpected tokens following angle in radians, in \"" +
                        valStr + "\"");
            else if (components[0].contains("."))
                throw new IllegalArgumentException("Only least-significant component of angle " +
                        "can be fractional, in \"" + valStr + "\"");

            // Since we now know our units are degrees, look for possible minutes or minutes and seconds
            token = components[1];
            // See if there is a units designator string suffix and handle accordingly
            suffixLen = 0;
            if (token.endsWith("'")) suffixLen = 1;
            else if (token.endsWith("M")) suffixLen = 1;
            else if (token.endsWith("MIN")) suffixLen = 3;
            else if (token.endsWith("MINS")) suffixLen = 4;
            else if (token.endsWith("MINUTES")) suffixLen = 7;
            // Remove the units suffix string (if any) and parse the numeric part for minutes
            token = token.substring(0, token.length() - suffixLen);
            min = Double.parseDouble(token);
            if (min < 0.0)
                throw new IllegalArgumentException("Illegal negative value for arc minutes of angle in \"" +
                        valStr + "\"");

            // Now look for possible additional value of arc seconds
            if (hasSec) {
                if (components[1].contains("."))
                    throw new IllegalArgumentException("Only least-significant component of angle " +
                            "can be fractional, in \"" + valStr + "\"");

                token = components[2];
                // See if there is a units designator string suffix and handle accordingly
                suffixLen = 0;
                if (token.endsWith("\"")) suffixLen = 1;
                else if (token.endsWith("S")) suffixLen = 1;
                else if (token.endsWith("SEC")) suffixLen = 3;
                else if (token.endsWith("SECS")) suffixLen = 4;
                else if (token.endsWith("SECONDS")) suffixLen = 7;
                // Remove the units suffix string (if any) and parse the numeric part for minutes
                token = token.substring(0, token.length() - suffixLen);
                sec = Double.parseDouble(token);
                if (sec < 0.0)
                    throw new IllegalArgumentException("Illegal negative value for arc seconds of angle in \"" +
                            valStr + "\"");
            } else
                sec = 0.0;
            // Initialize Angle instance using degree, minute, second form
            this.init(ang, min, sec);
        } else
            // Initialize Angle instance using decimal degrees or radians
            this.init(ang, units);

        if (negativeZeroAngle) {
            // swap the sign
            inRadians = -inRadians;
        }
    }

    /**
     * Default constructor makes an angle whose value is zero radians.
     */
    public Angle() {
        this.inRadians = 0.0;
    }

    /**
     * This accessor method returns the value of the Angle in decimal radians.
     *
     * @return value of this angle in radians
     */
    public double inRadians() {
        return this.inRadians;
    }

    /**
     * This accessor method returns the value of the Angle in decimal radians. It is
     * equalivent to the inRadians method, but it follows the standard JavaBean getter
     * signature, used by Frameworks like Hibernate and Spring.
     *
     * @return value of this angle in radians
     */
    double getRadians() {
        return this.inRadians;
    }

    /**
     * This settor method is used to update the value of this Angle by directly setting
     * its value in decimal radians. It is used by Frameworks like Hibernate and Spring
     *
     * @param radians new value in decimal radians for this Angle
     */
    void setRadians(double radians) {
        this.inRadians = radians;
    }

    /**
     * This accessor method returns the value of the Angle in decimal degrees.
     *
     * @return value of this angle in degrees
     */
    public double inDegrees() {
        return Math.toDegrees(this.inRadians);
    }

    /**
     * This accessor method returns the value of the Angle in the specified unitType.
     *
     * @param unitType type of units desired (one of Angle.RADIANS or Angle.DEGREES)
     * @return value of this angle in the decimal radians or degrees, depending on unitType specified
     * @throws IllegalArgumentException error if unitType is invalid
     */
    public double getValue(int unitType) {
        if (unitType == Angle.RADIANS) return this.inRadians();
        else if (unitType == Angle.DEGREES) return this.inDegrees();
        else throw new IllegalArgumentException("Invalid Angle unit type");
    }

    /**
     * This method computes and returns a new angle which is the result of adding the
     * specified angle to this one.
     *
     * @param that Angle to add to this one
     * @return a new Angle that is the result of adding input angle to this one
     */
    @NonNull
    public Angle add(Angle that) {
        return new Angle(this.inRadians + that.inRadians);
    }

    /**
     * This method computes and returns a new angle which is the the shortest angular
     * difference between the specified angle and this one.  The sign indicates which
     * direction (+ is counter-clockwise and - is clockwise, by convention) to go to
     * get from this Angle to that Angle.
     *
     * @param that - Angle to compare to this one
     * @return a new Angle that is the angular difference between input angle and this one
     */
    @NonNull
    public Angle difference(Angle that) {
        return new Angle(that.inRadians - this.inRadians);
    }

    // use 1e-8 for ~ 0.1 meter resolution distinction in radians (on earth's surface)
    private static final double EPSILON = 1e-8;

    /**
     * This method returns a hash code for this Angle object. The result depends only on
     * the angluar value in radians, so that equal Angles have equal hash codes.
     *
     * @return a hash code value for this object.
     */
    @Override
    public int hashCode() {
        // Note we're using approximate equals vs absolute equals on floating point number
        // so must ignore beyond ~7 decimal places in computing the hashCode, otherwise
        // we break the equals-hashCode contract. Changing EPSILON or equals(Angle) may
        // require changing the logic used here also.
        return (int) (inRadians * 1e+7);
    }

    /**
     * equals() method that determines equality of two double precision values that
     * works over a wide range of magnitudes, infinities, NaNs, and positive and
     * negative zero.
     * <p/>
     * Because floating point calculations may involve rounding, calculated float and double
     * values may not be accurate. For values that must be precise, such as monetary values,
     * consider using a fixed-precision type such as BigDecimal. For values that need not be
     * precise, consider comparing for equality within some range.
     * <p/>
     * See the Java Language Specification, section 4.2.4.
     * <p/>
     * See http://www.oreilly.com/catalog/javacook/errata/javacook.unconfirmed#Example_5-1
     * http://findbugs.sourceforge.net/bugDescriptions.html#FE_FLOATING_POINT_EQUALITY
     * http://bobcat.webappcabaret.net/javachina/faq/08.htm#math_Q19
     *
     * @param a the first double to compare with.
     * @param b the second double to compare with.
     * @return <code>true</code> if the objects are the same;
     *         <code>false</code> otherwise.
     */
    static boolean equals(double a, double b) {
        /*
        Comparing the difference of two numbers with epsilon requires that the
        magnitude of epsilon change with the magnitude of the numbers being
        compared. But the whole purpose of floating point is that the point floats,
        that is, the programmer need not keep track of the magnitude of the numbers
         */
        return a == b || Math.abs(a - b) < EPSILON * Math.max(Math.abs(a), Math.abs(b));
    }

    /**
     * equals() method that determines equality of two double precision values that
     * works over a wide range of magnitudes, infinities, NaNs, and positive and
     * negative zero.
     * <p/>
     * Because floating point calculations may involve rounding, calculated float and double
     * values may not be accurate. For values that must be precise, such as monetary values,
     * consider using a fixed-precision type such as BigDecimal. For values that need not be
     * precise, consider comparing for equality within some range.
     * <p/>
     * See the Java Language Specification, section 4.2.4.
     * <p/>
     * See http://www.oreilly.com/catalog/javacook/errata/javacook.unconfirmed#Example_5-1
     * http://findbugs.sourceforge.net/bugDescriptions.html#FE_FLOATING_POINT_EQUALITY
     * http://bobcat.webappcabaret.net/javachina/faq/08.htm#math_Q19
     *
     * @param a       the first double to compare with.
     * @param b       the second double to compare with.
     * @param epsilon the epsilon value for error threshold
     * @return <code>true</code> if the objects are the same;
     *         <code>false</code> otherwise.
     */
    static boolean equals(double a, double b, double epsilon) {
        return a == b || Math.abs(a - b) < epsilon * Math.max(Math.abs(a), Math.abs(b));
    }

    /**
     * This method is used to test whether another Angle is equal to this Angle by testing
     * that their values in radians differ in value by no more than epsilon (1e-7).
     *
     * @param that Angle to compare against this one.
     * @return true if the specified Angle is equal to this Angle.
     */
    public boolean equals(Angle that) {
        return that != null && equals(this.inRadians, that.inRadians);
    }

    // Inherited Javadoc
    @Override
    public boolean equals(Object that) {
        return that instanceof Angle && equals((Angle) that);
    }

    /**
     * This method returns a formatted String representation of the angle value in degree,
     * minute, seconds format, with the specified number of fractional digits of arc seconds.
     *
     * @param fractDigOfSec number of fraction digits of seconds to format
     * @return String formatted as integer degrees, minutes, and seconds with optional fractional part
     */
    @NonNull
    public String toString(int fractDigOfSec) {
        int sign = (this.inRadians < 0) ? -1 : +1;
        double rem = Math.abs(Math.toDegrees(this.inRadians));
        int deg = (int) rem;
        rem = (rem - (double) deg) * 60.0;
        int min = (int) rem;
        double sec = (rem - (double) min) * 60.0;

        String fract = (fractDigOfSec > 0) ? "." : "";
        while (fractDigOfSec-- > 0) fract += "0";
        DecimalFormat secFmt = new DecimalFormat("0" + fract);
        if (secFmt.format(sec).startsWith("60")) {
            min += 1;
            sec = 0.0;
        }
        if (min == 60) {
            deg += 1;
            min = 0;
        }
        return Integer.toString(sign * deg) + DEGSYM + " " + Integer.toString(min) +
                "' " + secFmt.format(sec) + "\"";
    }

    /**
     * This method returns a formatted String representation of the angle value in degree,
     * minute, seconds format, with no fractional digits of arc seconds.
     *
     * @return String formatted as integer degrees, minutes, and seconds
     */
    @Override
    public String toString() {
        return this.toString(0);
    }

    /**
     * @param that Angle to compare lengths with this one
     * @return a negative integer, zero, or a positive integer as this object is less than,
     *         equal to, or greater than the specified object.
     * @throws NullPointerException if that is null
     */
    public int compareTo(Angle that) {
        return (this.inRadians < that.inRadians) ? -1 :
                (this.inRadians > that.inRadians) ? +1 : 0;
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Longitude.java<|end_filename|>
/****************************************************************************************
 *  Longitude.java
 *
 *  Created: Dec 28, 2006
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

/**
 * The Longitude class extends Angle to provide parsing and formatting with direction indicator
 * suffix ('E' or 'W') from the Prime Meridian instead of numeric sign. Note that the normalized
 * form for Angles maps +180 degrees to its equivalent -180 degree form at the International
 * Date Line.
 */
public class Longitude extends Angle {

	private static final long serialVersionUID = 1L;

    /**
     * This constructor takes longitude angle value and units designator.
     *
     * @param value    Longitude angular value in specified units (Angle.DEGREES or Angle.RADIANS)
     * @param unitType Angle.DEGREES or Angle.RADIANS
     * @throws IllegalArgumentException error if units or value are out of range
     */
    public Longitude(double value, int unitType) {
        this.init(value, unitType);
    }

    /**
     * This constructor takes a longitude formatted as a string.  It looks for an 'E' or 'W'
     * direction indicator suffix character and sets the numeric sign accordingly (E is +, W is -).
     * Units are assumed to be in degrees unless marked as radians.
     *
     * @param lonStr longitude value string
     * @throws IllegalArgumentException error if syntax is incorrect or value is out of range
     */
    public Longitude(String lonStr) {
		if (lonStr == null) throw new IllegalArgumentException("lonStr cannot be null");
        lonStr = lonStr.trim().toUpperCase();
        int n = lonStr.length() - 1;
		if (n < 0) throw new IllegalArgumentException("Longitude must be non-empty value");
        char suffix = lonStr.charAt(n);
        char sign = ' ';
        if (suffix == 'E') sign = '+';
        else if (suffix == 'W') sign = '-';
        else if ((suffix == 'N') || (suffix == 'S'))
            throw new IllegalArgumentException("Longitude with N or S direction indicator");
        if (sign != ' ') {
            // Make sure there is no numeric sign prefix when direction char suffix is present
            char prefix = lonStr.charAt(0);
            if ((prefix == '+') || (prefix == '-'))
                throw new IllegalArgumentException("Longitude with E or W direction indicator " +
                        "should not have numeric sign prefix.");
            // Strip off the suffix and add the sign based on direction indicator
            lonStr = sign + lonStr.substring(0, n);
        }
        this.init(new Angle(lonStr).inRadians(), Angle.RADIANS);
    }

    /**
     * This constructor assumes Longitude angle units are in decimal radians.
     * Longitude value is normalized to the (-PI .. PI) range and by convention
     * +180 becomes -180.
     *
     * @param valueInRadians Longitude angular value in radians
     */
    public Longitude(double valueInRadians) {
        this.init(valueInRadians, Angle.RADIANS);
    }

    /**
     * This constructor takes numeric degree, minutes, seconds components for Longitude.
     *
     * @param degrees signed, integer degree component of Longitude angular value
     * @param minutes integer minutes component of Longitude angular value
     * @param seconds double seconds component of Longitude angular value
     * @throws IllegalArgumentException error if a component value is out of range
     */
    public Longitude(int degrees, int minutes, double seconds) {
        this.init(degrees, minutes, seconds);
    }

    /**
     * This constructor takes an Angle object for Longitude.
     *
     * @param angle Angle to use as basis for this Longitude
     */
    public Longitude(Angle angle) {
        this.init(angle.inRadians(), Angle.RADIANS);
    }

    /**
     * The default constructor makes a Longitude object at the Prime Meridian (value of zero radians).
     */
    public Longitude() {
        this.inRadians = 0.0;
    }

    /**
     * This accessor method returns the character for the direction from the
     * Prime Meridian of this Longitude ('E' or 'W').
     *
     * @return char indicating westerly ('W') or easterly ('E') direction
     */
    public char getDirectionFromPrimeMeridian() {
        return (this.inRadians < 0) ? 'W' : 'E';
    }

    /**
     * This method computes the distance in radians from this Longitude to the specified
     * eastLon Longitude by travelling eastward (positive increasing angular direction).
     *
     * @param eastLon destination Longitude around globe
     * @return positive radians needed to sweep east to get from this Longitude to eastLon
     */
    public double radiansEast(Longitude eastLon) {
        double diff = this.difference(eastLon).inRadians;
        return (diff < 0.0) ? Angle.TWO_PI + diff : diff;
    }

    /**
     * This predicate method determines whether this Longitude is within the test
     * interval specified from westLon to eastLon.
     *
     * @param westLon western most Longitude of the interval being tested
     * @param eastLon eastern most Longitude of the interval being tested
     * @return boolean indicating whether this Longitude is in the interval
     */
    public boolean inInterval(Longitude westLon, Longitude eastLon) {
        double westRad = westLon.inRadians;
        double eastRad = eastLon.inRadians;
        double testRad = this.inRadians;
        // If interval wraps around, break it into two parts and test each separately
        if (westRad <= eastRad) {
            return ((westRad <= testRad) && (testRad <= eastRad));
        } else {
            return (((westRad <= testRad) && (testRad < Math.PI)) ||
                    ((-Math.PI <= testRad) && (testRad <= eastRad)));
        }
    }

    /**
     * This method returns a string version of this longitude value in degree, minute,
     * seconds format, with direction character suffix.
     */
    public String toString(int fractDigOfSec) {
        String result = super.toString(fractDigOfSec);
        char dirFromMeridian = this.getDirectionFromPrimeMeridian();
        if (result.charAt(0) == '-') result = result.substring(1);
        // Normalize rounded +180 deg values to be -180 deg by display convention
        if ((result.substring(0, 3).equals("180")) && (dirFromMeridian == 'E'))
            dirFromMeridian = 'W';
        return result + " " + dirFromMeridian;
    }
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestGeodetic2DArc.java<|end_filename|>
package org.opensextant.geodesy.test;

import junit.framework.TestCase;
import org.opensextant.geodesy.Angle;
import org.opensextant.geodesy.Geodetic2DArc;
import org.opensextant.geodesy.Geodetic2DPoint;
import org.opensextant.geodesy.Latitude;
import org.opensextant.geodesy.Longitude;

/**
 * Simple tests for Geodetic2DArc
 * @author <NAME>, MITRE Corp.
 */
public class TestGeodetic2DArc extends TestCase {

	public void testCreate() {
		Geodetic2DPoint west = new Geodetic2DPoint(new Longitude(30, Angle.DEGREES),
				new Latitude(30, Angle.DEGREES));
		Geodetic2DPoint east = new Geodetic2DPoint(new Longitude(31, Angle.DEGREES),
				new Latitude(31, Angle.DEGREES));
		Geodetic2DArc arc = new Geodetic2DArc(west, east);
		assertNotNull(arc.getPoint1());
		assertNotNull(arc.getPoint2());
		assertNotNull(arc.toString());
		assertEquals(146647.2, arc.getDistanceInMeters(), 0.001);
		assertEquals(41, Math.round(arc.getForwardAzimuth().inDegrees()));

		Geodetic2DArc arc2 = new Geodetic2DArc(east, east);
		assertEquals(0.0, arc2.getDistanceInMeters(), 1e-6);
	}

	public void testEquals() {
		Geodetic2DPoint west = new Geodetic2DPoint(new Longitude(30, Angle.DEGREES),
				new Latitude(30, Angle.DEGREES));
		Geodetic2DPoint east = new Geodetic2DPoint(new Longitude(32, Angle.DEGREES),
				new Latitude(32, Angle.DEGREES));
		Geodetic2DArc arc = new Geodetic2DArc(west, east);

		Geodetic2DArc arc2 = new Geodetic2DArc(east, west);
		assertEquals(arc.getDistanceInMeters(), arc2.getDistanceInMeters(), 10-6);

		Geodetic2DArc arc3 = new Geodetic2DArc(west, arc.getDistanceInMeters(),
				arc.getForwardAzimuth());
		assertEquals(arc, arc3);
		assertEquals(0, arc.compareTo(arc3));

		// change point to make the object fail equality
		arc3.setPoint2(new Geodetic2DPoint(new Longitude(33, Angle.DEGREES),
				new Latitude(33, Angle.DEGREES)));
		assertFalse(arc.equals(arc3));
	}

	public void testNotEquals() {
		Geodetic2DPoint west = new Geodetic2DPoint(new Longitude(30, Angle.DEGREES),
				new Latitude(30, Angle.DEGREES));
		Geodetic2DPoint east = new Geodetic2DPoint(new Longitude(32, Angle.DEGREES),
				new Latitude(32, Angle.DEGREES));
		Geodetic2DArc arc = new Geodetic2DArc(west, east);
		assertFalse(arc.equals(null));
		Object other = this;
		assertFalse(arc.equals(other));
	}

}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Geodetic2DArc.java<|end_filename|>
/****************************************************************************************
 *  Geodetic2DArc.java
 *
 *  Created: Oct 31, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.io.Serializable;

/**
 * The Geodetic2DArc class represents two geodetic points on the surface of an
 * Ellipsoid, the distance in meters on the surface between them, and the forward
 * azimuth (Angle from North to use as heading from point 1 to point 2). It includes
 * a constructor that takes the two points and then calculates the distance and
 * azimuth, as well as a constructor that takes the first point, distance, and
 * azimuth and computes the second point. There are variations that assume the
 * WGS-84 Ellipsoid and ones that allow the user to specify which Ellipsoid to use.
 * This code was adapted from the open source GeoTools project, having been
 * translated to Java from original Fortran published by NOAA.
 */
public class Geodetic2DArc implements Serializable, Comparable<Geodetic2DArc> {
	
    private static final long serialVersionUID = 1L;

    private static final Ellipsoid WGS84 = Ellipsoid.getInstance("WGS 84");
    private static final String NOCONVERGE = "Failed to converge for arc from ";

    // Tolerance factors from the strictest to the most relaxed
    private static final double TOLERANCE_0 = 5.0e-15,
            TOLERANCE_1 = 5.0e-14,
            TOLERANCE_2 = 5.0e-13,
            TOLERANCE_3 = 7.0e-3;

    @NonNull private Ellipsoid ellip;
    @NonNull private Geodetic2DPoint point1;
    @NonNull private Geodetic2DPoint point2;
    @NonNull private Angle forwardAzimuth;
    private double distanceInMeters;

    private double semiMajorAxis, semiMinorAxis, eccentricitySquared;
    private double f, fo;
    private double A, B, C, D, E, F;
    private double maxOrthodromicDistance;
    private double T1, T2, T4, T6;
    private double a01, a02, a03, a21, a22, a23, a42, a43, a63;

    /**
     * The default constructor creates a degenerate Geodetic2DArc with both
     * points at (0,0).
     */
    public Geodetic2DArc() {
        Geodetic2DPoint pt = new Geodetic2DPoint();
        this.ellip = WGS84;
        init();

        this.point1 = pt;
        setPoint2(pt); // updates state (distance and azimuth)
    }

    /**
     * This constructor assumes the WGS-84 Ellipsoid model of the earth, takes two geodetic
     * points, and computes the forward azimuth and distance in meters between them.
     *
     * @param point1 first Geodetic2DPoint (fixed)
     * @param point2 second Geodetic2DPoint (re-settable)
	 * @throws NullPointerException if point1 or point2 are null
     */
    public Geodetic2DArc(Geodetic2DPoint point1, Geodetic2DPoint point2) {
        this(WGS84, point1, point2);
    }

    /**
     * This constructor takes the Ellipsoid and two geodetic points, and computes the
     * forward azimuth and distance in meters between them.
     *
     * @param ellip  Ellipsoid model of the earth to use
     * @param point1 first Geodetic2DPoint (fixed)
     * @param point2 second Geodetic2DPoint (re-settable)
	 * @throws NullPointerException if ellip, point1 or point2 are null
     */
    public Geodetic2DArc(Ellipsoid ellip, Geodetic2DPoint point1, Geodetic2DPoint point2) {
        this.ellip = ellip;
        init();

        this.point1 = point1;
        setPoint2(point2); // updates state (distance and azimuth)
    }

    /**
     * This constructor assumes the WGS-84 Ellipsoid model of the earth, takes the first
     * geodetic point, an Angle of forward azimuth, and a distance in meters, and then
     * computes the second geodetic point that lies at that distance along that azimuth.
     *
     * @param point1           first Geodetic2DPoint (fixed)
     * @param forwardAzimuth   Angle from North of azimuth from point 1 to point 2
     * @param distanceInMeters double distance in meters on surface of Ellipsoid
     * @throws IllegalArgumentException if error in argument values is detected
     */
    public Geodetic2DArc(Geodetic2DPoint point1, double distanceInMeters,
                         Angle forwardAzimuth) {
        this(WGS84, point1, distanceInMeters, forwardAzimuth);
    }

    /**
     * This constructor takes the Ellipsoid model of the earth, the first geodetic point,
     * an Angle of forward azimuth, and a distance in meters, and then computes the second
     * geodetic point that lies at that distance along that azimuth.
     *
     * @param ellip            Ellipsoid model of the earth to use
     * @param point1           first Geodetic2DPoint (fixed)
     * @param distanceInMeters double distance in meters on surface of Ellipsoid
     * @param forwardAzimuth   Angle from North of azimuth from point 1 to point 2
     * @throws IllegalArgumentException if error in argument values is detected
	 * @throws NullPointerException if ellip, point1 or forwardAzimuth are null
     */
    public Geodetic2DArc(Ellipsoid ellip, Geodetic2DPoint point1, double distanceInMeters,
                         Angle forwardAzimuth) {
        this.ellip = ellip;
        init();

        this.point1 = point1;
        setDistanceAndAzimuth(distanceInMeters, forwardAzimuth); // updates state (point2)
    }

    // Initialization of parameters based on the Ellipsoid
    private void init() {
        semiMajorAxis = ellip.getEquatorialRadius();
        semiMinorAxis = ellip.getPolarRadius();
        eccentricitySquared = ellip.getEccentricitySquared();

        // Calculate needed parameters based on the ellipsoid
        f = ellip.getFlattening();
        fo = 1.0 - f;
        double f2 = f * f;
        double f3 = f * f2;
        double f4 = f * f3;

        final double E2 = eccentricitySquared;
        final double E4 = E2 * E2;
        final double E6 = E4 * E2;
        final double E8 = E6 * E2;
        final double EX = E8 * E2;

        A = 1.0 + 0.75 * E2 + 0.703125 * E4 + 0.68359375 * E6 + 0.67291259765625 * E8 + 0.6661834716796875 * EX;
        B = 0.75 * E2 + 0.9375 * E4 + 1.025390625 * E6 + 1.07666015625 * E8 + 1.1103057861328125 * EX;
        C = 0.234375 * E4 + 0.41015625 * E6 + 0.538330078125 * E8 + 0.63446044921875 * EX;
        D = 0.068359375 * E6 + 0.15380859375 * E8 + 0.23792266845703125 * EX;
        E = 0.01922607421875 * E8 + 0.0528717041015625 * EX;
        F = 0.00528717041015625 * EX;

        maxOrthodromicDistance = semiMajorAxis * (1.0 - E2) * Math.PI * A - 1.0;

        T1 = 1.0;
        T2 = -0.25 * f * (1.0 + f + f2);
        T4 = 0.1875 * f2 * (1.0 + 2.25 * f);
        T6 = 0.1953125 * f3;

        final double a = f3 * (1.0 + 2.25 * f);
        a01 = -f2 * (1.0 + f + f2) / 4.0;
        a02 = 0.1875 * a;
        a03 = -0.1953125 * f4;
        a21 = -a01;
        a22 = -0.25 * a;
        a23 = 0.29296875 * f4;
        a42 = 0.03125 * a;
        a43 = 0.05859375 * f4;
        a63 = 5.0 * f4 / 768.0;
    }

    /**
     * Calculate the meridian arc length between two points in the same meridian
     * in the referenced ellipsoid.
     *
     * @param P1 The latitude of the first  point (in radians).
     * @param P2 The latitude of the second point (in radians).
     * @return Returned the meridian arc length between P1 and P2
     */
    private double getMeridianArcLengthRadians(final double P1, final double P2) {
        /*
         * Latitudes P1 and P2 in radians positive North and East.
         * Forward azimuths at both points returned in radians from North.
         *
         * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/inverse.for
         *         subroutine GPNARC
         *         version    200005.26
         *         written by Robert (<NAME>
         *
         * Ported from Fortran to Java by <NAME>.
         */
        double S1 = Math.abs(P1);
        double S2 = Math.abs(P2);
        double DA = (P2 - P1);
        // Check for a 90 degree lookup
        if (S1 > TOLERANCE_0 || S2 <= (Math.PI / 2 - TOLERANCE_0) || S2 >= (Math.PI / 2 + TOLERANCE_0)) {
            final double DB = Math.sin(P2 * 2.0) - Math.sin(P1 * 2.0);
            final double DC = Math.sin(P2 * 4.0) - Math.sin(P1 * 4.0);
            final double DD = Math.sin(P2 * 6.0) - Math.sin(P1 * 6.0);
            final double DE = Math.sin(P2 * 8.0) - Math.sin(P1 * 8.0);
            final double DF = Math.sin(P2 * 10.0) - Math.sin(P1 * 10.0);
            // Compute the S2 part of the series expansion
            S2 = -DB * B / 2.0 + DC * C / 4.0 - DD * D / 6.0 + DE * E / 8.0 - DF * F / 10.0;
        }
        // Compute the S1 part of the series expansion
        S1 = DA * A;
        // Compute the arc length
        return Math.abs(semiMajorAxis * (1.0 - eccentricitySquared) * (S1 + S2));
    }

    private void calcDistanceAndAzimuth() {
        final double lat1 = point1.getLatitude().inRadians();
        final double lat2 = point2.getLatitude().inRadians();

        Angle lon1 = point1.getLongitude();
        Angle lon2 = point2.getLongitude();
        final double dlon = lon1.difference(lon2).inRadians();

        /*
         * Solution of the geodetic inverse problem after T.Vincenty.
         * Modified Rainsford's method with Helmert's elliptical terms.
         * Effective in any azimuth and at any distance short of antipodal.
         *
         * Latitudes and longitudes in radians positive North and East.
         * Forward azimuths at both points returned in radians from North.
         *
         * Programmed for CDC-6600 by <NAME> NGS ROCKVILLE MD 18FEB75
         * Modified for IBM SYSTEM 360 by <NAME> NGS ROCKVILLE MD 7507
         * Ported from Fortran to Java by <NAME>.
         *
         * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/inverse.for
         *         subroutine GPNHRI
         *         version    200208.09
         *         written by robert (sid) safford
         */
        final double ss = Math.abs(dlon);
        if (ss < TOLERANCE_1) {
            distanceInMeters = getMeridianArcLengthRadians(lat1, lat2);
            forwardAzimuth = new Angle((lat2 > lat1) ? 0.0 : Math.PI);
            return;
        }
        /*
         * Compute the limit in longitude (alimit), it is equal
         * to twice  the distance from the equator to the pole,
         * as measured along the equator
         */
        //test for antinodal difference
        final double ESQP = eccentricitySquared / (1.0 - eccentricitySquared);
        final double alimit = Math.PI * fo;
        if (ss >= alimit &&
                lat1 < TOLERANCE_3 && lat1 > -TOLERANCE_3 &&
                lat2 < TOLERANCE_3 && lat2 > -TOLERANCE_3) {
            // Compute an approximate AZ
            final double CONS = (Math.PI - ss) / (Math.PI * f);
            double AZ = Math.asin(CONS);
            double AZ_TEMP, S, AO;
            int iter = 0;
            do {
                if (++iter > 8) {
                    String pts = point1.toString() + " to " + point2.toString();
                    throw new ArithmeticException(NOCONVERGE + pts);
                }
                S = Math.cos(AZ);
                final double C2 = S * S;
                // Compute new AO
                AO = T1 + T2 * C2 + T4 * C2 * C2 + T6 * C2 * C2 * C2;
                final double CS = CONS / AO;
                S = Math.asin(CS);
                AZ_TEMP = AZ;
                AZ = S;
            } while (Math.abs(S - AZ_TEMP) >= TOLERANCE_2);

            final double AZ1 = (dlon < 0.0) ? 2.0 * Math.PI - S : S;
            forwardAzimuth = new Angle(AZ1);
            // not sure what AZ2 is or why it is not used (could it be reverseAzimuth?)
            //final double AZ2 = 2.0 * Math.PI - AZ1;
            S = Math.cos(AZ1);

            // Equatorial - geodesic(S-s) SMS
            final double U2 = ESQP * S * S;
            final double U4 = U2 * U2;
            final double U6 = U4 * U2;
            final double U8 = U6 * U2;
            final double BO = 1.0 +
                    0.25 * U2 +
                    0.046875 * U4 +
                    0.01953125 * U6 +
                    -0.01068115234375 * U8;
            S = Math.sin(AZ1);
            final double SMS = semiMajorAxis * Math.PI * (1.0 - f * Math.abs(S) * AO - BO * fo);
            distanceInMeters = semiMajorAxis * ss - SMS;
            return;
        }

        // the reduced latitudes
        final double u1 = Math.atan(fo * Math.sin(lat1) / Math.cos(lat1));
        final double u2 = Math.atan(fo * Math.sin(lat2) / Math.cos(lat2));
        final double su1 = Math.sin(u1);
        final double cu1 = Math.cos(u1);
        final double su2 = Math.sin(u2);
        final double cu2 = Math.cos(u2);
        double xy, w, q2, q4, q6, r2, r3, sig, ssig, slon, clon, sinalf, ab = dlon;
        int kcount = 0;
        do {
            if (++kcount > 8) {
                String pts = point1.toString() + " to " + point2.toString();
                throw new ArithmeticException(NOCONVERGE + pts);
            }
            clon = Math.cos(ab);
            slon = Math.sin(ab);
            final double csig = su1 * su2 + cu1 * cu2 * clon;
            ssig = Math.sqrt(slon * cu2 * slon * cu2 + (su2 * cu1 - su1 * cu2 * clon) * (su2 * cu1 - su1 * cu2 * clon));
            sig = Math.atan2(ssig, csig);
            sinalf = cu1 * cu2 * slon / ssig;
            w = (1.0 - sinalf * sinalf);
            final double t4 = w * w;
            final double t6 = w * t4;

            // the coefficents of type a
            final double ao = f + a01 * w + a02 * t4 + a03 * t6;
            final double a2 = a21 * w + a22 * t4 + a23 * t6;
            final double a4 = a42 * t4 + a43 * t6;
            final double a6 = a63 * t6;

            // the multiple angle functions
            double qo = 0.0;
            if (w > TOLERANCE_0) {
                qo = -2.0 * su1 * su2 / w;
            }
            q2 = csig + qo;
            q4 = 2.0 * q2 * q2 - 1.0;
            q6 = q2 * (4.0 * q2 * q2 - 3.0);
            r2 = 2.0 * ssig * csig;
            r3 = ssig * (3.0 - 4.0 * ssig * ssig);

            // the longitude difference
            final double s = sinalf * (ao * sig + a2 * ssig * q2 + a4 * r2 * q4 + a6 * r3 * q6);
            double xz = dlon + s;
            xy = Math.abs(xz - ab);
            ab = dlon + s;
        } while (xy >= TOLERANCE_1);

        final double z = ESQP * w;
        final double bo = 1.0 + z * (1.0 / 4.0 + z * (-3.0 / 64.0 + z * (5.0 / 256.0 - z * (175.0 / 16384.0))));
        final double b2 = z * (-1.0 / 4.0 + z * (1.0 / 16.0 + z * (-15.0 / 512.0 + z * (35.0 / 2048.0))));
        final double b4 = z * z * (-1.0 / 128.0 + z * (3.0 / 512.0 - z * (35.0 / 8192.0)));
        final double b6 = z * z * z * (-1.0 / 1536.0 + z * (5.0 / 6144.0));

        // The distance in ellispoid axis units.
        distanceInMeters = semiMinorAxis * (bo * sig + b2 * ssig * q2 + b4 * r2 * q4 + b6 * r3 * q6);
        double az1 = (dlon < 0) ? Math.PI * (3 / 2) : Math.PI / 2;

        // now compute the az1 & az2 for latitudes not on the equator
        if ((Math.abs(su1) >= TOLERANCE_0) || (Math.abs(su2) >= TOLERANCE_0)) {
            final double tana1 = slon * cu2 / (su2 * cu1 - clon * su1 * cu2);
            final double sina1 = sinalf / cu1;

            // azimuths from north, longitudes positive east
            az1 = Math.atan2(sina1, sina1 / tana1);
        }
        forwardAzimuth = new Angle(az1);
    }

    private void calcPoint2() {
        final double azimuth = forwardAzimuth.inRadians();

        final double lat1 = point1.getLatitude().inRadians();
        final double lon1 = point1.getLongitude().inRadians();
        final double lat2;
        final double lon2;

        /*
        * Solution of the geodetic direct problem after T.Vincenty.
        * Modified Rainsford's method with Helmert's elliptical terms.
        * Effective in any azimuth and at any distance short of antipodal.
        *
        * Latitudes and longitudes in radians positive North and East.
        * Forward azimuths at both points returned in radians from North.
        *
        * Programmed for CDC-6600 by <NAME> NGS ROCKVILLE MD 18FEB75
        * Modified for IBM SYSTEM 360 by <NAME> NGS ROCKVILLE MD 7507
        * Ported from Fortran to Java by <NAME>.
        *
        * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/forward.for
        *         subroutine DIRECT1
        */
        double TU = fo * Math.sin(lat1) / Math.cos(lat1);
        double SF = Math.sin(azimuth);
        double CF = Math.cos(azimuth);
        double BAZ = (CF != 0) ? Math.atan2(TU, CF) * 2.0 : 0;
        double CU = 1 / Math.sqrt(TU * TU + 1.0);
        double SU = TU * CU;
        double SA = CU * SF;
        double C2A = 1.0 - SA * SA;
        double X = Math.sqrt((1.0 / fo / fo - 1) * C2A + 1.0) + 1.0;
        X = (X - 2.0) / X;
        double C = 1.0 - X;
        C = (X * X / 4.0 + 1.0) / C;
        double D = (0.375 * X * X - 1.0) * X;
        TU = distanceInMeters / fo / semiMajorAxis / C;
        double Y = TU;
        double SY, CY, CZ, E;
        do {
            SY = Math.sin(Y);
            CY = Math.cos(Y);
            CZ = Math.cos(BAZ + Y);
            E = CZ * CZ * 2.0 - 1.0;
            C = Y;
            X = E * CY;
            Y = E + E - 1.0;
            Y = (((SY * SY * 4.0 - 3.0) * Y * CZ * D / 6.0 + X) * D / 4.0 - CZ) * SY * D + TU;
        } while (Math.abs(Y - C) > TOLERANCE_1);
        BAZ = CU * CY * CF - SU * SY;
        C = fo * Math.sqrt(SA * SA + BAZ * BAZ);
        D = SU * CY + CU * SY * CF;
        lat2 = Math.atan2(D, C);
        C = CU * CY - SU * SY * CF;
        X = Math.atan2(SY * SF, C);
        C = ((-3.0 * C2A + 4.0) * f + 4.0) * C2A * f / 16.0;
        D = ((E * CY * C + CZ) * SY * C + Y) * SA;
        lon2 = lon1 + X - (1.0 - C) * D * f;
        point2 = new Geodetic2DPoint(new Longitude(lon2), new Latitude(lat2));
    }

    // Valid range checker for distance in meters
    private void validateDistance(double distanceInMeters) {
        if (distanceInMeters < 0.0 || distanceInMeters > maxOrthodromicDistance)
            throw new IllegalArgumentException("Distance is out of legal range (0 .. " +
                    maxOrthodromicDistance + ")");
    }

    /**
     * This method returns the fixed geodetic point 1 of this Geodetic2DArc object.
     *
     * @return Geodetic2DPoint point 1 (fixed)
     */
	@NonNull
    public Geodetic2DPoint getPoint1() {
        return point1;
    }

    /**
     * This method returns the re-settable geodetic point 2 of this Geodetic2DArc object.
     *
     * @return Geodetic2DPoint point 1 (re-settable)
     */
	@NonNull
    public Geodetic2DPoint getPoint2() {
        return point2;
    }

    /**
     * This method is used to re-set the second geodetic point of this Geodetic2DArc
     * object, resulting in the recalculation of forward azimuth and distance.
     *
     * @param point2 second Geodetic2DPoint (re-settable)
	 * @throws NullPointerException if point2 is null
     */
    public void setPoint2(Geodetic2DPoint point2) {
        this.point2 = point2;
        calcDistanceAndAzimuth();
    }

    /**
     * This method returns the Angle of forward azimuth from point 1 and point 2.
     *
     * @return Angle from North of forward azimuth (shortest distance direction from point1)
     */
	@NonNull
    public Angle getForwardAzimuth() {
        return forwardAzimuth;
    }

    /**
     * This method is used to set the Angle of forward azimuth from point 1, resulting in
     * point 2 to be recalculated at the current distance in meters.
     *
     * @param forwardAzimuth Angle from North of forward azimuth
	 * @throws NullPointerException if forwardAzimuth is null
     */
    public void setForwardAzimuth(Angle forwardAzimuth) {
        this.forwardAzimuth = forwardAzimuth;
        calcPoint2();
    }

    /**
     * This method returns the shortest distance in meters between the geodetic points
     * 1 and 2 of this Geodetic2DArc object.
     *
     * @return double shortest distance in meters on Ellipsoid between points 1 and 2
     */
    public double getDistanceInMeters() {
        return distanceInMeters;
    }

    /**
     * This method is used to set a new distance in meters for this Geodetic2DArc object,
     * resulting in the calculation of a new second point.
     *
     * @param distanceInMeters double distance from point 1 to new point 2
     * @throws IllegalArgumentException if distance is out of legal range
     */
    public void setDistanceInMeters(double distanceInMeters) {
        validateDistance(distanceInMeters);
        this.distanceInMeters = distanceInMeters;
        calcPoint2();
    }

    /**
     * This method is used to simultaneously set a new distance in meters and forward
     * azimuth for this Geodetic2DArc object, resulting in the calculation of a new second
     * point.
     *
     * @param distanceInMeters double distance from point 1 to new point 2
     * @param forwardAzimuth   Angle from North of forward azimuth
     * @throws IllegalArgumentException if distance is out of legal range
	 * @throws NullPointerException if forwardAzimuth is null
     */
    public void setDistanceAndAzimuth(double distanceInMeters, Angle forwardAzimuth)
            throws IllegalArgumentException {
        validateDistance(distanceInMeters);
        this.distanceInMeters = distanceInMeters;
        this.forwardAzimuth = forwardAzimuth;
        calcPoint2();
    }

    /**
     * Compares this Geodetic2DArc with the specified Geodetic2DArc for order. Returns a
     * negative integer, zero, or a positive integer as this arc's length is less than,
     * equal to, or greater than the length of the specified arc.
     *
     * @param that Geodetic2DArc to compare lengths with this one
     * @return a negative integer, zero, or a positive integer as this object is less than,
     *         equal to, or greater than the specified object.
	 * @throws NullPointerException if that is null
     */
    public int compareTo(Geodetic2DArc that) {
        Double thisLen = this.getDistanceInMeters();
        Double thatLen = that.getDistanceInMeters();
        return thisLen.compareTo(thatLen);
    }

    /**
     * This method returns a hash code for this Geodetic2DArc object.
     *
     * @return a hash code value for this object.
     */
    @Override
    public int hashCode() {
        return 31 * point1.hashCode() + point2.hashCode();
    }

    /**
     * This method is used to compare this arc to another Geodetic2DArc object.
     *
     * @param arc Geodetic2DArc to compare to this one
     * @return boolean indicating whether this arc is equal to the specified arc
     */
    public boolean equals(Geodetic2DArc arc) {
        return arc != null && point1.equals(arc.getPoint1()) && point2.equals(arc.getPoint2());
    }

    /**
     * This method is used to compare this arc to another Geodetic2DArc object.
     *
     * @param that Geodetic2DArc to compare to this one
     * @return boolean indicating whether this arc is equal to the specified object
     */
    @Override
    public boolean equals(Object that) {
        return that instanceof Geodetic2DArc && equals((Geodetic2DArc) that);
    }

    /**
     * The toString method formats the arc for printing.
     *
     * @return String representation of points, distance, and heading
     */
    @Override
    public String toString() {
        return "(" + point2 + " is " + distanceInMeters + "m from " +
                point1 + " at heading " + forwardAzimuth + ")";
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/PolarStereographic.java<|end_filename|>
/****************************************************************************************
 *  PolarStereographic.java
 *
 *  Created: Jan 18, 2007
 *
 *  @author <NAME> (Re-written in Java based closely on GeoTrans 2.4 C code)
 *
 *  POLAR STEREOGRAPHIC originated from :
 *      U.S. Army Topographic Engineering Center
 *      Geospatial Information Division
 *      7701 Telegraph Road
 *      Alexandria, VA  22310-3864
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.io.Serializable;

/**
 * The PolarStereographic class implements a map projection of the same name, and is used by
 * the Universal Polar Stereographic (UPS) coordinate encoding system for points near the
 * North or South Poles.<p/>
 *
 * TODO: implement equals() and hashCode() methods.
 */
public class PolarStereographic implements Serializable {
	private static final long serialVersionUID = 1L;
	
	// Constants
    private static final double PI_OVER_2 = Math.PI / 2.0;
    private static final double PI_OVER_4 = Math.PI / 4.0;
    private static final double EPSILON = 1.0e-2;

    // Instance variables
    @NonNull private Ellipsoid ellipsoid = Ellipsoid.getInstance("WGS 84");

    // Ellipsoid Parameters, default to WGS 84
    private double Polar_a = 6378137.0;                 // Semi-major axis of ellipsoid in meters
    private double Polar_e = 0.08181919084262188000;    // Eccentricity of ellipsoid
    private double Polar_eh = .040909595421311;         // Eccentricity over 2.0 (e half)
    private double Polar_2a = 12756274.0;               // 2.0 * Polar_a Semi-major axis

    private double Polar_tc = 1.0;
    private double Polar_e4 = 1.0033565552493;
    private double Polar_a_mc = 6378137.0;              // Polar_a * mc

    // Polar Stereographic projection Parameters
    private double Polar_Origin_Lat = PI_OVER_2;        // Latitude of true scale (init at N Pole)
    private final boolean Southern_Hemisphere;          // Depends on sign of latitude of true scale at construction
    private boolean True_Scale_At_Pole = true;          // Calculation method differs when true scale is at pole
    private double Polar_Origin_Long = 0.0;             // Longitude down from pole (northings decrease along)

    // Maximum variance for easting and northing values for WGS 84.
    private double Polar_Delta_Easting = 12713601.0;
    private double Polar_Delta_Northing = 12713601.0;

    // Utiltiy Function
    private double POLAR_POW(double esSin) {
        return (Math.pow((1.0 - esSin) / (1.0 + esSin), Polar_eh));
    }

    // Initialization method to update derivative parameters when ellipsoid or origin has changed
    private void init() {
        True_Scale_At_Pole = Math.abs(Math.abs(Polar_Origin_Lat) - PI_OVER_2) <= 1.0e-10;
        if (True_Scale_At_Pole) {
            double ope = 1.0 + Polar_e;    // One Plus Eccentricity
            double ome = 1.0 - Polar_e;    // One Minus Eccentricity
            Polar_e4 = Math.sqrt(Math.pow(ope, ope) * Math.pow(ome, ome));
        } else {
            double esSin = Polar_e * Math.sin(Polar_Origin_Lat);
            double mc = Math.cos(Polar_Origin_Lat) / Math.sqrt(1.0 - (esSin * esSin));
            Polar_a_mc = Polar_a * mc;
            Polar_tc = Math.tan(PI_OVER_4 - (Polar_Origin_Lat / 2.0)) / POLAR_POW(esSin);
        }
        Topocentric2DPoint en =
                toPolarStereographic(new Longitude(Polar_Origin_Long), new Latitude(0.0));
        Polar_Delta_Northing = Math.abs(en.getNorthing()) + EPSILON;
        Polar_Delta_Easting = Polar_Delta_Northing;
    }

    /**
     * This constructor makes an instance of PolarStereographic (initialized with WGS84 ellipsoid,
     * origin at North Pole Prime Meridian).  Ellipsoid can be changed with call to setEllipsoid.
     *
     * @param latOfTrueScale Latitude of True Scale (sign determines hemisphere of projection)
     * @throws NullPointerException if latOfTrueScale is null
     */
    public PolarStereographic(Latitude latOfTrueScale) {
        Polar_Origin_Lat = latOfTrueScale.inRadians;
        Southern_Hemisphere = (latOfTrueScale.inRadians < 0.0);
        if (Southern_Hemisphere) {
            Polar_Origin_Lat = -Polar_Origin_Lat;
            Polar_Origin_Long = -Polar_Origin_Long;
        }
        init();
    }

    /**
     * This method is used to set the Ellipsoid model parameters.  This is only needed if
     * the default WGS84 ellipsoid is not being used, since WGS84 values are the default.
     *
     * @param ellip Ellipsoid object to use
     * @throws NullPointerException if ellip is null
     */
    public void setEllipsoid(Ellipsoid ellip) {
        // Only update ellipsoid and other parameters if ellip has changed
        if (ellipsoid != ellip) {
            ellipsoid = ellip;

            Polar_a = ellip.getEquatorialRadius();
            Polar_e = ellip.getEccentricity();

            Polar_eh = Polar_e / 2.0;
            Polar_2a = Polar_a * 2.0;
            init();
        }
    }

    /**
     * This method will return the currently set Ellipsoid for this PolarStereographic instance.
     *
     * @return ellipsoid model of the earth currently being used by this PolarStereographic object
     */
    @NonNull
    public Ellipsoid getEllipsoid() {
        return ellipsoid;
    }

    /**
     * This method converts geodetic (latitude and longitude) coordinates
     * to Polar Stereographic projection (easting and northing) coordinates, according
     * to the current ellipsoid and Polar Stereographic projection coordinates.
     *
     * @param lon Longitude to project
     * @param lat Latitude to project
     * @return Topocentric2DPoint object containing Easting (x coord) and Northing (y coord) in meters
     * @throws IllegalArgumentException error if latitude is in different hemisphere from origin
     */
    @NonNull
    public Topocentric2DPoint toPolarStereographic(Longitude lon, Latitude lat)
            throws IllegalArgumentException {

        double latRad = lat.inRadians;
        double lonRad = lon.inRadians;
        if ((latRad < 0.0 && !Southern_Hemisphere) || (latRad > 0.0 && Southern_Hemisphere)) {
            throw new IllegalArgumentException("Latitude in different hemisphere " +
                    "from Polar Stereographic origin");
        }
        double easting = 0.0;
        double northing = 0.0;
        // Don't project if we're sufficiently near a pole
        if (Math.abs(Math.abs(latRad) - PI_OVER_2) >= 1.0e-10) {
            if (Southern_Hemisphere) {
                lonRad *= -1.0;
                latRad *= -1.0;
            }

            double t = Math.tan(PI_OVER_4 - (latRad / 2.0)) / POLAR_POW(Polar_e * Math.sin(latRad));
            // Compute differently if latitude of true scale is sufficiently away from pole
            double rho = (True_Scale_At_Pole) ? (Polar_2a * t / Polar_e4) : (Polar_a_mc * t / Polar_tc);
            double dlam = Angle.normalize(lonRad - Polar_Origin_Long);
            easting = rho * Math.sin(dlam);
            northing = -rho * Math.cos(dlam);

            if (Southern_Hemisphere) {
                easting *= -1.0;
                northing *= -1.0;
            }
        }
        return new Topocentric2DPoint(easting, northing);
    }

    /**
     * This method converts Polar Stereographic projection (easting and northing)
     * coordinates to geodetic (latitude and longitude) coordinates, according to the current
     * ellipsoid and Polar Stereographic projection parameters.
     *
     * @param easting  X in meters
     * @param northing Y in meters
     * @return Geodetic2DPoint object containing un-projected Longitude and Latitude
     * @throws IllegalArgumentException error if Easting or Northing are out of legal range
     */
    @NonNull
    public Geodetic2DPoint toGeodetic(double easting, double northing) {
        double lonRad;
        double latRad;

        if ((easting < -Polar_Delta_Easting) || (Polar_Delta_Easting < easting)) {
            throw new IllegalArgumentException("Polar Stereographic Easting value is out of range");
        }
        if ((northing < -Polar_Delta_Northing) || (Polar_Delta_Northing < northing)) {
            throw new IllegalArgumentException("Polar Stereographic Northing value is out of range");

        }
        double temp = Math.sqrt((easting * easting) + (northing * northing));
        if ((temp < -Polar_Delta_Easting) || (Polar_Delta_Easting < temp) ||
                (temp < -Polar_Delta_Northing) || (Polar_Delta_Northing < temp)) {
            throw new IllegalArgumentException("Polar Stereographic point is outside of projection area");
        }

        if ((northing == 0.0) && (easting == 0.0)) {
            latRad = PI_OVER_2;
            lonRad = Polar_Origin_Long;
        } else {
            if (Southern_Hemisphere) {
                northing *= -1.0;
                easting *= -1.0;
            }
            double rho = Math.sqrt((easting * easting) + (northing * northing));
            // Compute differently if latitude of true scale is sufficiently away from pole
            double t = rho * ((True_Scale_At_Pole) ? (Polar_e4 / Polar_2a) : (Polar_tc / Polar_a_mc));
            double Phi = PI_OVER_2 - (2.0 * Math.atan(t));
            double tempPhi = 0.0;
            while (Math.abs(Phi - tempPhi) > 1.0e-10) {
                tempPhi = Phi;
                Phi = PI_OVER_2 - (2.0 * Math.atan(t * POLAR_POW(Polar_e * Math.sin(Phi))));
            }
            latRad = Phi;
            lonRad = Angle.normalize(Polar_Origin_Long + Math.atan2(easting, -northing));

            // force distorted values to 90, -90 degrees
            if (latRad > PI_OVER_2) latRad = PI_OVER_2;
            else if (latRad < -PI_OVER_2) latRad = -PI_OVER_2;
            // force distorted values to 180, -180 degrees
            if (lonRad > Math.PI) lonRad = Math.PI;
            else if (lonRad < -Math.PI) lonRad = -Math.PI;
        }
        if (Southern_Hemisphere) {
            latRad *= -1.0;
            lonRad *= -1.0;
        }

        return new Geodetic2DPoint(new Longitude(lonRad), new Latitude(latRad));
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Ellipsoid.java<|end_filename|>
/****************************************************************************************
 *  Ellipsoid.java
 *
 *  Created: Dec 29, 2006
 *
 *  @author <NAME> (derived from API found at
 *    "http://www.inf.ufsc.br/grafos/represen/java-or/drasys/or/geom/geo/Ellipsoid.html")
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.io.Serializable;
import java.util.Map;
import java.util.HashMap;
import java.util.Set;

/**
 * The Ellipsoid class maintains an object pool of Singleton Ellipsoid objects (used as
 * models of a slightly flattened spherical earth), indexed by String names.  The user
 * can create additional named Ellipsoids if needed, or retrieve one of the statically
 * defined ones that are common or in historic use by geographic information and mapping
 * systems. There is an accessor method to get the list of available cached Ellipsoid
 * names, and there are accessor methods for the commonly used mathematical parameters
 * of an Ellipsoid model. <p/>
 * 
 * TODO: implement equals() and hashCode() methods.
 */
public class Ellipsoid implements Serializable {
	private static final long serialVersionUID = 1L;
	
	private final String name;
    private final double a;       // Radius in m of circle on equitorial plane (parameter "a")
    private final double b;       // Distance from ellipsoid center to a pole (parameter "b")
    private final double f;       // Flattening ratio f = (a - b) / a
    private final double es;      // Eccentricity of ellipsoid, squared
    private final double ebs;     // Second Eccentricity of ellipsoid, squared

    // Cache of frequently used singleton Ellipsoids (comment out if not needed by app)
    private static final Map<String, Ellipsoid> namedEllipsoids;

    static {
        namedEllipsoids = new HashMap<String, Ellipsoid>();

        new Ellipsoid("WGS 84", 6378137.000, 1.0 / 298.25722356299722);
        new Ellipsoid("WGS 72", 6378135.000, 1.0 / 298.26);
        new Ellipsoid("WGS 66", 6378145.000, 1.0 / 298.25);
        new Ellipsoid("WGS 60", 6378165.000, 1.0 / 298.30);
        new Ellipsoid("South American 1969", 6378160.000, 1.0 / 298.25);
        new Ellipsoid("Krassovsky 1940", 6378245.000, 1.0 / 298.30);
        new Ellipsoid("Hough 1956", 6378270.000, 1.0 / 297.00);
        new Ellipsoid("GRS 1980", 6378137.000, 1.0 / 298.257222101);
        new Ellipsoid("GRS 1975", 6378140.000, 1.0 / 298.257);
        new Ellipsoid("GRS 1967", 6378160.000, 1.0 / 298.247167427);
        new Ellipsoid("Fischer 1968", 6378150.000, 1.0 / 298.30);
        new Ellipsoid("Fischer 1960", 6378166.000, 1.0 / 298.30);
        new Ellipsoid("Everest 1830", 6377276.345, 1.0 / 300.8017);
        new Ellipsoid("Clarke 1880", 6378249.145, 1.0 / 293.465);
        new Ellipsoid("Clarke 1866", 6378206.400, 1.0 / 294.9786982);
        new Ellipsoid("Bessel 1841", 6377397.155, 1.0 / 299.1528128);
        new Ellipsoid("Bessel 1841 Namibia", 6377483.865, 1.0 / 299.1528128);
        new Ellipsoid("Airy 1830", 6377563.396, 1.0 / 299.3249646);
    }

    /**
     * This is the constructor for a custom ellipsoid to be created and cached in the
     * object pool. If the name has not been used, the new Ellipsoid will be cached as
     * a singleton for future access. If already defined, an IllegalArgumentException
     * is thrown. Commonly used Ellipsoids are defined in a static block of this class
     * at class load time, so the constructor is typically not needed by client programs.
     * Use the getInstance method to retrieve by name instead.
     *
     * @param name             String name used to get this singleton Ellipsoid instance
     * @param equatorialRadius radius in meters of circle on equitorial plane, a
     * @param flattening       ratio of axes difference to equatorial radius, f = (a - b) / a
     * @throws IllegalArgumentException error if named ellipoid has been previously defined
     */
    public Ellipsoid(String name, double equatorialRadius, double flattening)
            throws IllegalArgumentException {
        // User can't redefine initial named ellipsoids using standard names
        if (namedEllipsoids.get(name) != null)
            throw new IllegalArgumentException("Ellipsoid " + name + " is already defined");

        // Do error checking on arguments
        if (equatorialRadius <= 0.0)
            throw new IllegalArgumentException("Ellipsoid semi-major axis (equatorialRadius)" +
                    " must be greater than zero");

        double inv_f = 1.0 / flattening;
        if ((inv_f < 250.0) || (350.0 < inv_f))
            throw new IllegalArgumentException("Ellipsoid's inverse flattening (1/f)" +
                    " must be between 250 and 350");

        this.name = name;
        a = equatorialRadius;
        f = flattening;
        b = equatorialRadius * (1.0 - flattening);

        es = (2.0 * f) - (f * f);
        ebs = (1.0 / (1.0 - es)) - 1.0;

        namedEllipsoids.put(name, this);
    }

    /**
     * This method returns a Set of String names of defined Ellipsoids available in the
     * object pool.
     *
     * @return the Set of String names for currently defined Ellipsoid singletons
     */
    @NonNull
    public static Set<String> getEllipsoidNames() {
        return namedEllipsoids.keySet();
    }

    /**
     * This method will retrieve the named singleton Ellipsoid, if defined.
     *
     * @param name String used to identify the Ellipsoid (e.g. "WGS 84")
     * @return Ellipsoid singleton with given name, retrieved from cache
     */
    public static Ellipsoid getInstance(String name) {
        return namedEllipsoids.get(name);
    }

    /**
     * This accessor method returns the name of this Ellipsoid.
     *
     * @return name String used to identify this ellipsoid
     */
    public String getName() {
        return name;
    }

    /**
     * This accessor method returns the equatorialRadius, ellipsoid parameter 'a'
     * (semi-major axis).
     *
     * @return radius in meters of circle on equitorial plane
     */
    public double getEquatorialRadius() {
        return a;
    }

    /**
     * This accessor method returns the polarRadius, ellipsoid parameter 'b'
     * (semi-minor axis).
     *
     * @return distance in meters from center to pole
     */
    public double getPolarRadius() {
        return b;
    }

    /**
     * This method returns the radius of curvature of this Ellipsoid, in meters,
     * at the given Latitude.  It is also called the radius of curvature in the prime
     * vertical.
     *
     * @param lat the Latitude for which the ellipsoid's radius is desired
     * @return the radius of the ellipsoid at the given Latitude, in meters.
     */
    public double getRadius(Latitude lat) {
        double sinLat = Math.sin(lat.inRadians());
        return a / Math.sqrt(1.0 - es * sinLat * sinLat);
    }

    /**
     * This accessor method returns the ellipsoid flattening parameter,
     * 'f' = (a - b) / a.
     *
     * @return flattening ratio of axes difference to equatorial radius
     */
    public double getFlattening() {
        return f;
    }

    /**
     * This accessor method returns the ellipsoid eccentricity parameter,
     * 'e' = sqrt(1 - (b^2 / a^2)).
     *
     * @return eccentricity of the ellipsoid
     */
    public double getEccentricity() {
        return Math.sqrt(es);
    }


    /**
     * This accessor method returns the ellipsoid eccentricity value squared.
     *
     * @return eccentricity of the ellipsoid, squared
     */
    public double getEccentricitySquared() {
        return es;
    }

    /**
     * This accessor method returns the ellipsoid second eccentricity value squared.
     *
     * @return second eccentricity of the ellipsoid, squared
     */
    public double get2ndEccentricitySquared() {
        return ebs;
    }

    /**
     * This method calculates the orthodromic distance in meters between two points, which
     * is the shortest distance between two geodetic points on the surface of this ellipsoid.
     *
     * @param p1 geodetic point 1
     * @param p2 geodetic point 2
     * @return distance in meters
     */
    public double orthodromicDistance(Geodetic2DPoint p1, Geodetic2DPoint p2) {
        Geodetic2DArc arc = new Geodetic2DArc(this, p1, p2);
        return arc.getDistanceInMeters();

        /*
        // Alternate method (appears slightly faster) - don't know which method is most accurate
        // Adapted from <NAME>: from Atronomical Algorithms, <NAME>, 2nd edition.
        double lon1 = p1.getLongitude().inRadians();
        double lat1 = p1.getLatitude().inRadians();
        double lon2 = p2.getLongitude().inRadians();
        double lat2 = p2.getLatitude().inRadians();

        double f = (lat1 + lat2) / 2.0;
        double g = (lat1 - lat2) / 2.0;
        double l = (lon1 - lon2) / 2.0;

        double sf = Math.sin(f);
        double sg = Math.sin(g);
        double sl = Math.sin(l);

        double s, c, w, r, d, h1, h2;
        sg = sg * sg;
        sl = sl * sl;
        sf = sf * sf;
        s = sg * (1.0 - sl) + (1.0 - sf) * sl;
        c = (1.0 - sg) * (1.0 - sl) + sf * sl;

        if ((c == 0) || (s == 0)) {
            return 0.0;
        } else {
            w = Math.atan(Math.sqrt(s / c));
            r = Math.sqrt(s * c) / w;
            d = 2.0 * w * a;
            h1 = (3.0 * r - 1.0) / 2.0 / c;
            h2 = (3.0 * r + 1.0) / 2.0 / s;
            return d * (1.0 + this.f * (h1 * sf * (1.0 - sg) - h2 * (1.0 - sf) * sg));
        }
        */
    }

    /**
     * This method returns the name of this Ellipsoid.
     *
     * @return name String used to identify this Ellipsoid in the object pool.
     */
    public String toString() {
        return getName();
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/GeoPoint.java<|end_filename|>
/****************************************************************************************
 *  GeoPoint.java
 *
 *  Created: Mar 23, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including
 *  the warranty of non-infringement and the implied warranties of merchantability and
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any
 *  damages suffered by you as a result of using the Program.  In no event will the
 *  Copyright owner be liable for any special, indirect or consequential damages or
 *  lost profits even if the Copyright owner has been advised of the possibility of
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

/**
 * A GeoPoint is defined by its interface as any object that can be accurately
 * converted to a Geodetic3DPoint point given a FrameOfReference. GeoPoints can
 * therefore be compared for nearness on the surface of an Ellipsoid model of
 * the Earth, or nearness in space as measured by Euclidean distance. This
 * interface enables the three different point representations used in this
 * geodesy package (Geodetic, Geocentric, and Topocentric) to be abstracted
 * into a single logical type.
 */
public interface GeoPoint {
    public Geodetic3DPoint toGeodetic3D(FrameOfReference fRef);
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestGeodetic3DBounds.java<|end_filename|>
/****************************************************************************************
 *  TestGeodetic3DBounds.java
 *
 *  Created: Mar 29, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantibility and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy.test;

import junit.framework.TestCase;
import junit.framework.TestSuite;
import junit.textui.TestRunner;

import org.opensextant.geodesy.*;

import java.util.Random;

public class TestGeodetic3DBounds extends TestCase {

	private static final double MIN_ELEV = 100;
	private static final double MAX_ELEV = 2500;

	private final Geodetic3DPoint west = new Geodetic3DPoint(new Longitude(161, 54, 44),
			new Latitude(-85, 41, 54), MIN_ELEV);
	
	private final Geodetic3DPoint east = new Geodetic3DPoint(new Longitude(99, 8, 8),
            new Latitude(79, 39, 57), MAX_ELEV);

    public void testBBox() {
        // bbox1: (161 deg 54' 44" E, 85 deg 41' 54" S) .. (99 deg 8' 8" E, 79 deg 39' 57" N)
        // bbox2: (91 deg 4' 4" E, 89 deg 57' 12" S) .. (0 deg 13' 54" W, 87 deg 50' 13" N)
        // together??: (161 deg 54' 44" E, 89 deg 57' 12" S) .. (99 deg 8' 8" E, 87 deg 50' 13" N)

		double minElev = MIN_ELEV;
		double maxElev = MAX_ELEV;

        Geodetic3DBounds bbox1 = new Geodetic3DBounds(west, east);

		assertTrue(bbox1.toString().indexOf("161") > 0);
        assertTrue(bbox1.contains(east));

        Geodetic3DPoint cPt = bbox1.getCenter();
        assertTrue(cPt.getElevation() >= minElev);
        assertTrue(cPt.getElevation() <= maxElev);

        // expand the volume of the bounding box
        minElev -= 10;
        maxElev += 10;
        Geodetic3DPoint west2 = new Geodetic3DPoint(new Longitude(91, 4, 4),
                new Latitude(-89, 57, 12), maxElev);
        Geodetic3DPoint east2 = new Geodetic3DPoint(new Longitude(-1, 13, 54),
                new Latitude(87, 50, 13), minElev);
        Geodetic3DBounds bbox2 = new Geodetic3DBounds(west2, east2);

        bbox1.include(bbox2);

        cPt = bbox1.getCenter();
        assertTrue(cPt.getElevation() >= minElev);
        assertTrue(cPt.getElevation() <= maxElev);

        assertTrue(bbox1.contains(bbox2));
        assertTrue(bbox1.contains(east2));

        Geodetic3DPoint outside = new Geodetic3DPoint(new Longitude(172, 54, 44),
                new Latitude(89, 41, 54), minElev);
        assertFalse(bbox1.contains(outside));

        Geodetic3DPoint ptTooHigh = new Geodetic3DPoint(cPt.getLongitude(), cPt.getLatitude(), maxElev * 2);
        assertFalse(bbox1.contains(ptTooHigh));

		Geodetic3DBounds readonlyCopy = new UnmodifiableGeodetic3DBounds(bbox1);
		assertEquals(bbox1, readonlyCopy);
		assertTrue(bbox1.equals((Geodetic2DBounds)readonlyCopy));
		try {
			readonlyCopy.grow(100);
			fail("readonly bounds expected to throw UnsupportedOperationException");
		} catch (UnsupportedOperationException e) {
			// expected
		}
		try {
			readonlyCopy.include(outside);
			fail("readonly bounds expected to throw UnsupportedOperationException");
		} catch (UnsupportedOperationException e) {
		}
		try {
			readonlyCopy.include(bbox1);
			fail("readonly bounds expected to throw UnsupportedOperationException");
		} catch (UnsupportedOperationException e) {
		}
		try {
			readonlyCopy.include((Geodetic2DPoint)outside);
			fail("readonly bounds expected to throw UnsupportedOperationException");
		} catch (UnsupportedOperationException e) {
		}
		try {
			readonlyCopy.setEastLon(east2.getLongitude());
			fail("readonly bounds expected to throw UnsupportedOperationException");
		} catch (UnsupportedOperationException e) {
			// expected
		}
		assertEquals(bbox1, readonlyCopy);
		assertEquals(bbox1.hashCode(), readonlyCopy.hashCode());
    }

	public void testContains() {
        // bbox: (161 deg 54' 44" E, 85 deg 41' 54" S) .. (99 deg 8' 8" E, 79 deg 39' 57" N)

        //double minElev = 100;
        //double maxElev = 2500;
        Geodetic3DBounds bbox = new Geodetic3DBounds(west);
		Geodetic2DPoint pt = west;
		assertEquals(bbox.maxElev, west.getElevation(), 1e-5);
		assertEquals(bbox.minElev, west.getElevation(), 1e-5);
		assertTrue(bbox.contains(pt)); // test as 3-d point
		pt = new Geodetic2DPoint(west.getLongitude(), west.getLatitude());
		assertTrue(bbox.contains(pt)); // tested as 2-d point in this context
		pt = new Geodetic3DPoint(west.getLongitude(), west.getLatitude(), MIN_ELEV + MAX_ELEV);
		assertFalse(bbox.contains(pt)); // tested as 3-d point with wrong elevation

		bbox.include((Geodetic2DPoint)east); // should include as 3d point
		assertTrue(bbox.contains((Geodetic2DPoint)east)); // test as 3-d point
		assertEquals(bbox.maxElev, MAX_ELEV, 1e-5);
		assertEquals(bbox.minElev, MIN_ELEV, 1e-5);

		Geodetic3DBounds bbox2 = new Geodetic3DBounds();
		final Geodetic2DBounds bounds2d = bbox; // handle as 2-d bbox
		bbox2.include(bounds2d); // should include as 3D bbox
		assertEquals(bbox.maxElev, MAX_ELEV, 1e-5);
		assertEquals(bbox.minElev, MIN_ELEV, 1e-5);
		assertTrue(bbox2.contains(bounds2d)); // should test as 3-d point

		bbox2 = new Geodetic3DBounds(bounds2d, 50.0, 50000.0);
		assertEquals(bbox2.maxElev, 50000.0, 1e-5);
		assertTrue(bbox2.contains(bbox));
	}

	public void test2dEquals() {
		Geodetic3DBounds bbox = new Geodetic3DBounds(east);
		bbox.maxElev = 0.0;
		bbox.minElev = 0.0;
		Geodetic2DBounds bounds2d = new Geodetic2DBounds(bbox);
		assertTrue(bbox.contains(bounds2d)); // should test as 2-d point
		assertEquals(bbox, bounds2d);   // should test as 2-d point

		bbox.include(new Geodetic2DPoint(west.getLongitude(), west.getLatitude()));
		assertTrue(bbox.contains(bounds2d));
	}

	public void testNullBBoxCompare() {
        Geodetic3DBounds bbox1 = new Geodetic3DBounds();
		Geodetic3DBounds bbox2 = null;
		assertFalse(bbox1.equals(bbox2));
	}

    public void testRandomBBox() {
        /*
        * sometimes this test fails if random bbox wraps world
        * so we run test multiple times where a single successful run
        * means the test passed. If fails 8 consecutive times then got a real problem.
        */
        AssertionError ex = null;
        for (int i = 1; i <= 8; i++) {
            try {
                realTestRandomBBox();
                return; // test successful
            } catch (WrappedAssertionException e) {
                System.out.println("*** warning: failed at testRandomBBox: " + i);
                AssertionError ae = e.getAssertionError();
                ae.printStackTrace(System.out);
                if (ex == null) ex = ae; // save first failed test result
            }
        }
        if (ex != null) {
            // this means we failed all attempts so we really failed
            throw ex;
        }
    }

    private void realTestRandomBBox() throws WrappedAssertionException {
        Random r = new Random();
        FrameOfReference f = new FrameOfReference();
        Geodetic3DPoint pt = TestGeoPoint.randomGeodetic3DPoint(r);
        Geodetic3DBounds bbox1 = new Geodetic3DBounds(pt.toGeodetic3D(f));
        for (int i = 0; i < 10; i++) {
            pt = TestGeoPoint.randomGeodetic3DPoint(r);
            bbox1.include(pt.toGeodetic3D(f));
        }
        assertTrue(bbox1.contains(pt));
        Geodetic3DBounds bbox2 = new Geodetic3DBounds(pt.toGeodetic3D(f));
        for (int i = 0; i < 10; i++) {
            pt = TestGeoPoint.randomGeodetic3DPoint(r);
            bbox2.include(pt.toGeodetic3D(f));
        }
        assertTrue(bbox2.contains(pt));
        bbox1.include(bbox2);
        assertTrue(bbox1.contains(bbox2));
        try {
            assertTrue("bbox contains point", bbox1.contains(pt)); // sometimes this test fails: see comment below
        } catch (AssertionError e) {
            // wrap this assertion failure so we can re-run the test again
            throw new WrappedAssertionException(e);
        }

        /*
         try {
            assertTrue(bbox1.contains(pt));
        } catch(Throwable e) {
            // once in a while there is an null assertion error here.
            e.printStackTrace();
            System.out.println(bbox1);
            System.out.println("center=" + bbox1.getCenter());
            System.out.println(pt);
            fail();
        }
        somehow the east/west lon get switched or we're wrapping the world
        example: west lon = -61 east lin = -81
            (61 deg 53' 58" W, 89 deg 37' 35" S) .. (81 deg 9' 34" W, 89 deg 33' 41" N) .. (1802.0m, 997780.0m)
            center=(108 deg 28' 14" E, 0 deg 1' 57" S) @ 499791m
            (64 deg 49' 50" W, 80 deg 39' 21" S) @ 853675m
         */
    }

    /**
     * Main method for running class tests.
     *
     * @param args standard command line arguments - ignored.
     */
    public static void main(String[] args) {
        TestSuite suite = new TestSuite(TestGeodetic3DBounds.class);
        new TestRunner().doRun(suite);
    }

    final static class WrappedAssertionException extends Exception {
        private final AssertionError ae;
        
        public WrappedAssertionException(AssertionError e) {
            super(e);
            this.ae = e;
        }

        public AssertionError getAssertionError() {
            return ae;
        }
    }
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestGeodetic3DPoint.java<|end_filename|>
package org.opensextant.geodesy.test;

import junit.framework.TestCase;
import junit.framework.TestSuite;
import junit.textui.TestRunner;

import org.opensextant.geodesy.*;

import java.util.Random;

public class TestGeodetic3DPoint extends TestCase {

    private final Random r = new Random();

    /**
     * This method is used to generate a random Geodetic3DPoint
     * that is near the surface of the earth (between -12,000 and +26,000 meters from the
     * surface of the Elliposoid earth model).
     *
     * @param r a Random object to use for selecting the coordinate components
     * @return a random Geodetic3DPoint suitable for testing conversions with
     */
    public static Geodetic3DPoint randomGeoPoint(Random r) {
        double maxBelow = 12000.0; // deeper than max ocean depth (Marianas Trench)
        double maxAbove = 26000.0; // higher than max altitude of SR-71 spy plane

        // Define a random Geodetic point
        double lonDeg = (r.nextDouble() * 360.0) - 180.0;  // lonDeg
        Longitude lon = new Longitude(lonDeg, Angle.DEGREES);
        double latDeg = (r.nextDouble() * 180.0) - 90.0;  // latDeg
        Latitude lat = new Latitude(latDeg, Angle.DEGREES);
        double elev = r.nextDouble() * ((r.nextBoolean()) ? maxAbove : -maxBelow);
        return new Geodetic3DPoint(lon, lat, elev);
    }

    /**
     * This method performs pair-wise format conversions to test for consistency
     * among the GeoPoint representations, using the default FrameOfReference.
     */
    public void testConversions() {
        FrameOfReference f = new FrameOfReference();
        for (int i = 0; i < 1000; i++) {
            GeoPoint a1 = randomGeoPoint(r);
            GeoPoint a2 = f.toTopocentric(a1);
            GeoPoint a3 = f.toGeocentric(a1);
            GeoPoint a4 = f.toGeodetic(a2);
            assertTrue(f.proximallyEquals(a1, a2));
            assertTrue(f.proximallyEquals(a1, a3));
            assertTrue(f.proximallyEquals(a1, a4));
        }
    }

	/**
     * This method performs format conversions with string constructors in degrees and dms
     */
    public void testRandomEquals() {
		for (int i = 0; i < 1000; i++) {
            Geodetic3DPoint a1 = randomGeoPoint(r);
			Geodetic3DPoint a2 = new Geodetic3DPoint(a1.getLongitude(), a1.getLatitude(), a1.getElevation());
			// note by making equals() work in round off errors (such as in this case) using Angle.equals() vs phi1=phi2 && lamb1==lamb2
			// but break contract in hashCode such that a.equals(b) -> true but hashCode(a) may not equal hashCode(b)
			assertEquals(a1, a2);
			assertEquals(a1.hashCode(), a2.hashCode());

			// for symmetric tests to work elevation must be non-zero
			final double elevation = a1.getElevation();
			if (elevation == 0.0 || Math.abs(elevation) < 1e-8) a1.setElevation(1234.5);

			// test symmetric equals tests a.equals(b) -> b.equals(a)
			Geodetic2DPoint pt2 = new Geodetic2DPoint(a1.getLongitude(), a1.getLatitude());
			assertFalse(pt2.equals(a1));
			assertFalse(a1.equals(pt2));
			a1.setElevation(0);
			assertEquals(pt2, a1); // pt2.equals(al) -> a1.equals(pt2)
			assertEquals(a1, pt2);
		}
    }

	private final Geodetic3DPoint a = new Geodetic3DPoint(new Longitude(-1), new Latitude(1), 31);
	private final Geodetic3DPoint b = new Geodetic3DPoint(new Longitude(-1), new Latitude(1), 31);
	private final Geodetic3DPoint c = new Geodetic3DPoint(new Longitude(-1.001), new Latitude(1), 31);
	private final Geodetic3DPoint d = new Geodetic3DPoint(new Longitude(-1), new Latitude(1), 31.001);

	public void testEquals() {
		// test equality with known geo-points
		assertEquals(a, b);
		assertFalse(a.equals(c));
		assertFalse(a.equals(d));
		assertFalse(c.equals(d));

		// test equality with any possible round off errors
		Geodetic3DPoint a2 = new Geodetic3DPoint(new Longitude(Math.toRadians(a.getLongitude().inDegrees())),
									new Latitude(Math.toRadians(a.getLatitude().inDegrees())), a.getElevation());
		assertEquals(a, a2);
		assertEquals(a.hashCode(), a2.hashCode());

		// approximate equals test for elevations up to 3 decimal places
		Geodetic3DPoint a3 = new Geodetic3DPoint(a.getLongitude(), a.getLatitude(), a.getElevation() + 10e-6);
		assertEquals(a, a3);
		assertEquals(a.hashCode(), a3.hashCode());
	}

	public void testNullCompare() {
		Geodetic3DPoint other = null;
		assertFalse(a.equals(other));
		Geodetic2DPoint p2 = null;
		assertFalse(a.equals(p2));
	}

    /**
     * Test 2d vs 3d points with 3d elavations at/close to 0
     */
    public void testMismatchedEquals() {
        Geodetic2DPoint p1 = makePoint(-81.9916466079043, 29.9420387052815, 5000.1);
		Geodetic2DPoint p2 = makePoint(-81.9916466079043, 29.9420387052815, 0.0);
        if (p1.equals(p2)) fail("different elevation but equals() == true");
        if (p2.equals(p1)) fail("different elevation but equals() == true");
		Geodetic2DPoint p3 = makePoint(-81.9916466079043, 29.9420387052815);
        assertEquals("2d with elev=0 and 3d point same lat/lon", p2, p3);
        assertEquals("2d with elev=0 and 3d point same lat/lon", p3, p2);
		Geodetic2DPoint p4 = makePoint(-81.9916466079043, 29.9420387052815, 1e-6);
        //assertEquals(p3, new Geodetic2DPoint(p4.getLongitude(), p4.getLatitude()));
        assertEquals("3d with elev=1e-4 and 3d point same lat/lon", p2, p4);
        assertEquals("3d with elev=1e-4 and 3d point same lat/lon", p4, p2);
        assertEquals("3d with elev=1e-4 and 3d point same lat/lon", p3, p4);
    }

	public void testHashCode() {
		assertEquals(a.hashCode(), b.hashCode());
		assertTrue(a.hashCode() != c.hashCode()); // this test isn't required by equals-hashCode contract but by how the hashCode is computed
	}

    public void testToString() {
        Geodetic3DPoint p = TestGeoPoint.randomGeodetic3DPoint(r);
        int prevLen = 0;

        String base = p.toString();
        String prefix = base.substring(0, base.indexOf(' '));      // e.g. "75 deg"
        char letter = base.charAt(base.indexOf(')') - 1);          // e.g. "S"

        // System.out.println("base=" + base + " pre=" + prefix + " letter=" + letter);

        for(int i=0; i < 6; i++) {
            /*
           length of string should increase two digits for each fractional # of digits in lat & lon
           0 (133 deg 47' 15" W, 73 deg 17' 33" S) @ 733484m
           1 (133 deg 47' 14.8" W, 73 deg 17' 32.7" S) @ 733484.0m
           2 (133 deg 47' 14.81" W, 73 deg 17' 32.67" S) @ 733484.00m
           3 (133 deg 47' 14.806" W, 73 deg 17' 32.675" S) @ 733484.000m
           4 (133 deg 47' 14.8060" W, 73 deg 17' 32.6749" S) @ 733484.0000m
           5 (133 deg 47' 14.80597" W, 73 deg 17' 32.67487" S) @ 733484.00000m
            ...
            */
            String s = p.toString(i);
            System.out.println(i + " " + s);
            int len = s.length();
            assertTrue(len >= prevLen + 2);
            // simple tests: 1) each should start with same degrees and 2) end with same hemisphere
            assertTrue(s.startsWith(prefix));
            assertEquals(letter, s.charAt(s.indexOf(')') - 1));
            prevLen = len;
        }
    }

    private static Geodetic2DPoint makePoint(double lon, double lat, double elev) {
		return new Geodetic3DPoint(
				new Longitude(lon, Angle.DEGREES),
				new Latitude(lat, Angle.DEGREES),
				elev);
	}

	private static Geodetic2DPoint makePoint(double lon, double lat) {
		return new Geodetic2DPoint(
				new Longitude(lon, Angle.DEGREES),
				new Latitude(lat, Angle.DEGREES));
	}

    /**
     * Main method for running class tests.
     *
     * @param args standard command line arguments - ignored.
     */
    public static void main(String[] args) {
        TestSuite suite = new TestSuite(TestGeodetic3DPoint.class);
        new TestRunner().doRun(suite);
    }
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestISO8601DateTime.java<|end_filename|>
/************************************************************************************
 * TestISO8601DateTime.java 12/22/10 11:08 AM psilvey
 *
 *  (C) Copyright MITRE Corporation 2010
 *
 * The program is provided "as is" without any warranty express or implied, including
 * the warranty of non-infringement and the implied warranties of merchantability and
 * fitness for a particular purpose.  The Copyright owner will not be liable for any
 * damages suffered by you as a result of using the Program. In no event will the
 * Copyright owner be liable for any special, indirect or consequential damages or
 * lost profits even if the Copyright owner has been advised of the possibility of
 * their occurrence.
 *
 ***********************************************************************************/
package org.opensextant.geodesy.test;

import junit.framework.TestCase;
import org.junit.Test;
import org.junit.runner.JUnitCore;
import org.opensextant.geodesy.ISO8601DateTimeInterval;
import org.opensextant.geodesy.ISO8601DateTimePoint;

import java.util.Random;

public class TestISO8601DateTime extends TestCase {
    private final static Random rand = new Random();

    private final static String[] partials = {
                "1970",
                "1970-01",
                "1970-01-01",
                "1970-01-01T00",
                "1970-01-01T00:00",
                "1970-01-01T00:00:00",
                "1970-01-01T00:00:00.000",
                "1970-01-01T00:00:00.000Z"
        };

    @Test
    public static void testTimePoint() throws Exception {
        ISO8601DateTimePoint t1;
        ISO8601DateTimePoint t2;

        for (int trials = 0; trials < 1000; trials++) {
            // Create a random point in time (within about 8000 years of the Epoch)
            long rTime = (long) (100.0 * 365.25 * 24.0 *
                    60.0 * 60.0 * 1000.0 * rand.nextDouble());
            t1 = new ISO8601DateTimePoint(rTime);
            String isoDateTime = t1.toString();
            t2 = new ISO8601DateTimePoint(isoDateTime);
            assertEquals(t1, t2);
            assertEquals(isoDateTime, t2.toString());
			assertEquals(t1.hashCode(), t2.hashCode());
        }

        ISO8601DateTimePoint lastTime = null;
        for (String isoStr : partials) {
            t1 = new ISO8601DateTimePoint(isoStr);
            // all times should be equivalent
            if (lastTime != null) assertEquals(lastTime, t1);
            lastTime = t1;
            System.out.println(isoStr + " means " + t1.toString());
        }
    }

	@Test
    public static void testDefaultTimeInterval() {
		final long time = System.currentTimeMillis();
		ISO8601DateTimeInterval t1 = new ISO8601DateTimeInterval();
		final long startTime = t1.getStartTimeInMillis();
		assertTrue(startTime >= time);
		assertEquals(startTime, t1.getEndTimeInMillis());
	}

	@Test
	public static void testCompareTo() throws Exception {
        ISO8601DateTimePoint t1 = new ISO8601DateTimePoint();
        ISO8601DateTimePoint t2 = new ISO8601DateTimePoint(t1.getStartTimeInMillis());
		assertEquals(0, t1.compareTo(t2));
		ISO8601DateTimePoint t3 = new ISO8601DateTimePoint(t1.getStartTimeInMillis() + 1000);
		assertEquals(-1, t1.compareTo(t3));
		ISO8601DateTimePoint t4 = new ISO8601DateTimePoint(t1.getStartTimeInMillis() - 1000);
		assertEquals(1, t1.compareTo(t4));
	}

    @Test
    public static void testTimeInterval() throws Exception {
        ISO8601DateTimeInterval t1;
        ISO8601DateTimeInterval t2;

        for (int trials = 0; trials < 1000; trials++) {
            // Create random points in time (within about 8000 years of the Epoch)
            long rt1 = (long) (100.0 * 365.25 * 24.0 *
                    60.0 * 60.0 * 1000.0 * rand.nextDouble());
            long rt2 = (long) (100.0 * 365.25 * 24.0 *
                    60.0 * 60.0 * 1000.0 * rand.nextDouble());
            // Make sure they are in the correct temporal order
            if (rt1 > rt2) {
                long temp = rt1;
                rt1 = rt2;
                rt2 = temp;
            }

            t1 = new ISO8601DateTimeInterval(rt1, rt2);
            String isoDateTime = t1.toString();
            t2 = new ISO8601DateTimeInterval(isoDateTime);
            assertEquals(t1, t2);
            assertEquals(isoDateTime, t2.toString());
			assertEquals(t1.hashCode(), t2.hashCode());
        }

		t1 = new ISO8601DateTimeInterval();
		ISO8601DateTimePoint pt = new ISO8601DateTimePoint();
		assertFalse(t1.equals(pt)); // t1.class != pt.class

        for (String isoStr : partials) {
            t1 = new ISO8601DateTimeInterval(isoStr);
            System.out.println(isoStr + " means " + t1.toString());
        }
    }

	@Test
	public void testInvalidCreation() {
		try {
			new ISO8601DateTimeInterval(2000, 1000);
			fail("Expected to throw IllegalArgumentException");
		} catch (IllegalArgumentException e) {
			// expected invalid time order
		}

		try {
			new ISO8601DateTimeInterval("2011-06-21T18:06:37.039Z--2011-01-01T18:06:37.039Z");
			fail("Expected to throw IllegalArgumentException");
		} catch (IllegalArgumentException e) {
			// expected invalid time order
		}
	}

	@Test
	public static void testIntervalCompareTo() {
        ISO8601DateTimeInterval t1 = new ISO8601DateTimeInterval();
		ISO8601DateTimePoint tp = new ISO8601DateTimePoint(t1.getStartTimeInMillis());
		assertEquals(0, t1.compareTo(tp));

		ISO8601DateTimeInterval t2 = new ISO8601DateTimeInterval(t1.getStartTimeInMillis(), t1.getEndTimeInMillis());

		assertEquals(0, t1.compareTo(t2));

		t2.setEndTimeInMillis(t2.getEndTimeInMillis() + 1000);
		assertEquals(-1, t1.compareTo(t2));

		t2.setStartTimeInMillis(t2.getStartTimeInMillis() - 5000);
		t2.setEndTimeInMillis(t2.getEndTimeInMillis() - 2000);
		assertEquals(1, t1.compareTo(t2));
	}

	@Test
    public static void testNullCompare() throws Exception {
		ISO8601DateTimePoint t1 = new ISO8601DateTimePoint();
		ISO8601DateTimePoint t2 = null;
		assertFalse(t1.equals(t2));
	}

	@Test
    public static void testNullTimeInterval() {
		ISO8601DateTimeInterval t1 = new ISO8601DateTimeInterval();
        ISO8601DateTimeInterval t2 = null;
		assertFalse(t1.equals(t2));
	}

    /**
     * Main method for running class tests.
     *
     * @param args standard command line arguments - ignored.
     */
    public static void main(String[] args) {
        JUnitCore.runClasses(TestISO8601DateTime.class);
    }
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestUTM.java<|end_filename|>
/****************************************************************************************
 *  TestUTM.java
 *
 *  Created: Mar 28, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy.test;

import junit.framework.TestSuite;
import junit.textui.TestRunner;

import java.util.Random;

import org.junit.Test;
import org.opensextant.geodesy.*;

import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertFalse;
import static org.junit.Assert.assertTrue;
import static org.junit.Assert.fail;

public class TestUTM {

    /**
     * This method tests UTM methods to determine Longitude Zones and their limits
     */
	@Test
	public void testLonZones() {
        int lonZone;
        char latBand = 'N';
        double lonDeg;
        // Test Lon Zones by longitude degrees (1° increment lons at lat band 'N')
        for (lonDeg = -180.0; lonDeg < +180.0; lonDeg += 1.0) {
            lonZone = UTM.getLonZone(lonDeg, latBand);
            assertTrue((UTM.minLonDegrees(lonZone, latBand) <= lonDeg) &&
                       (lonDeg < UTM.maxLonDegrees(lonZone, latBand)));
        }
        // Test each Lon Zones for Central Meridian inclusion
        for (lonZone = 1; lonZone <= 60; lonZone += 1) {
            lonDeg = UTM.getCentralMeridian(lonZone, latBand).inDegrees();
            assertTrue((UTM.minLonDegrees(lonZone, latBand) <= lonDeg) &&
                       (lonDeg < UTM.maxLonDegrees(lonZone, latBand)));
        }
        // Test random longitude degrees and random lon zones
        Random r = new Random();
        for (int i = 0; i < 1000; i++) {
            lonDeg = (r.nextDouble() * 360.0) - 180.0;  // lonDeg
            lonZone = UTM.getLonZone(lonDeg, latBand);
            assertTrue((UTM.minLonDegrees(lonZone, latBand) <= lonDeg) &&
                       (lonDeg < UTM.maxLonDegrees(lonZone, latBand)));
            lonZone = r.nextInt(60) + 1;
            lonDeg = UTM.getCentralMeridian(lonZone, latBand).inDegrees();
            assertTrue((UTM.minLonDegrees(lonZone, latBand) <= lonDeg) &&
                       (lonDeg < UTM.maxLonDegrees(lonZone, latBand)));
        }
    }

    /**
     * This method tests UTM methods to determine Latitude Bands and their limits
     */
	@Test
	public void testLatBands() {
        char latBand;
        double latDeg;
        // Test the latitude band char mappings and band ranges
        // (note that 84° is considered in band 'X', hence use of '<=' for max)
        for (latDeg = -90.0; latDeg < 90.0; latDeg += 1.0) {
            try {
                latBand = UTM.getLatBand(latDeg);
                assertTrue((UTM.minLatDegrees(latBand) <= latDeg) &&
                           (latDeg <= UTM.maxLatDegrees(latBand)));
            } catch (Exception ex) {
                //System.err.println(ex.toString());
            }
        }
        // Test random latitude degrees and lat band constraints
        Random r = new Random();
        for (int i = 0; i < 1000; i++) {
            latDeg = (164.0 * r.nextDouble()) - 80.0;
            latBand = UTM.getLatBand(latDeg);
            assertTrue((UTM.minLatDegrees(latBand) <= latDeg) &&
                       (latDeg < UTM.maxLatDegrees(latBand)));
        }
    }

	@Test
	public void testInvalidCreation() {
		try {
			// Hemisphere 'X', should be 'N' or 'S'
			new UTM(31, 'X', 353305.0, 7100467.0);
			fail("Expected to throw IllegalArgumentException");
		} catch (IllegalArgumentException e) {
			// expected
		}

		try {
			// outside UTM longitudinal zone valid range (1 to 60)
			new UTM(0, 'N', 353305.0, 7100467.0);
			fail("Expected to throw IllegalArgumentException");
		} catch (IllegalArgumentException e) {
			// expected
		}

		try {
			// outside easting valid range: 100,000 to 900,000 meters
			new UTM(31, 'N', 0.0, 7100467.0);
			fail("Expected to throw IllegalArgumentException");
		} catch (IllegalArgumentException e) {
			// expected
		}

		try {
			// outside northing valid range (0 to 10,000,000 meters)
			new UTM(31, 'N', 353305.0, -1.0);
			fail("Expected to throw IllegalArgumentException");
		} catch (IllegalArgumentException e) {
			// expected
		}
	}

	@Test
	public void testEquals() {
		Geodetic2DPoint g1 = new Geodetic2DPoint(
				new Longitude(-79, 23, 13.7),
				new Latitude(43, 38, 33.24));
		UTM u1 = new UTM(g1);

		assertTrue(u1.equals(u1));

		UTM u2 = new UTM(u1.getLonZone(), u1.getHemisphere(), u1.getEasting(), u1.getNorthing());
		assertEquals(u1, u2);
		assertEquals(u1.hashCode(), u2.hashCode());
		assertEquals(u1.getLatitude(), u2.getLatitude());
		assertEquals(u1.getLongitude(), u2.getLongitude());

		UTM u3 = new UTM(u1.getEllipsoid(), u1.getLonZone(), u1.getHemisphere(), u1.getEasting(), u1.getNorthing());
		assertEquals(u1, u3);
		assertEquals(u1.hashCode(), u3.hashCode());

		UTM u4 = new UTM(u1.getEllipsoid(), g1);
		assertEquals(u1, u4);
		assertEquals(u1.hashCode(), u4.hashCode());

		UTM u5 = new UTM(u1.getEllipsoid(), g1.getLongitude(), g1.getLatitude());
		assertEquals(u1, u5);
		assertEquals(u1.hashCode(), u5.hashCode());

		UTM u6 = new UTM(u1.getLonZone(), u1.getHemisphere(), u1.getEasting(), u1.getNorthing() + 500);
		assertFalse(u1.equals(u6));

		UTM u7 = null;
		assertFalse(u1.equals(u7));

		Object other = g1;
		assertFalse(u1.equals(other));
	}

    /**
     * This method tests both some known point cases as well as some randomly
     * selected Geodetic points, for consistency in round trip conversions at
     * 2 fractional digits of precision (hundreths of arc seconds and hundreths
     * of meters)
     */
	@Test
    public void testProjections() {
        // Test Case : Toronto's CNN Tower
        //   WGS 84 Geodetic2DPoint (lon-lat): (79° 23' 13.70" W, 43° 38' 33.24" N)
        //   Zone 17, Band T (hemisphere 'N'), 630084m east, 4833439m north
        Geodetic2DPoint g1, g2;
        UTM u1, u2;
        int fractDig = 2; // Fractional digits to print in toString conversions

        // utm.toGeodetic(17, 'N', 630084.30, 4833438.55) => (79° 23' 13.70" W, 43° 38' 33.24" N)
        //  MGRS? 17TPJ 30084 33439
        g1 = new Geodetic2DPoint(
                new Longitude(-79, 23, 13.700468),
                new Latitude(43, 38, 33.240046));
        u1 = new UTM(g1);
        u2 = new UTM(u1.getLonZone(), u1.getHemisphere(), u1.getEasting(), u1.getNorthing());
        g2 = u2.getGeodetic();
        assertEquals(g1.toString(fractDig), g2.toString(fractDig));

        u1 = new UTM(17, 'N', 630084.30, 4833438.55);
        g1 = u1.getGeodetic();
        u2 = new UTM(g1);
        assertEquals(u1.toString(fractDig), u2.toString(fractDig));

        // Another test case
        //   NAD 83 Geodetic2DPoint (lon-lat): (97° 44' 25.19" W, 30° 16' 28.82" N)
        //   UTM equivalent: Zone 14 R 621160.98, 3349893.53
        g1 = new Geodetic2DPoint(new Longitude(-97, 44, 25.19), new Latitude(30, 16, 28.82));
        u1 = new UTM(g1);
        u2 = new UTM(u1.getLonZone(), u1.getHemisphere(), u1.getEasting(), u1.getNorthing());
        g2 = u2.getGeodetic();
        assertEquals(g1.toString(fractDig), g2.toString(fractDig));

        u1 = new UTM(14, 'N', 621160.98, 3349893.53);
        g1 = u1.getGeodetic();
        u2 = new UTM(g1);
        assertEquals(u1.toString(fractDig), u2.toString(fractDig));

        // The following test case is for a round off error that causes lonZone shift
        // if not detected
        g1 = new Geodetic2DPoint(new Longitude(0, 0, 0.0), new Latitude(64, 0, 0.0));
        u1 = new UTM(g1);
        u2 = new UTM(u1.getLonZone(), u1.getHemisphere(), u1.getEasting(), u1.getNorthing());
        g2 = u2.getGeodetic();
        assertEquals(g1.toString(fractDig), g2.toString(fractDig));

        u1 = new UTM(31, 'N', 353305.0, 7100467.0);
        g1 = u1.getGeodetic();
        u2 = new UTM(g1);
        assertEquals(u1.toString(fractDig), u2.toString(fractDig));

        // Do some random Geodetic point conversion round trips
        Random r = new Random();
        for (int i = 0; i < 1000; i++) {
            double lonDeg = (r.nextDouble() * 360.0) - 180.0;
            double latDeg = (164.0 * r.nextDouble()) - 80.0;
            g1 = new Geodetic2DPoint(
                    new Longitude(lonDeg, Angle.DEGREES), new Latitude(latDeg, Angle.DEGREES));
            u1 = new UTM(g1);
            u2 = new UTM(u1.getLonZone(), u1.getHemisphere(), u1.getEasting(), u1.getNorthing());
            g2 = u2.getGeodetic();
            assertEquals(g1.toString(fractDig), g2.toString(fractDig));
        }
    }

	@Test
	public void testLonZone() {
		assertEquals(31, UTM.getLonZone(360, 'N'));
		assertEquals(31, UTM.getLonZone(-360, 'N'));
	}

	@Test(expected = IllegalArgumentException.class)
	public void testMinNorthing() {
		UTM.minNorthing('I');
		// UTM latitude band character ('C' to 'X", not including 'I' or 'O')
	}

	@Test(expected = IllegalArgumentException.class)
	public void testMaxNorthing() {
		UTM.maxNorthing(1, 'O');
		// UTM latitude band character ('C' to 'X", not including 'I' or 'O')
	}

	@Test(expected = IllegalArgumentException.class)
	public void testLatZoneLowerRange() {
		UTM.getLatBand(-81);
		// expected out of legal range (-80 deg to 84 deg) for UTM
	}

	@Test(expected = IllegalArgumentException.class)
	public void testLatZoneUpperRange() {
		UTM.getLatBand(90);
		// expected out of legal range (-80 deg to 84 deg) for UTM
	}

	@Test(expected = IllegalArgumentException.class)
	public void testValidateLonZoneLowerRange() {
		UTM.validateLonZone(0);
	}

	@Test(expected = IllegalArgumentException.class)
	public void testValidateLonZoneUpperRange() {
		UTM.validateLonZone(61);
	}

	@Test
	public void testToString() {
		UTM u1 = new UTM(14, 'N', 621160.08, 3349893.03);
		// tp.toString(0) WGS 84 UTM 14 N hemisphere 621160m E, 3349893m N
		// tp.toString(4) WGS 84 UTM 14 N hemisphere 621160.0800m E, 3349893.0300m N

        String base = u1.toString();
		assertEquals(base, u1.toString(0)); // toString() same as toString(0)
		int prevLen = base.length();
        String prefix = base.substring(0, base.indexOf("m E"));   // WGS 84 UTM 14 N hemisphere 621160
		// System.out.printf("prefix=[%s]%n", prefix);

		for(int i=1; i < 6; i++) {
            String s = u1.toString(i);
			// System.out.println(s);
            int len = s.length();
            assertTrue(len >= prevLen + 2);
            // simple tests: each should start with same easting whole number
			assertEquals(prefix, s.substring(0, prefix.length()));
            // assertTrue("prefix fails to match target", s.startsWith(prefix));
            prevLen = len;
		}
	}

    /**
     * Main method for running class tests.
     *
     * @param args standard command line arguments - ignored.
     */
    public static void main(String[] args) {
        TestSuite suite = new TestSuite(TestUTM.class);
        new TestRunner().doRun(suite);
    }
}


<|start_filename|>library/src/test/java/org/opensextant/geodesy/test/TestTransverseMercator.java<|end_filename|>
/****************************************************************************************
 *  TestTransverseMercator.java
 *
 *  Created: Mar 28, 2007
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantibility and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy.test;

import junit.framework.TestCase;
import junit.framework.TestSuite;
import junit.textui.TestRunner;

import org.opensextant.geodesy.*;

import java.util.Random;

public class TestTransverseMercator extends TestCase {

    /**
     * This method tests TransverseMercator projections along the Prime Meridian,
     * the Equator, and over randomly selected points within a typical UTM zone.
     */
    public void testProjections() {
        // Circumference of the earth at the Equator is 40,075.16 kilometers
        // so 1 arc minute along the Equator should equal 1855.3314815 meters
        // (note that 1852m in the official length for a standard nautical mile)
        // 6' == 0.1 degree is 11,131.989 meters
        Longitude lon;
        Latitude lat;
        Topocentric2DPoint en;
        TransverseMercator tm = new TransverseMercator(false);

        // Longitude tests vary central meridian around the globe at the Equator
        // and test for a consistent offset at 1° east of the central meridian.
        // Eastings should all be nearly 111325.18096 m, and Northings 0 m.
        lat = new Latitude(0.0, Angle.DEGREES);
        for (int cm = -181; cm <= +180; cm++) {
            tm.setCentralMeridian(new Longitude(cm, Angle.DEGREES));
            lon = new Longitude(cm + 1.0, Angle.DEGREES);
            en = tm.toTransverseMercator(lon, lat);
            assertEquals(en.getEasting(), 111325.18096, 0.00001);
            assertEquals(en.getNorthing(), 0.0);
        }

        // Latitude tests vary the Latitude from -90° to +88° by steps of 1°
        // Northings should differ from previous value by no less than 110574m and
        // by no more than 111693m. This is just a sanity range check, since precisely
        // correct values would require a table to compare against
        int prevNorthing = -10001405;
        lon = new Longitude(0.0, Angle.DEGREES);
        tm.setCentralMeridian(lon);
        for (int ol = -90; ol <= +88; ol++) {
            lat = new Latitude(ol + 1.0, Angle.DEGREES);
            en = tm.toTransverseMercator(lon, lat);
            int currNorthing = (int) en.getNorthing();
            assertEquals(en.getEasting(), 0.0);
            assertEquals(currNorthing - prevNorthing, 111134, 560);
            prevNorthing = currNorthing;
        }

        // Random point tests project and unproject to test for equality in transform
        // This is based on default Central Meridian (at the Prime Meridian), so
        // Longitude values are restricted to be within 5° of the Prime Meridian, and
        // Latitude values are restricted to be between -89° and +89°.  Accuracy should
        // be to hundredths of an arc second (2 fractional digits in toString format).
        Geodetic2DPoint p1, p2;
        Random r = new Random();
        int fractDig = 2;
        for (int i = 0; i < 1000; i++) {
            lon = new Longitude(r.nextDouble() * 5.0 *
                    ((r.nextBoolean()) ? -1 : +1), Angle.DEGREES);
            lat = new Latitude(r.nextDouble() * 89.0 *
                    ((r.nextBoolean()) ? -1 : +1), Angle.DEGREES);
            p1 = new Geodetic2DPoint(lon, lat);
            en = tm.toTransverseMercator(lon, lat);
            p2 = tm.toGeodetic(en.getEasting(), en.getNorthing());
			assertEquals("Geodetic2DPoint.toString()", p1.toString(fractDig), p2.toString(fractDig));
        }

    }

    /**
     * Main method for running class tests.
     *
     * @param args standard command line arguments - ignored.
     */
    public static void main(String[] args) {
        TestSuite suite = new TestSuite(TestTransverseMercator.class);
        new TestRunner().doRun(suite);
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/UTM.java<|end_filename|>
/************************************************************************************
 *  UTM.java
 *
 *  Created: Dec 21, 2006
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 * The program is provided "as is" without any warranty express or implied, including
 * the warranty of non-infringement and the implied warranties of merchantability and
 * fitness for a particular purpose.  The Copyright owner will not be liable for any
 * damages suffered by you as a result of using the Program. In no event will the
 * Copyright owner be liable for any special, indirect or consequential damages or
 * lost profits even if the Copyright owner has been advised of the possibility of
 * their occurrence.
 *
 ***********************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.util.Map;
import java.util.HashMap;
import java.io.Serializable;
import java.text.DecimalFormat;

/**
 * <P>
 * The UTM (Universal Transverse Mercator) class contains methods to parse and format
 * between UTM coordinate strings and their geodetic (longitude and latitude)
 * equivalents. A UTM object is defined only in terms of the Ellipsoid data model
 * against which projections are made. The default constructor uses the WGS 84
 * ellipsoid.
 * <P/><P>
 * Notes:
 * <P/><P>
 * Projection: Transverse Mercator (Gauss-Kruger type) in zones 6 deg wide.
 * <P/><P>
 * Longitude of Origin: Central meridian (CM) of each projection zone (degrees of 3,
 * 9, 15, 21, 27, 33, 39, 45, 51, 57, 63, 69, 75, 81, 87, 93, 99, 105, 111, 117, 123,
 * 129, 135, 141, 147, 153, 159, 165, 171, 177, E and W). <P>
 * <P/><P>
 * Latitude of Origin: 0 deg (the Equator).
 * <P/> <P>
 * Unit: Meter.
 * <P/><P>
 * False Northing: 0 meters at the Equator for the Northern Hemisphere;
 * 10,000,000 meters at the Equator for the Southern Hemisphere. 
 * <P/><P>
 * False Easting: 500,000 meters at the Central Meridian (CM) of each zone. 
 * <P/><P>
 * Scale Factor at the Central Meridian: 0.9996. 
 * <P/><P>
 * Latitude Limits of System: From 80 deg S to 84 deg N.
 * </P><P>
 * Overlap: On large-scale maps and trig lists, the data for each zone, datum, or
 * ellipsoid overlaps the adjacent zone, datum, or ellipsoid a minimum of 40
 * kilometers. The UTM grid extends to 80.5 deg S and 84.5 deg N, providing a
 * 30-minute overlap with the UPS grid.
 * <P/><P>
 * Limits of Projection Zones: The zones are bounded by meridians, the longitudes
 * of which are multiples of 6 deg east and west of the prime meridian.
 * </P><P>
 * Universal Transverse Mercator (UTM) coordinates define two dimensional,
 * horizontal, positions. The sixty UTM zone numbers designate 6 degree wide
 * longitudinal strips extending from 80 degrees South latitude to 84 degrees North
 * latitude. UTM zone characters are letters which designate 8 degree zones
 * extending north and south from the equator. Beginning at 80 deg south and
 * proceeding northward, twenty bands are lettered C through X, omitting I and O.
 * These bands are all 8 deg wide except for band X which is 12 deg wide (between
 * 72-84 deg N).
 * </P><P>
 * There are special UTM zones between 0 degrees and 36 degrees longitude above 72
 * degrees latitude and a special zone 32 between 56 degrees and 64 degrees north
 * latitude: </P>
 * <P>
 * UTM Zone 32 has been widened to 9 deg (at the expense of zone 31) between
 * latitudes 56 deg and 64 deg (band V) to accommodate southwest Norway. Thus zone
 * 32 it extends westwards to 3 deg E in the North Sea.
 * </P><P>
 * Similarly, between 72 deg and 84 deg (band X), zones 33 and 35 have been widened
 * to 12 deg to accommodate Svalbard. To compensate for these 12 deg wide zones,
 * zones 31 and 37 are widened to 9 deg and zones 32, 34, and 36 are eliminated.
 * Thus the W and E boundaries of zones are 31: 0 - 9 deg E, 33: 9 - 21 deg E,
 * 35: 21 - 33 deg E and 37: 33 - 42 deg E.
 * </P>
 */
public class UTM implements Serializable {
    private static final long serialVersionUID = 1L;

    // Static Class Constants
    public static final double MAX_NORTH_LATDEG = +84.0;
    public static final double MIN_SOUTH_LATDEG = -80.0;
    private static final double MIN_EASTING = 100000.0;
    private static final double MAX_EASTING = 900000.0;
    private static final double MIN_NORTHING = 0.0;
    private static final double MAX_NORTHING = 10000000.0;
    private static final double SCALE = 0.9996;
    private static final Latitude EQUATOR = new Latitude(0.0);

    public static final double FALSE_EASTING = 500000.0;
    public static final double FALSE_NORTHING = 10000000.0;

    private static final double ROUNDING_POS = 1e+6;
    private static final double ROUNDING_NEG = 1e-6;

    // Cached Central Meridian values to avoid object creation and computation at run-time
    private static final Longitude[] CENTRAL_MERIDIAN = new Longitude[61];

    static {
        for (int lonZone = 1; lonZone <= 60; lonZone++) {
            CENTRAL_MERIDIAN[lonZone] =
                    new Longitude((6.0 * lonZone + ((lonZone >= 31) ? -183.0 : 177.0)),
                            Angle.DEGREES);
        }
    }

    // Initialize min northings for lat bands
    // These tables are used to help assign and validate MGRS square identifiers for northings
    private static final Map<Character, Integer> MIN_NORTHINGS =
            new HashMap<Character, Integer>();
    private static final Map<Character, Integer> MAX_NORTHINGS =
            new HashMap<Character, Integer>();

    static {
        Longitude lon;
        Longitude centLon = new Longitude(-3.0, Angle.DEGREES);
        Longitude edgeLon = new Longitude(-6.0, Angle.DEGREES);

        for (char latBand = 'C'; latBand <= 'X'; latBand++) {
            if ((latBand == 'I') || (latBand == 'O')) continue;

            // min Northings are inclusive (northing must be greater than or equal to this min)
            // min Northings for Northern Hemisphere are along a central meridian, but for
            // the Southern Hemisphere, they occur at the edges of a lonZone (+/- 3 deg from CM)
            lon = (latBand < 'N') ? edgeLon : centLon;
            Latitude minLat = new Latitude(UTM.minLatDegrees(latBand), Angle.DEGREES);
            double minNorthing = new UTM(new Geodetic2DPoint(lon, minLat)).getNorthing();
            MIN_NORTHINGS.put(latBand, (int) Math.round(minNorthing));

            // max Northings are exclusive (northing must be less than this max)
            // max Northings for Southern Hemisphere are along a central meridian, but for
            // the Northern Hemisphere, they occur at the edges of a lonZone (typically +/- 3 deg
            // from CM). Special case at equator for southern hemisphere latBand 'M' is
            // FALSE_NORTHING
            double maxNorthing;
            if (latBand == 'M') {
                maxNorthing = UTM.FALSE_NORTHING;
            } else {
                lon = (latBand < 'N') ? centLon : edgeLon;
                Latitude maxLat = new Latitude(UTM.maxLatDegrees(latBand), Angle.DEGREES);
                maxNorthing = new UTM(new Geodetic2DPoint(lon, maxLat)).getNorthing();
            }
            MAX_NORTHINGS.put(latBand, (int) Math.round(maxNorthing));
        }
    }

    private static final Longitude V31_CENTRAL_MERIDIAN = new Longitude(1.5, Angle.DEGREES);
    private static final Longitude X31_CENTRAL_MERIDIAN = new Longitude(4.5, Angle.DEGREES);
    private static final Longitude V32_CENTRAL_MERIDIAN = new Longitude(7.5, Angle.DEGREES);
    private static final Longitude X37_CENTRAL_MERIDIAN = new Longitude(37.5, Angle.DEGREES);

    private static final int V31_MAX_NORTHING = (int) Math.round(
            new UTM(new Longitude(0.0, Angle.DEGREES), new Latitude(63, 59, 59.99)).getNorthing());
    private static final int V32_MAX_NORTHING = (int) Math.round(
            new UTM(new Longitude(3.0, Angle.DEGREES), new Latitude(63, 59, 59.99)).getNorthing());
    private static final int X31_X37_MAX_NORTHING = (int) Math.round(
            new UTM(new Longitude(0.0, Angle.DEGREES), new Latitude(83, 59, 59.99)).getNorthing());
    private static final int X33_X35_MAX_NORTHING = (int) Math.round(
            new UTM(new Longitude(21.0, Angle.DEGREES), new Latitude(83, 59, 59.99)).getNorthing());

    // Static Class Methods

    /**
     * This method determines the UTM longitudinal zone number for a given lon (in decimal degrees).
     *
     * @param lonDeg  degrees of longitude to map to UTM longitudinal zone number
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return int for longitudinal zone number (1 to 60)
     */
    public static int getLonZone(double lonDeg, char latBand) {
        // Round value to correct for accumulation of numerical error in UTM projection
        lonDeg = Math.floor(Math.round(lonDeg * ROUNDING_POS) * ROUNDING_NEG);

        // Normalize the longitude in degrees if necessary
        while (lonDeg < -180.0) lonDeg += 360.0;
        while (lonDeg >= +180.0) lonDeg -= 360.0;
        // Compute nominal zone value for most cases
        int lonZone = 1 + (int) ((lonDeg + 180.0) / 6.0);
        // Five zones have special boundaries to be checked, override nominal if necessary
        if (latBand == 'V') {
            if ((3.0 <= lonDeg) && (lonDeg < 12.0)) lonZone = 32;
        } else if (latBand == 'X') {
            if ((0.0 <= lonDeg) && (lonDeg < 9.0)) lonZone = 31;
            else if ((9.0 <= lonDeg) && (lonDeg < 21.0)) lonZone = 33;
            else if ((21.0 <= lonDeg) && (lonDeg < 33.0)) lonZone = 35;
            else if ((33.0 <= lonDeg) && (lonDeg < 42.0)) lonZone = 37;
        }
        return lonZone;
    }

    /**
     * This method determines the minimum longitude (in degrees) for a given lon zone
     * and lat band.
     *
     * @param lonZone UTM Longitudinal Zone (1 to 60)
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return minimum longitude (in degrees) for this UTM cell (lon zone and lat band)
     */
    public static double minLonDegrees(int lonZone, char latBand) {
        double lonDeg = (lonZone * 6.0) - 186.0;
        // Four zones have min lon overrides
        if ((latBand == 'V') && (lonZone == 32)) lonDeg = 3.0;
        else if (latBand == 'X') {
            if (lonZone == 33) lonDeg = 9.0;
            else if (lonZone == 35) lonDeg = 21.0;
            else if (lonZone == 37) lonDeg = 33.0;
        }
        return lonDeg;
    }

    /**
     * This accessor method gets the appropriate central meridian for the given UTM lon zone
     * and lat band.
     *
     * @param lonZone UTM Longitudinal Zone (1 to 60)
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return Longitude of appropriate central meridian for the given UTM cell (zone and band)
     */
    public static Longitude getCentralMeridian(int lonZone, char latBand) {
        Longitude cm = CENTRAL_MERIDIAN[lonZone];
        if (latBand == 'V') {
            if (lonZone == 31) cm = V31_CENTRAL_MERIDIAN;
            else if (lonZone == 32) cm = V32_CENTRAL_MERIDIAN;
        } else if (latBand == 'X') {
            if (lonZone == 31) cm = X31_CENTRAL_MERIDIAN;
            else if (lonZone == 37) cm = X37_CENTRAL_MERIDIAN;
        }
        return cm;
    }

    /**
     * This method determines the maximum longitude (in degrees) for a given lon zone & lat band.
     *
     * @param lonZone UTM Longitudinal Zone (1 to 60)
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return maximum longitude (in degrees) for this UTM cell (lon zone and lat band)
     */
    public static double maxLonDegrees(int lonZone, char latBand) {
        double lonDeg = (lonZone * 6.0) - 180.0;
        // Four zones have max lon overrides
        if ((latBand == 'V') && (lonZone == 31)) lonDeg = 3.0;
        else if (latBand == 'X') {
            if (lonZone == 31) lonDeg = 9.0;
            else if (lonZone == 33) lonDeg = 21.0;
            else if (lonZone == 35) lonDeg = 33.0;
        }
        return lonDeg;
    }

    /**
     * This method determines the UTM latitudinal band char for a given lat (in decimal degrees).
     *
     * @param latDeg degrees of latitude to map to UTM latitude band character
     * @return character representing latitude band ('C' to 'X', skipping 'I' and 'O')
     * @throws IllegalArgumentException error if latitude value is out of valid range
     */
    public static char getLatBand(double latDeg) {
        // validate that latitude is within proper range (allow half degree overlap with UPS)
        if ((latDeg < -80.5) || (84.5 < latDeg))
            throw new IllegalArgumentException("Latitude value '" + latDeg +
                    "' is out of legal range (-80 deg to 84 deg) for UTM");

        // Round value to correct for accumulation of numerical error in UTM projection
        latDeg = Math.floor(Math.round(latDeg * ROUNDING_POS) * ROUNDING_NEG);

        // Restrict calculated index to 20 8 deg-wide bands between -80 deg and +80 deg
        int i;
        if (latDeg < -80.0) i = 0;                  // correct for extra 'C' band range
        else if (80.0 <= latDeg) i = 19;            // correct for extra 'X' band range
        else i = (int) (10.0 + (latDeg / 8.0));     // index values for all others

        char latBand = (char) ('C' + i);
        if (i > 10) latBand += 2;
        else if (i > 5) latBand += 1;               // adjust for missing 'I' and 'O'
        return latBand;
    }

    /**
     * This method returns the inclusive minimum latitude (in decimal degrees)
     * for the given UTM latitude band char (extra half degree of allowance not included).
     *
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return minimum latitude for the given band, in decimal degrees
     * @throws IllegalArgumentException - error if lat band char is out of range
     */
    public static double minLatDegrees(char latBand) {
        validateLatBand(latBand);

        int i = (int) latBand - (int) 'C';
        if (latBand > 'N') i -= 2;
        else if (latBand > 'H') i -= 1;             // adjust for missing 'I' and 'O'
        return (8.0 * (i - 10.0));
    }

    /**
     * This method returns the exclusive maximum latitude (in decimal degrees)
     * for the given UTM latitude band char (extra half degree of allowance not included).
     *
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return maximum latitude for the given band, in decimal degrees
     * @throws IllegalArgumentException - error if lat band char is out of range
     */
    public static double maxLatDegrees(char latBand) {
        validateLatBand(latBand);

        int i = (int) latBand - (int) 'C';
        if (latBand > 'N') i -= 2;
        else if (latBand > 'H') i -= 1;             // adjust for missing 'I' and 'O'
        double w = (latBand == 'X') ? 12.0 : 8.0;   // Band 'X' has extra 4.0 deg N
        return ((8.0 * (i - 10.0)) + w);
    }

    /**
     * This method returns the minimum northing value (in meters) for the specified
     * latitude band.
     *
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return minimum northing value for the specified latitude band
     * @throws IllegalArgumentException - error if the latitude band is invalid
     */
    public static int minNorthing(char latBand) {
        Integer minNorthing = UTM.MIN_NORTHINGS.get(latBand);
        if (minNorthing == null)
            throw new IllegalArgumentException("Invalid latitude band '" + latBand + "'");
        return minNorthing;
    }

    /**
     * This method returns the maximum northing value (in meters) for the specified
     * latitude band.
     *
     * @param lonZone UTM Longitudinal Zone (1 to 60)
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return maximum northing value for the specified latitude band
     * @throws IllegalArgumentException error if the latitude band is invalid
     */
    public static int maxNorthing(int lonZone, char latBand) {
        Integer maxNorthing = UTM.MAX_NORTHINGS.get(latBand);
        if (maxNorthing == null)
            throw new IllegalArgumentException("Invalid latitude band '" + latBand + "'");
        // See if we need to override for special exception cells
        if (latBand == 'V') {
            if (lonZone == 31) maxNorthing = V31_MAX_NORTHING;
            else if (lonZone == 32) maxNorthing = V32_MAX_NORTHING;
        } else if (latBand == 'X') {
            if ((lonZone == 31) || (lonZone == 37)) maxNorthing = X31_X37_MAX_NORTHING;
            else if ((lonZone == 33) || (lonZone == 35)) maxNorthing = X33_X35_MAX_NORTHING;
        }
        return maxNorthing;
    }

    /**
     * This method determines the hemisphere ('N' or 'S') for the given lat band.
     *
     * @param latBand UTM latitude band character ('C' to 'X", not including 'I' or 'O')
     * @return hemisphere character ('N' for Northern, or 'S' for Southern)
     * @throws IllegalArgumentException error if lat band char is out of range
     */
    public static char getHemisphere(char latBand) {
        validateLatBand(latBand);
        return (latBand < 'N') ? 'S' : 'N';
    }

    /**
     * This method tests longitudinal zone to see if it is valid (by itself).
     *
     * @param lonZone UTM Longitudinal Zone (1 to 60)
     * @throws IllegalArgumentException error if UTM Longitudinal Zone is invalid
     */
    public static void validateLonZone(int lonZone) {
        if ((lonZone < 1) || (60 < lonZone))
            throw new IllegalArgumentException("UTM longitudinal zone '" + lonZone +
                    "' is outside of valid range (1 to 60)");
    }

    /**
     * This method tests latitudinal zone to see if it is valid (by itself).
     *
     * @param latBand UTM Latitudinal Band ('C' to 'X', but not 'I' or 'O')
     * @throws IllegalArgumentException error if UTM Latitudinal Band is invalid
     */
    public static void validateLatBand(char latBand) {
        if ((latBand < 'C') || ('X' < latBand) || (latBand == 'I') || (latBand == 'O'))
            throw new IllegalArgumentException("UTM latitudinal band '" + latBand +
                    "' is not valid");
    }

    /**
     * This method tests longitudinal zone and latitudinal band to make sure they are
     * consistent together.
     *
     * @param lonZone UTM Longitudinal Zone (1 to 60)
     * @param latBand UTM Latitudinal Band ('C' to 'X', but not 'I' or 'O')
     * @throws IllegalArgumentException error if Zone and Band combination is invalid
     */
    public static void validateZoneAndBand(int lonZone, char latBand)
            throws IllegalArgumentException {
        validateLonZone(lonZone);
        validateLatBand(latBand);
        if ((latBand == 'X') && ((lonZone == 32) || (lonZone == 34) || (lonZone == 36)))
            throw new IllegalArgumentException("Invalid longitude zone '" +
                    lonZone + "' in latitude band 'X'");
    }

    /*
     * This method tests hemisphere character to see if it is valid (by itself).
     *
     * @param hemisphere UTM Hemisphere character ('N' for North, 'S' for South)
     * @throws IllegalArgumentException error if hemisphere character is invalid
     */
    private static void validateHemisphere(char hemisphere) {
        if ((hemisphere != 'N') && (hemisphere != 'S'))
            throw new IllegalArgumentException("Invalid hemisphere '" +
                    hemisphere + "', should be 'N' or 'S'");
    }

    /*
     * This method tests easting value to see if it is within allowed positive numerical range.
     *
     * @param easting UTM Easting value in meters
     * @throws IllegalArgumentException error if easting value is out of range
     */
    private static void validateEasting(double easting) {
        if ((easting < MIN_EASTING) || (MAX_EASTING < easting))
            throw new IllegalArgumentException("Easting value '" + easting +
                    "' is outside of valid range (100,000 to 900,000 meters)");
    }

    /*
     * This method tests northing value to see if it is within allowed positive numerical range.
     *
     * @param northing UTM Northing value in meters
     * @throws IllegalArgumentException error if northing value is out of range
     */
    private static void validateNorthing(double northing) {
        if ((northing < MIN_NORTHING) || (MAX_NORTHING < northing))
            throw new IllegalArgumentException("Northing value '" + northing +
                    "' is outside of valid range (0 to 10,000,000 meters)");
    }

    // *************************** End of Static Definitions ********************************

    // Instance Variables
    @NonNull
    private final TransverseMercator tm;  // Transverse Mercator projection object (keeper of ellipsoid)
    private int lonZone;            // UTM Longitudinal Zone number (1 to 60)
    private char latBand;           // UTM Lat Band char('C' to 'X', not including 'I' or 'O')
    private char hemisphere;        // Hemisphere char ('N' for Northern, 'S' for Southern)
    private double easting;         // meters E of false easting origin (relative to lonZone's CM)
    private double northing;        // meters N of false northing origin (equator for 'N' hemisphere)

    @NonNull
    private Geodetic2DPoint lonLat; // Longitude and Latitude coordinates for this UTM position

    /**
     * This constructor takes an ellipsoid.  If WGS 84 is desired, instead use one of
     * the constructors without the ellipsoid parameter, since WGS 84 is the default.
     * The specified Geodetic2DPoint (lon-lat) point is converted to its UTM equivalent.
     *
     * @param ellip  Ellipsoid data model for earth
     * @param lonLat Geodetic2DPoint coordinates (lon-lat) to be converted to UTM
     * @throws NullPointerException if ellip or lonLat are null
     */
    public UTM(Ellipsoid ellip, Geodetic2DPoint lonLat) {
        tm = new TransverseMercator(true);    // Distortion warning will cause exception
        tm.setOriginLatitude(EQUATOR);
        tm.setScaleFactor(SCALE);
        tm.setEllipsoid(ellip);
        this.lonLat = lonLat;
        fromGeodetic();
    }

    /**
     * This constructor takes an ellipsoid.  If WGS 84 is desired, instead use one of
     * the constructors without the ellipsoid parameter, since WGS 84 is the default.
     * The specified longitude and latitude values are converted to a UTM equivalent.
     *
     * @param ellip Ellipsoid data model for earth
     * @param lon   Longitude of point to convert to UTM
     * @param lat   Latitude of point to convert to UTM
     * @throws NullPointerException if ellip is null
     */
    public UTM(Ellipsoid ellip, Longitude lon, Latitude lat) {
        tm = new TransverseMercator(true);    // Distortion warning will cause exception
        tm.setOriginLatitude(EQUATOR);
        tm.setScaleFactor(SCALE);
        tm.setEllipsoid(ellip);
        lonLat = new Geodetic2DPoint(lon, lat);
        fromGeodetic();
    }

    /**
     * This constructor takes an ellipsoid.  If WGS 84 is desired, instead use one of
     * the constructors without the ellipsoid parameter, since WGS 84 is the default.
     * Convert UTM parameters to equivalent geodetic (lon-lat) position, compute lat band.
     *
     * @param ellip      Ellipsoid data model for earth
     * @param lonZone    UTM longitudinal zone (1 to 60)
     * @param hemisphere character 'N' for Northern or 'S' for Southern hemisphere
     * @param easting    positive meters east of false adjusted central meridian for lonZone
     * @param northing   positive meters north of false adjusted origin (equator for 'N' hemisphere)
     * @throws NullPointerException     if ellip is null
     * @throws IllegalArgumentException input parameter error(s)
     */
    public UTM(Ellipsoid ellip, int lonZone, char hemisphere, double easting, double northing) {
        tm = new TransverseMercator(true);    // Distortion warning will cause exception
        tm.setOriginLatitude(EQUATOR);
        tm.setScaleFactor(SCALE);
        tm.setEllipsoid(ellip);
        this.lonZone = lonZone;
        this.hemisphere = hemisphere;
        this.easting = easting;
        this.northing = northing;
        toGeodetic(true);                     // Validate easting and northing values
    }

    /**
     * This constructor assumes the WGS 84 Ellipsoid (default for TransverseMercator).
     * The specified Geodetic2DPoint (lon-lat) point is converted to its UTM equivalent.
     *
     * @param lonLat Geodetic2DPoint coordinates (lon-lat) to be converted to UTM
     * @throws NullPointerException if lonLat is null
     */
    public UTM(Geodetic2DPoint lonLat) {
        tm = new TransverseMercator(true);
        tm.setOriginLatitude(EQUATOR);
        tm.setScaleFactor(SCALE);
        this.lonLat = lonLat;
        fromGeodetic();
    }

    /**
     * This constructor assumes the WGS 84 Ellipsoid (default for TransverseMercator).
     * The specified longitude and latitude values are converted to a UTM equivalent.
     *
     * @param lon Longitude of point to convert to UTM
     * @param lat Latitude of point to convert to UTM
     * @throws NullPointerException if lat or lon are null
     */
    public UTM(Longitude lon, Latitude lat) {
        tm = new TransverseMercator(true);    // Distortion warning will cause exception
        tm.setOriginLatitude(EQUATOR);
        tm.setScaleFactor(SCALE);
        lonLat = new Geodetic2DPoint(lon, lat);
        fromGeodetic();
    }

    /**
     * This constructor assumes the WGS 84 Ellipsoid (default for TransverseMercator).
     * Convert UTM parameters to equivalent geodetic (lon-lat) position, compute lat band,
     * and allows projection distortion warnings to optionally cause an Exception.
     *
     * @param lonZone             UTM longitudinal zone (1 to 60)
     * @param hemisphere          character 'N' for Northern or 'S' for Southern hemisphere
     * @param easting             positive meters east of false adjusted central meridian for lonZone
     * @param northing            positive meters north of false adjusted origin (equator for 'N' hemisphere)
     * @param distortionException boolean indicating whether distortion exceptions should be thrown
     * @throws IllegalArgumentException input parameter error(s)
     */
    public UTM(int lonZone, char hemisphere, double easting, double northing,
               boolean distortionException) {
        tm = new TransverseMercator(distortionException);    // Distortion warning exception?
        tm.setOriginLatitude(EQUATOR);
        tm.setScaleFactor(SCALE);
        this.lonZone = lonZone;
        this.hemisphere = hemisphere;
        this.easting = easting;
        this.northing = northing;
        toGeodetic(distortionException);                     // if distortion OK, skip validation too
    }

    /**
     * This constructor assumes the WGS 84 Ellipsoid (default for TransverseMercator).
     * Convert UTM parameters to equivalent geodetic (lon-lat) position, compute lat band.
     *
     * @param lonZone    UTM longitudinal zone (1 to 60)
     * @param hemisphere character 'N' for Northern or 'S' for Southern hemisphere
     * @param easting    positive meters east of false adjusted central meridian for lonZone
     * @param northing   positive meters north of false adjusted origin (equator for 'N' hemisphere)
     * @throws IllegalArgumentException input parameter error(s)
     */
    public UTM(int lonZone, char hemisphere, double easting, double northing) {
        this(lonZone, hemisphere, easting, northing, true);
    }

    /*
     * This method converts UTM projection (zone, hemisphere, easting and northing)
     * coordinates to geodetic (longitude and latitude) coordinates, based on the
     * current ellipsoid model.
     *
     * @throws IllegalArgumentException input parameter error(s)
     */
    private void toGeodetic(boolean validate) {
        // Validate input parameters
        validateLonZone(lonZone);
        validateHemisphere(hemisphere);
        if (validate) {
            validateEasting(easting);
            validateNorthing(northing);
        }

        // set nominal central meridian & adjust false values to regain signed offsets
        Longitude cm = CENTRAL_MERIDIAN[lonZone];
        double e = easting - FALSE_EASTING;
        double n = (hemisphere == 'S') ? northing - FALSE_NORTHING : northing;

        // Un-project to geodetic coordinates, assume no special zone override necessary
        tm.setCentralMeridian(cm);
        lonLat = tm.toGeodetic(e, n);

        // Determine lat band and validate cell combo
        // (lon zones 32, 34, & 36 are illegal in band 'X')
        latBand = getLatBand(lonLat.getLatitude().inDegrees());
//        if (validate) 
        validateZoneAndBand(lonZone, latBand);

        // Determine if special zone override makes it necessary to un-project again
        // with adjustments
        if (hemisphere == 'N') {
            // Use the un-projected northing latitude band to determine special case
            // lon zone CM values. Un-project again if central meridian has changed due
            // to special zones
            Longitude ocm = getCentralMeridian(lonZone, latBand);
            if (ocm != cm) {
                tm.setCentralMeridian(ocm);
                lonLat = tm.toGeodetic(e, n);
            }
        }
    }

    /*
     * This method converts from Geodetic2DPoint coordinates (lon-lat) to UTM parameters.
     */
    private void fromGeodetic() {
        Longitude lon = lonLat.getLongitude();
        Latitude lat = lonLat.getLatitude();

        latBand = UTM.getLatBand(lat.inDegrees());
        lonZone = UTM.getLonZone(lon.inDegrees(), latBand);

        hemisphere = getHemisphere(latBand);

        tm.setCentralMeridian(getCentralMeridian(lonZone, latBand));
        Topocentric2DPoint en = tm.toTransverseMercator(lon, lat);

        easting = en.getEasting() + FALSE_EASTING;
        double n = en.getNorthing();
        // Correct northing for latitudes barely south of the equator
        if ((n < 0.0) && (hemisphere == 'N')) n = 0.0;
        northing = (n < 0.0) ? n + FALSE_NORTHING : n;
    }

    /**
     * This method is used to get the currently set Ellipsoid earth model for making
     * projections during conversion to and from geodetic coordinates.
     *
     * @return the currently set Ellipsoid earth model being used by this UTM object
     */
    public Ellipsoid getEllipsoid() {
        return tm.getEllipsoid();
    }

    /**
     * This method returns the geodetic (lon-lat) point for this UTM position.
     *
     * @return Geodetic2DPoint object (lon-lat) coordinates
     */
    @NonNull
    public Geodetic2DPoint getGeodetic() {
        return lonLat;
    }

    /**
     * This method returns the Longitude part of the geodetic coordinates for this UTM object.
     *
     * @return Longitude coordinate of the geodetic (lon-lat) point for this UTM position
     */
    @NonNull
    public Longitude getLongitude() {
        return lonLat.getLongitude();
    }

    /**
     * This method returns the Latitude part of the geodetic coordinates for this UTM object.
     *
     * @return Latitude coordinate of the geodetic (lon-lat) point for this UTM position
     */
    @NonNull
    public Latitude getLatitude() {
        return lonLat.getLatitude();
    }

    /**
     * This method returns the UTM longitudinal zone (1 to 60) for this UTM object.
     *
     * @return int longitude zone number (1 to 60)
     */
    public int getLonZone() {
        return lonZone;
    }

    /**
     * This method returns the UTM latitudinal band character ('C' to 'X', not
     * including 'I' or 'O') for this UTM object.
     *
     * @return char for latitudinal band ('C' to 'X', not including 'I' or 'O')
     */
    public char getLatBand() {
        return latBand;
    }

    /**
     * This method returns the UTM hemisphere character ('N' for Northern,
     * 'S' for Southern) for this UTM object.
     *
     * @return char for hemisphere ('N' for Northern, 'S' for Southern)
     */
    public char getHemisphere() {
        return hemisphere;
    }

    /**
     * This method returns the UTM easting value (number of meters east of the central
     * meridian for the longitudinal zone, false adjusted to make this number always
     * positive).
     *
     * @return double easting value for this UTM object in positive meters
     */
    public double getEasting() {
        return easting;
    }

    /**
     * This method returns the UTM northing value (number of meters north of the equator or
     * falsely adjusted equator, to make this number always positive)
     *
     * @return double northing value for this UTM object in positive meters
     */
    public double getNorthing() {
        return northing;
    }

    @Override
    public boolean equals(Object o) {
        if (this == o) return true;
        if (o == null || getClass() != o.getClass()) return false;

        UTM utm = (UTM) o;

        if (Double.compare(utm.easting, easting) != 0) return false;
        if (hemisphere != utm.hemisphere) return false;
        if (latBand != utm.latBand) return false;
        if (lonZone != utm.lonZone) return false;
        if (Double.compare(utm.northing, northing) != 0) return false;
        if (!lonLat.equals(utm.lonLat)) return false;
        // TODO: following fails since TransverseMercator does not implement Object.equals(Object) so skip it
        // if (tm != null ? !tm.equals(utm.tm) : utm.tm != null) return false;
        return true;
    }

    @Override
    public int hashCode() {
        // int result;
        long temp;
        // TODO: following fails since TransverseMercator does not implement Object.hashCode() so skip it
        //result = tm != null ? tm.hashCode() : 0;
        // result = 31 * result + lonZone;
        int result = lonZone;
        result = 31 * result + (int) latBand;
        result = 31 * result + (int) hemisphere;
        temp = easting != +0.0d ? Double.doubleToLongBits(easting) : 0L;
        result = 31 * result + (int) (temp ^ (temp >>> 32));
        temp = northing != +0.0d ? Double.doubleToLongBits(northing) : 0L;
        result = 31 * result + (int) (temp ^ (temp >>> 32));
        result = 31 * result + lonLat.hashCode();
        return result;
    }

    /**
     * This method formats the UTM object as a printable string, using the specified number of
     * fractional digits for the easting and northing offsets.
     *
     * @param fractDig number of digits past the decimal point for easting and northing
     * @return String formatted UTM coordinate
     */
    public String toString(int fractDig) {
        StringBuilder fract = new StringBuilder("0");
        if (fractDig > 0) {
            fract.append('.');
            while (fractDig-- > 0) fract.append('0');
        }
        DecimalFormat offFmt = new DecimalFormat(fract.toString());
        return (getEllipsoid().getName() + " UTM " + lonZone + " " + hemisphere +
                " hemisphere " + offFmt.format(easting) + "m E, " +
                offFmt.format(northing) + "m N");
    }

    /**
     * This method formats the UTM object to a printable string, with easting and northing offsets
     * displayed to the nearest integer number of meters.
     *
     * @return String formatted UTM coordinate
     */
    public String toString() {
        return toString(0);
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Geodetic2DCircle.java<|end_filename|>
/******************************************************************************
 * Geodetic2DCircle.java Nov 3, 2009 8:20:19 PM psilvey
 *
 * (C) Copyright MITRE Corporation 2009
 *
 * The program is provided "as is" without any warranty express or implied,
 * including the warranty of non-infringement and the implied warranties of
 * merchantability and fitness for a particular purpose.  The Copyright
 * owner will not be liable for any damages suffered by you as a result of
 * using the Program.  In no event will the Copyright owner be liable for
 * any special, indirect or consequential damages or lost profits even if
 * the Copyright owner has been advised of the possibility of their
 * occurrence.
 *
 *****************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.util.ArrayList;

/**
 * The Geodetic2DCircle class represents an circle on the surface of the earth
 * (itself modeled as an Ellipsoid). Similar to the Geodetic2DPoint class,
 * instances of this object are lightweight in the sense that they represent
 * circles in terms of simple parameters: center and radius. The center point
 * is specified as a Geodetic2DPoint, and the radius in meters. Just as with
 * the Geodetic2DPoint class, the reference Ellipsoid earth model is not
 * explicitly carried with the objects, but may be needed for some uses of the
 * circle.
 */
public class Geodetic2DCircle {

    @NonNull private Geodetic2DPoint center;

    private double radius;

    /**
     * Default constructor makes a geodetic circle at the central meridian on
     * the equator (0, 0), with a radius of 0.
     */
    public Geodetic2DCircle() {
        this.center = new Geodetic2DPoint();
        this.radius = 0.0;
    }

    /**
     * This constructor takes center Geodetic2DPoint and the radius in meters.
     *
     * @param center Geodetic2DPoint of the center of this circle
     * @param radius length in meters of the radius of this circle
	 * @throws NullPointerException if center is null
     */
    public Geodetic2DCircle(Geodetic2DPoint center, double radius) {
		setCenter(center);
        this.radius = radius;
    }

    /**
     * Getter method for accessing the center of this circle object.
     *
     * @return Geodetic2DPoint at the center of this circle
     */
	@NonNull
    public Geodetic2DPoint getCenter() {
        return center;
    }

    /**
     * Setter method to assign the center Geodetic2DPoint of this circle object.
     *
     * @param center Geodetic2DPoint to make the center of this circle
	 * @throws NullPointerException if center is null
     */
    public void setCenter(Geodetic2DPoint center) {
		if (center == null) throw new NullPointerException();
        this.center = center;
    }

    /**
     * Getter method for accessing the radius of this circle object.
     *
     * @return double number of meters of the radius of this circle
     */
    public double getRadius() {
        return radius;
    }

    /**
     * Setter method to assign the radius of this circle object.
     *
     * @param radius double number of meters of this circle
     */
    public void setRadius(double radius) {
        this.radius = radius;
    }

    /**
     * Convenience method to iterate over boundary points of circle at nPoints
     * resolution.
     *
     * @param nPoints int number of points on boundary to use (1st is due South),
	 * 					should be greater than 0
     * @return Iterable collection of nPoints Geodetic2DPoints on circle boundary
	 * @throws ArithmeticException if nPoints = 0
     */
	@NonNull
    public Iterable<Geodetic2DPoint> boundary(int nPoints) {
        ArrayList<Geodetic2DPoint> ptList = new ArrayList<Geodetic2DPoint>(nPoints);
        Angle compassDirection = new Angle(-Math.PI); // Start from South direction
        Angle inc = new Angle(2.0 * Math.PI / (double) nPoints); // degree increments
        Geodetic2DArc arc = new Geodetic2DArc(center, radius, compassDirection);
        for (int i = 0; i < nPoints; i++) {
            ptList.add(arc.getPoint2());
            compassDirection = compassDirection.add(inc);
            arc.setForwardAzimuth(compassDirection);
        }
        return ptList;
    }

    /**
     * This method returns a hash code for this Geodetic2DCircle object.
     * Circles that have equal parameters have the same hash code, since
     * they are considered equal objects.
     *
     * @return a hash code value for this object.
     */
    @Override
    public int hashCode() {
        return 31 * center.hashCode() + Double.valueOf(radius).hashCode();
    }

    /**
     * This method is used to compare this circle to another Geodetic2DCircle
     * object. Circles that have equal parameters are considered equal objects.
     *
     * @param circle Geodetic2DCircle to compare to this one
     * @return boolean indicating if this circle is equal to specified one
     */
    public boolean equals(Geodetic2DCircle circle) {
        return circle != null && this.radius == circle.radius &&
				this.getCenter().equals(circle.getCenter());
    }

    /**
     * This method is used to compare this circle to another Geodetic2DCircle
     * object. Circles that have equal parameters are considered equal objects.
     *
     * @param that Geodetic2DCircle to compare to this one
     * @return boolean indicating if this circle is equal to specified one
     */
    @Override
    public boolean equals(Object that) {
        return that instanceof Geodetic2DCircle &&
                equals((Geodetic2DCircle) that);
    }

    /**
     * The toString method formats the circle for printing.
     *
     * @return String representation of center point + radius
     */
    @Override
    public String toString() {
        return "(within " + radius + " meters of " + center + ")";
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/NotAnMGRSBoxException.java<|end_filename|>
/***************************************************************************
 * $Id$
 * 
 * (C) Copyright MITRE Corporation 2008
 *
 * The program is provided "as is" without any warranty express or implied,
 * including the warranty of non-infringement and the implied warranties of
 * merchantability and fitness for a particular purpose.  The Copyright
 * owner will not be liable for any damages suffered by you as a result of
 * using the Program.  In no event will the Copyright owner be liable for
 * any special, indirect or consequential damages or lost profits even if
 * the Copyright owner has been advised of the possibility of their
 * occurrence.
 *
 ***************************************************************************/

package org.opensextant.geodesy;

/**
 * A specialized exception for MGRS conversions.
 * 
 * @author jgibson
 */
public class NotAnMGRSBoxException extends IllegalArgumentException {
	private static final long serialVersionUID = 1L;

    /**
     * Constructs an instance of {@code NotAnMGRSBoxException} with the
	 * specified detail message.
	 * 
     * @param msg the detail message.
     */
    public NotAnMGRSBoxException(final String msg) {
        super(msg);
    }
}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Topocentric2DPoint.java<|end_filename|>
/****************************************************************************************
 *  Topocentric2DPoint.java
 *
 *  Created: Mar 14, 2007
 *
 *  @author <NAME>, based on work by <NAME>, <NAME>, and <NAME>
 *
 *  (C) Copyright MITRE Corporation 2006
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

import edu.umd.cs.findbugs.annotations.NonNull;

import java.io.Serializable;
import java.text.DecimalFormat;

/**
 * The Topocentric2DPoint class is used to represent a point in Cartesian coordinate
 * space, relative to a plane that is tangent to the surface of an Ellipsoid at a
 * particular Geodetic Latitude and Longitude point, called the topocentric origin.
 * To be correctly interpreted, one needs a particular reference Topocentric origin
 * point and Ellipsoid model of the earth (see FrameOfReference class for ways to
 * define and use this in coordinate conversions). The coordinate system used here is
 * consistent with the one defined by Nato Standards for Airborne Ground Survelliance
 * (AGS) data. In particular, the x-axis is positive in the Easterly direction, the
 * y-axis is positive in the Northly direction. All units of length are specified in
 * meters. Note that there are other commonly used systems of topocentric reference
 * that make different assumptions about the directional orientation of the axes.
 * The x and y values are called easting and northing here, but they may not actually
 * follow the corresponding compass directions, depending on how close you are to the
 * polar regions. For example, the MGRS system uses a PolarStereographic projection
 * for the poles, and the easting and northing values of a topocentric point there
 * would follow the MGRS grid layed over that projection.
 */
public class Topocentric2DPoint implements GeoPoint, Serializable {
	private static final long serialVersionUID = 1L;
	
	protected double easting;     // x-axis
    protected double northing;    // y-axis

    /**
     * The constructor takes double precision floating point values for easting and
     * northing distances in meters from the topographic origin.  The origin itself
     * is a contextual variable set in the FrameOfReference object and used there by
     * the coordinate transform methods that need it.
     *
     * @param easting  number of meters along x-axis from topographic origin
     * @param northing number of meters along y-axis from topographic origin
     */
    public Topocentric2DPoint(double easting, double northing) {
        this.easting = easting;
        this.northing = northing;
    }

    /**
     * This accessor method is used to get the number of meters from the topographic
     * origin, along the x-axis or easting direction.
     *
     * @return number of meters along x-axis from topographic origin (easting)
     */
    public double getEasting() {
        return easting;
    }

    /**
     * This settor method is used to update the number of meters from the topographic
     * origin, along the x-axis or easting direction.
     *
     * @param easting number of meters along x-axis from topographic origin
     */
    public void setEasting(double easting) {
        this.easting = easting;
    }

    /**
     * This accessor method is used to get the number of meters from the topographic
     * origin, along the y-axis or northing direction.
     *
     * @return number of meters along y-axis from topographic origin (northing)
     */
    public double getNorthing() {
        return northing;
    }

    /**
     * This settor method is used to update the number of meters from the topographic
     * origin, along the y-axis or northing direction.
     *
     * @param northing number of meters along y-axis from topographic origin
     */
    public void setNorthing(double northing) {
        this.northing = northing;
    }

    /**
     * This method returns a hash code for this Topocentric2DPoint object. The
     * result is the exclusive OR of the component values to maintain the general
     * contract for the hashCode method, which states that equal objects must have
     * equal hash codes.
     *
     * @return a hash code value for this object.
     */
    public int hashCode() {
        return 31 * Double.valueOf(easting).hashCode() +
                Double.valueOf(northing).hashCode();
    }

    /**
     * This method is used to test whether two points are equal in the sense that have
     * the same coordinate value.
     *
     * @param that Topocentric2DPoint point to compare against this one.
     * @return true if specified Topocentric2DPoint point is equal in value to this
     *         Topocentric2DPoint point.
     */
    public boolean equals(Topocentric2DPoint that) {
        return that != null && this.easting == that.easting &&
				this.northing == that.northing;
    }

    /**
     * This method is used to test whether two points are equal in the sense that have
     * the same coordinate value.
     *
     * @param that Topocentric2DPoint point to compare against this one.
     * @return true if specified Topocentric2DPoint point is equal in value to this
     *         Topocentric2DPoint point.
     */
    public boolean equals(Object that) {
        return that instanceof Topocentric2DPoint && equals((Topocentric2DPoint) that);
    }

    /**
     * This method converts this Topocentric2DPoint to a Topocentric3DPoint, by assuming
     * an elevation (z-axis) value of zero meters from the Topographic point of origin.
     *
     * @return the equivalent Topocentric3DPoint
     */
    @NonNull
    public Topocentric3DPoint toTopocentric3D() {
        return new Topocentric3DPoint(easting, northing, 0.0);
    }

    /**
     * This method abstracts this Topocentric2DPoint object as a general purpose GeoPoint.
     *
     * @param fRef the FrameOfReference in which to interpret this coordinate.
     * @return the equivalent Geodetic3DPoint
     */
    @NonNull
    public Geodetic3DPoint toGeodetic3D(FrameOfReference fRef) {
        return fRef.toGeodetic(this.toTopocentric3D());
    }

    /**
     * This method formats the Topocentric2DPoint point as a parenthezised String containing
     * easting and northing values with the number of fractional digits of precision specified.
     *
     * @param fractDig number of fractional digits to display
     * @return String representation of this coordinate suitable for debugging
     */
    public String toString(int fractDig) {
        String fract = (fractDig > 0) ? "." : "";
        while (fractDig-- > 0) fract += "0";
        DecimalFormat fmt = new DecimalFormat("0" + fract);
        return "(" + fmt.format(easting) + "m East, " + fmt.format(northing) + "m North)";
    }

    /**
     * This method formats the Topocentric2DPoint point as a parenthezised String containing
     * easting and northing values in integer meters.
     *
     * @return String representation of this coordinate suitable for debugging
     */
    public String toString() {
        return toString(0);
    }

}


<|start_filename|>library/src/main/java/org/opensextant/geodesy/Latitude.java<|end_filename|>
/****************************************************************************************
 *  Latitude.java
 *
 *  Created: Dec 28, 2006
 *
 *  @author <NAME>
 *
 *  (C) Copyright MITRE Corporation 2007
 *
 *  The program is provided "as is" without any warranty express or implied, including 
 *  the warranty of non-infringement and the implied warranties of merchantability and 
 *  fitness for a particular purpose.  The Copyright owner will not be liable for any 
 *  damages suffered by you as a result of using the Program.  In no event will the 
 *  Copyright owner be liable for any special, indirect or consequential damages or 
 *  lost profits even if the Copyright owner has been advised of the possibility of 
 *  their occurrence.
 *
 ***************************************************************************************/
package org.opensextant.geodesy;

/**
 * The Latitude class extends Angle to restrict its value to be between -90&deg and +90&deg,
 * and to provide parsing and formatting with hemisphere indicator suffix ('N' or 'S')                 ser
 * instead of numeric sign.
 */
public class Latitude extends Angle {
	
	private static final long serialVersionUID = 1L;

    private static final double PI_OVER_2 = Math.PI / 2.0;

    // Base initializer is extended to perform additional range restriction test
    protected void init(double value, int unitType) {
        super.init(value, unitType);
        // Latitudes are Angles restricted between -90 and +90 degrees
        if (this.inRadians < -PI_OVER_2 || PI_OVER_2 < this.inRadians)
            throw new IllegalArgumentException("Latitude value exceeds pole value");
    }

    /**
     * This constructor takes latitude angle value and units designator.
     *
     * @param value    Latitude angular value in specified units (Angle.DEGREES or Angle.RADIANS)
     * @param unitType Angle.DEGREES or Angle.RADIANS
     * @throws IllegalArgumentException error if units or value are out of range
     */
    public Latitude(double value, int unitType) {
        this.init(value, unitType);
    }

    /**
     * This constructor takes a latitude formatted as a string.  It looks for a 'N' or 'S'
     * direction indicator suffix character and sets the numeric sign accordingly (N is +, S is -).
     * Units are assumed to be in degrees unless marked as radians.
     *
     * @param latStr latitude value string
     * @throws IllegalArgumentException error if syntax is incorrect or value is out of range
     */
    public Latitude(String latStr) {
		if (latStr == null) throw new IllegalArgumentException("latStr cannot be null");
        latStr = latStr.trim().toUpperCase();
        int n = latStr.length() - 1;
		if (n < 0) throw new IllegalArgumentException("Latitude must be non-empty value");
        char suffix = latStr.charAt(n);
        char sign = ' ';
        if (suffix == 'N') sign = '+';
        else if (suffix == 'S') sign = '-';
        else if ((suffix == 'E') || (suffix == 'W'))
            throw new IllegalArgumentException("Latitude with E or W direction indicator ");
        if (sign != ' ') {
            // Make sure there is no numeric sign prefix when direction char suffix is present
            char prefix = latStr.charAt(0);
            if ((prefix == '+') || (prefix == '-'))
                throw new IllegalArgumentException("Latitude with N or S direction indicator " +
                        "should not have numeric sign prefix.");
            // Strip off the suffix and add the sign based on direction indicator
            latStr = sign + latStr.substring(0, n);
        }
        this.init(new Angle(latStr).inRadians(), Angle.RADIANS);
    }

    /**
     * This constructor assumes Latitude angle units are in decimal radians.
     *
     * @param valueInRadians Latitude angular value in radians
     * @throws IllegalArgumentException error if value is out of range (-PI/2 to +PI/2 radians)
     */
    public Latitude(double valueInRadians) {
        this.init(valueInRadians, Angle.RADIANS);
    }

    /**
     * This constructor takes numeric degree, minutes, seconds components for Latitude.
     *
     * @param degrees signed, integer degree component of Latitude angular value
     * @param minutes integer minutes component of Latitude angular value
     * @param seconds double seconds component of Latitude angular value
     * @throws IllegalArgumentException error if a component value is out of range
     */
    public Latitude(int degrees, int minutes, double seconds) {
        this.init(degrees, minutes, seconds);
    }

    /**
     * This constructor takes an Angle object for Latitude.
     *
     * @param angle Angle to use as basis for this Latitude (must be between -90 and +90 degrees)
     */
    public Latitude(Angle angle) {
        this.init(angle.inRadians(), Angle.RADIANS);
    }

    /**
     * The default constructor makes a Latitude object at the Equator (value of zero radians).
     */
    public Latitude() {
        this.inRadians = 0.0;
    }

    /**
     * This accessor method returns the character for the hemisphere of this Latitude ('N' or 'S').
     *
     * @return char indicating southern ('S') or northern ('N') hemisphere
     */
    public char getHemisphere() {
        return (this.inRadians < 0.0) ? 'S' : 'N';
    }

    /**
     * This predicate method determines whether this Latitude is within the test
     * interval specified from southLat to northLat.
     *
     * @param southLat southern most Latitude of the interval being tested
     * @param northLat northern most Latitude of the interval being tested
     * @return boolean indicating whether this Latitude is in the test interval
     */
    public boolean inInterval(Latitude southLat, Latitude northLat) {
        return (southLat.inRadians <= this.inRadians) && (this.inRadians <= northLat.inRadians);
    }

    /**
     * This method returns a string version of this latitude value in degree, minute,
     * seconds format, with hemisphere character suffix.
     */
    public String toString(int fractDigOfSec) {
        String result = super.toString(fractDigOfSec);
        char hemisphere = this.getHemisphere();
        if (result.charAt(0) == '-') result = result.substring(1);
        return result + " " + hemisphere;
    }
}
