001/*
002 * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
003 * Copyright (C) 2006 - JScience (http://jscience.org/)
004 * All rights reserved.
005 *
006 * Permission to use, copy, modify, and distribute this software is
007 * freely granted, provided that this notice is preserved.
008 */
009package org.jscience.geography.coordinates;
010
011import javax.measure.Measure;
012import javax.measure.converter.UnitConverter;
013import javax.measure.quantity.*;
014import static javax.measure.unit.SI.*;
015import javax.measure.unit.Unit;
016
017import javolution.context.ObjectFactory;
018import javolution.xml.XMLFormat;
019import javolution.xml.stream.XMLStreamException;
020
021import static org.jscience.geography.coordinates.crs.ReferenceEllipsoid.WGS84;
022
023import org.jscience.geography.coordinates.crs.GeocentricCRS;
024import org.jscience.mathematics.vector.DimensionException;
025import org.jscience.mathematics.vector.Float64Vector;
026import org.opengis.referencing.cs.CoordinateSystem;
027
028/**
029 * This class represents the {@link GeocentricCRS geocentric} Earth-Centered, 
030 * Earth-Fixed (ECEF) cartesian coordinates used in GPS/GLONASS.
031 *
032 * @author Paul D. Anderson
033 * @version 3.0, February 18, 2006
034 */
035public final class XYZ extends Coordinates<GeocentricCRS<XYZ>> {
036
037    /**
038     * Holds the coordinate reference system for all instances of this class.
039     */
040    public static final GeocentricCRS<XYZ> CRS = new GeocentricCRS<XYZ>() {
041
042        @Override
043        protected XYZ coordinatesOf(AbsolutePosition position) {
044            double latitude = position.latitudeWGS84.doubleValue(RADIAN);
045            double longitude = position.longitudeWGS84.doubleValue(RADIAN);
046            double height = (position.heightWGS84 != null) ?
047                position.heightWGS84.doubleValue(METRE) : 0.0;
048
049            double cosLat = Math.cos(latitude);
050            double sinLat = Math.sin(latitude);
051            double cosLon = Math.cos(longitude);
052            double sinLon = Math.sin(longitude);
053
054            double roc = WGS84.verticalRadiusOfCurvature(latitude);
055            double x = (roc + height) * cosLat * cosLon;
056            double y = (roc + height) * cosLat * sinLon;
057            double z = ((1.0 - WGS84.getEccentricitySquared()) * roc + height)
058                    * sinLat;
059
060            return XYZ.valueOf(x, y, z, METRE);
061        }
062
063        @Override
064        protected AbsolutePosition positionOf(XYZ coordinates,
065                AbsolutePosition position) {
066            final double x = coordinates._x;
067            final double y = coordinates._y;
068            final double z = coordinates._z;
069
070            final double longitude = Math.atan2(y, x);
071
072            final double latitude;
073            final double xy = Math.hypot(x, y);
074            // conventional result if xy == 0.0...
075            if (xy == 0.0) {
076                latitude = (z >= 0.0) ? Math.PI / 2.0 : -Math.PI / 2.0;
077            } else {
078                final double a = WGS84.getSemimajorAxis().doubleValue(METRE);
079                final double b = WGS84.getsSemiminorAxis().doubleValue(METRE);
080                final double ea2 = WGS84.getEccentricitySquared();
081                final double eb2 = WGS84.getSecondEccentricitySquared();
082                final double beta = Math.atan2(a * z, b * xy);
083                double numerator = z + b * eb2 * cube(Math.sin(beta));
084                double denominator = xy - a * ea2 * cube(Math.cos(beta));
085                latitude = Math.atan2(numerator, denominator);
086            }
087
088            final double height = xy / Math.cos(latitude)
089                    - WGS84.verticalRadiusOfCurvature(latitude);
090            position.latitudeWGS84 = Measure.valueOf(latitude, RADIAN);
091            position.longitudeWGS84 = Measure.valueOf(longitude, RADIAN);
092            position.heightWGS84 = Measure.valueOf(height, METRE);
093            return position;
094        }
095
096        @Override
097        public CoordinateSystem getCoordinateSystem() {
098            return GeocentricCRS.XYZ_CS;
099        }
100    };
101
102    private static double cube(final double x) {
103        return x * x * x;
104    }
105
106    /**
107     * Holds the x position in meters.
108     */
109    private double _x;
110
111    /**
112     * Holds the y position in meters.
113     */
114    private double _y;
115
116    /**
117     * Holds the z position in meters.
118     */
119    private double _z;
120
121    /**
122     * Returns the spatial position corresponding to the specified coordinates.
123     *
124     * @param x the x value stated in the specified unit.
125     * @param y the y value stated in the specified unit.
126     * @param z the z value stated in the specified unit.
127     * @param unit the length unit in which the coordinates are stated.
128     * @return the corresponding 3D position.
129     */
130    public static XYZ valueOf(double x, double y, double z, Unit<Length> unit) {
131        XYZ xyz = FACTORY.object();
132        if (unit == METRE) {
133            xyz._x = x;
134            xyz._y = y;
135            xyz._z = z;
136        } else {
137            UnitConverter toMeter = unit.getConverterTo(METRE);
138            xyz._x = toMeter.convert(x);
139            xyz._y = toMeter.convert(y);
140            xyz._z = toMeter.convert(z);
141        }
142        return xyz;
143    }
144
145    private static final ObjectFactory<XYZ> FACTORY = new ObjectFactory<XYZ>() {
146
147        @Override
148        protected XYZ create() {
149            return new XYZ();
150        }
151    };
152
153    private XYZ() {
154    }
155
156    /**
157     * Returns the spatial position corresponding to the specified 
158     * 3-dimensional vector.
159     *
160     * @param vector the 3-dimensional vector holding the x/y/z coordinates.
161     * @param unit the length unit in which the coordinates are stated.
162     * @return the corresponding 3D position.
163     */
164    public static XYZ valueOf(Float64Vector vector, Unit<Length> unit) {
165        if (vector.getDimension() != 3)
166            throw new DimensionException("3-dimensional vector expected");
167        return XYZ.valueOf(vector.getValue(0), vector.getValue(1), vector
168                .getValue(2), unit);
169    }
170
171    /**
172     * Returns the x coordinate value as <code>double</code>
173     *
174     * @param unit the length unit of the x coordinate value to return.
175     * @return the x coordinate stated in the specified unit.
176     */
177    public double xValue(Unit<Length> unit) {
178        return (unit == METRE) ? _x : METRE.getConverterTo(unit).convert(_x);
179    }
180
181    /**
182     * Returns the y coordinate value as <code>double</code>
183     *
184     * @param unit the length unit of the x coordinate value to return.
185     * @return the z coordinate stated in the specified unit.
186     */
187    public double yValue(Unit<Length> unit) {
188        return (unit == METRE) ? _y : METRE.getConverterTo(unit).convert(_y);
189    }
190
191    /**
192     * Returns the z coordinate value as <code>double</code>
193     *
194     * @param unit the length unit of the x coordinate value to return.
195     * @return the z coordinate stated in the specified unit.
196     */
197    public double zValue(Unit<Length> unit) {
198        return (unit == METRE) ? _z : METRE.getConverterTo(unit).convert(_z);
199    }
200
201    /**
202     * Returns the x/y/z coordinates value as a 3-dimensional vector.
203     *
204     * @param unit the length unit of the vector coordinates.
205     * @return a vector holding the x/y/z coordinates stated in the 
206     *         specified unit.
207     */
208    public Float64Vector toVector(Unit<Length> unit) {
209        if (unit == METRE)
210            return Float64Vector.valueOf(_x, _y, _z);
211        UnitConverter cvtr = METRE.getConverterTo(unit);
212        return Float64Vector.valueOf(cvtr.convert(_x), cvtr.convert(_y), cvtr
213                .convert(_z));
214    }
215
216    @Override
217    public GeocentricCRS<XYZ> getCoordinateReferenceSystem() {
218        return CRS;
219    }
220
221    // OpenGIS Interface.
222    public int getDimension() {
223        return 3;
224    }
225
226    // OpenGIS Interface.
227    public double getOrdinate(int dimension) throws IndexOutOfBoundsException {
228        if (dimension == 0) {
229            Unit<?> u = GeocentricCRS.XYZ_CS.getAxis(0).getUnit();
230            return METRE.getConverterTo(u).convert(_x);
231        } else if (dimension == 1) {
232            Unit<?> u = GeocentricCRS.XYZ_CS.getAxis(1).getUnit();
233            return METRE.getConverterTo(u).convert(_y);
234        } else if (dimension == 2) {
235            Unit<?> u = GeocentricCRS.XYZ_CS.getAxis(2).getUnit();
236            return METRE.getConverterTo(u).convert(_z);
237        } else {
238            throw new IndexOutOfBoundsException();
239        }
240    }
241
242    @Override
243    public XYZ copy() {
244        return XYZ.valueOf(_x, _y, _z, METRE);
245    }
246
247    // Default serialization.
248    //
249
250    static final XMLFormat<XYZ> XML = new XMLFormat<XYZ>(XYZ.class) {
251
252        @Override
253        public XYZ newInstance(Class<XYZ> cls, InputElement xml)
254                throws XMLStreamException {
255            return FACTORY.object();
256        }
257
258        public void write(XYZ xyz, OutputElement xml) throws XMLStreamException {
259            xml.setAttribute("x", xyz._x);
260            xml.setAttribute("y", xyz._y);
261            xml.setAttribute("z", xyz._z);
262        }
263
264        public void read(InputElement xml, XYZ xyz) throws XMLStreamException {
265            xyz._x = xml.getAttribute("x", 0.0);
266            xyz._y = xml.getAttribute("y", 0.0);
267            xyz._z = xml.getAttribute("z", 0.0);
268        }
269    };
270
271    private static final long serialVersionUID = 1L;
272
273}