45 | | /** |
46 | | * create a new ellipsoid and precompute its parameters |
47 | | * |
48 | | * @param a ellipsoid long axis (in meters) |
49 | | * @param b ellipsoid short axis (in meters) |
50 | | */ |
51 | | public Ellipsoid(double a, double b) { |
52 | | this.a = a; |
53 | | this.b = b; |
54 | | e2 = (a*a - b*b) / (a*a); |
55 | | e = Math.sqrt(e2); |
56 | | } |
| 45 | /** |
| 46 | * square of the second eccentricity |
| 47 | */ |
| 48 | public final double eb2; |
| 49 | |
| 50 | /** |
| 51 | * create a new ellipsoid and precompute its parameters |
| 52 | * |
| 53 | * @param a ellipsoid long axis (in meters) |
| 54 | * @param b ellipsoid short axis (in meters) |
| 55 | */ |
| 56 | public Ellipsoid(double a, double b) { |
| 57 | this.a = a; |
| 58 | this.b = b; |
| 59 | e2 = (a*a - b*b) / (a*a); |
| 60 | e = Math.sqrt(e2); |
| 61 | eb2 = e2 / (1.0 - e2); |
| 62 | } |
| 63 | |
| 64 | /** |
| 65 | * Returns the <i>radius of curvature in the prime vertical</i> |
| 66 | * for this reference ellipsoid at the specified latitude. |
| 67 | * |
| 68 | * @param phi The local latitude (radians). |
| 69 | * @return The radius of curvature in the prime vertical (meters). |
| 70 | */ |
| 71 | public double verticalRadiusOfCurvature(final double phi) { |
| 72 | return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi)))); |
| 73 | } |
| 74 | |
| 75 | private static double sqr(final double x) { |
| 76 | return x * x; |
| 77 | } |
| 78 | |
| 79 | /** |
| 80 | * Returns the meridional arc, the true meridional distance on the |
| 81 | * ellipsoid from the equator to the specified latitude, in meters. |
| 82 | * |
| 83 | * @param phi The local latitude (in radians). |
| 84 | * @return The meridional arc (in meters). |
| 85 | */ |
| 86 | public double meridionalArc(final double phi) { |
| 87 | final double sin2Phi = Math.sin(2.0 * phi); |
| 88 | final double sin4Phi = Math.sin(4.0 * phi); |
| 89 | final double sin6Phi = Math.sin(6.0 * phi); |
| 90 | final double sin8Phi = Math.sin(8.0 * phi); |
| 91 | // TODO . calculate 'f' |
| 92 | //double f = 1.0 / 298.257222101; // GRS80 |
| 93 | double f = 1.0 / 298.257223563; // WGS84 |
| 94 | final double n = f / (2.0 - f); |
| 95 | final double n2 = n * n; |
| 96 | final double n3 = n2 * n; |
| 97 | final double n4 = n3 * n; |
| 98 | final double n5 = n4 * n; |
| 99 | final double n1n2 = n - n2; |
| 100 | final double n2n3 = n2 - n3; |
| 101 | final double n3n4 = n3 - n4; |
| 102 | final double n4n5 = n4 - n5; |
| 103 | final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5)); |
| 104 | final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5); |
| 105 | final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5)); |
| 106 | final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5); |
| 107 | final double ep = (315.0 / 512.0) * a * (n4n5); |
| 108 | return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi; |
| 109 | } |
| 110 | |
| 111 | /** |
| 112 | * Returns the <i>radius of curvature in the meridian<i> |
| 113 | * for this reference ellipsoid at the specified latitude. |
| 114 | * |
| 115 | * @param phi The local latitude (in radians). |
| 116 | * @return The radius of curvature in the meridian (in meters). |
| 117 | */ |
| 118 | public double meridionalRadiusOfCurvature(final double phi) { |
| 119 | return verticalRadiusOfCurvature(phi) |
| 120 | / (1.0 + eb2 * sqr(Math.cos(phi))); |
| 121 | } |
| 122 | |
| 123 | /** |
| 124 | * Returns isometric latitude of phi on given first eccentricity (e) |
| 125 | * @param phi The local latitude (radians). |
| 126 | * @return isometric latitude of phi on first eccentricity (e) |
| 127 | */ |
| 128 | public double latitudeIsometric(double phi, double e) { |
| 129 | double v1 = 1-e*Math.sin(phi); |
| 130 | double v2 = 1+e*Math.sin(phi); |
| 131 | return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2)); |
| 132 | } |
| 133 | |
| 134 | /** |
| 135 | * Returns isometric latitude of phi on first eccentricity (e) |
| 136 | * @param phi The local latitude (radians). |
| 137 | * @return isometric latitude of phi on first eccentricity (e) |
| 138 | */ |
| 139 | public double latitudeIsometric(double phi) { |
| 140 | double v1 = 1-e*Math.sin(phi); |
| 141 | double v2 = 1+e*Math.sin(phi); |
| 142 | return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2)); |
| 143 | } |
| 144 | |
| 145 | /* |
| 146 | * Returns geographic latitude of isometric latitude of first eccentricity (e) |
| 147 | * and epsilon precision |
| 148 | */ |
| 149 | public double latitude(double latIso, double e, double epsilon) { |
| 150 | double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2; |
| 151 | double lati = lat0; |
| 152 | double lati1 = 1.0; // random value to start the iterative processus |
| 153 | while(Math.abs(lati1-lati)>=epsilon) { |
| 154 | lati = lati1; |
| 155 | double v1 = 1+e*Math.sin(lati); |
| 156 | double v2 = 1-e*Math.sin(lati); |
| 157 | lati1 = 2*Math.atan(Math.pow(v1/v2,e/2)*Math.exp(latIso))-Math.PI/2; |
| 158 | } |
| 159 | return lati1; |
| 160 | } |