Index: Ellipsoid.java
===================================================================
--- Ellipsoid.java (revision 1919)
+++ Ellipsoid.java (working copy)
@@ -11,49 +11,154 @@
* the reference ellipsoids
*/
class Ellipsoid {
- /**
- * Clarke's elipsoid (NTF system)
- */
- public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
- /**
- * Hayford's ellipsoid (ED50 system)
- */
- public static final Ellipsoid hayford =
- new Ellipsoid(6378388.0, 6356911.9461);
- /**
- * WGS84 elipsoid
- */
- public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.314);
+ /**
+ * Clarke's ellipsoid (NTF system)
+ */
+ public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
+ /**
+ * Hayford's ellipsoid (ED50 system)
+ */
+ public static final Ellipsoid hayford =
+ new Ellipsoid(6378388.0, 6356911.9461);
+ /**
+ * WGS84 ellipsoid
+ */
+ public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.314);
- /**
- * half long axis
- */
- public final double a;
- /**
- * half short axis
- */
- public final double b;
- /**
- * first excentricity
- */
- public final double e;
- /**
- * first excentricity squared
- */
- public final double e2;
+ /**
+ * half long axis
+ */
+ public final double a;
+ /**
+ * half short axis
+ */
+ public final double b;
+ /**
+ * first eccentricity
+ */
+ public final double e;
+ /**
+ * first eccentricity squared
+ */
+ public final double e2;
- /**
- * create a new ellipsoid and precompute its parameters
- *
- * @param a ellipsoid long axis (in meters)
- * @param b ellipsoid short axis (in meters)
- */
- public Ellipsoid(double a, double b) {
- this.a = a;
- this.b = b;
- e2 = (a*a - b*b) / (a*a);
- e = Math.sqrt(e2);
- }
+ /**
+ * square of the second eccentricity
+ */
+ public final double eb2;
+
+ /**
+ * create a new ellipsoid and precompute its parameters
+ *
+ * @param a ellipsoid long axis (in meters)
+ * @param b ellipsoid short axis (in meters)
+ */
+ public Ellipsoid(double a, double b) {
+ this.a = a;
+ this.b = b;
+ e2 = (a*a - b*b) / (a*a);
+ e = Math.sqrt(e2);
+ eb2 = e2 / (1.0 - e2);
+ }
+
+ /**
+ * Returns the radius of curvature in the prime vertical
+ * for this reference ellipsoid at the specified latitude.
+ *
+ * @param phi The local latitude (radians).
+ * @return The radius of curvature in the prime vertical (meters).
+ */
+ public double verticalRadiusOfCurvature(final double phi) {
+ return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
+ }
+
+ private static double sqr(final double x) {
+ return x * x;
+ }
+
+ /**
+ * Returns the meridional arc, the true meridional distance on the
+ * ellipsoid from the equator to the specified latitude, in meters.
+ *
+ * @param phi The local latitude (in radians).
+ * @return The meridional arc (in meters).
+ */
+ public double meridionalArc(final double phi) {
+ final double sin2Phi = Math.sin(2.0 * phi);
+ final double sin4Phi = Math.sin(4.0 * phi);
+ final double sin6Phi = Math.sin(6.0 * phi);
+ final double sin8Phi = Math.sin(8.0 * phi);
+ // TODO . calculate 'f'
+ //double f = 1.0 / 298.257222101; // GRS80
+ double f = 1.0 / 298.257223563; // WGS84
+ final double n = f / (2.0 - f);
+ final double n2 = n * n;
+ final double n3 = n2 * n;
+ final double n4 = n3 * n;
+ final double n5 = n4 * n;
+ final double n1n2 = n - n2;
+ final double n2n3 = n2 - n3;
+ final double n3n4 = n3 - n4;
+ final double n4n5 = n4 - n5;
+ final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
+ final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
+ final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
+ final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
+ final double ep = (315.0 / 512.0) * a * (n4n5);
+ return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
+ }
+
+ /**
+ * Returns the radius of curvature in the meridian
+ * for this reference ellipsoid at the specified latitude.
+ *
+ * @param phi The local latitude (in radians).
+ * @return The radius of curvature in the meridian (in meters).
+ */
+ public double meridionalRadiusOfCurvature(final double phi) {
+ return verticalRadiusOfCurvature(phi)
+ / (1.0 + eb2 * sqr(Math.cos(phi)));
+ }
+
+ /**
+ * Returns isometric latitude of phi on given first eccentricity (e)
+ * @param phi The local latitude (radians).
+ * @return isometric latitude of phi on first eccentricity (e)
+ */
+ public double latitudeIsometric(double phi, double e) {
+ double v1 = 1-e*Math.sin(phi);
+ double v2 = 1+e*Math.sin(phi);
+ return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
+ }
+
+ /**
+ * Returns isometric latitude of phi on first eccentricity (e)
+ * @param phi The local latitude (radians).
+ * @return isometric latitude of phi on first eccentricity (e)
+ */
+ public double latitudeIsometric(double phi) {
+ double v1 = 1-e*Math.sin(phi);
+ double v2 = 1+e*Math.sin(phi);
+ return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
+ }
+
+ /*
+ * Returns geographic latitude of isometric latitude of first eccentricity (e)
+ * and epsilon precision
+ */
+ public double latitude(double latIso, double e, double epsilon) {
+ double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
+ double lati = lat0;
+ double lati1 = 1.0; // random value to start the iterative processus
+ while(Math.abs(lati1-lati)>=epsilon) {
+ lati = lati1;
+ double v1 = 1+e*Math.sin(lati);
+ double v2 = 1-e*Math.sin(lati);
+ lati1 = 2*Math.atan(Math.pow(v1/v2,e/2)*Math.exp(latIso))-Math.PI/2;
+ }
+ return lati1;
+ }
}
+