001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.data.projection.proj; 003 004import static java.lang.Math.PI; 005import static java.lang.Math.abs; 006import static java.lang.Math.asin; 007import static java.lang.Math.atan; 008import static java.lang.Math.cos; 009import static java.lang.Math.exp; 010import static java.lang.Math.log; 011import static java.lang.Math.pow; 012import static java.lang.Math.sin; 013import static java.lang.Math.sqrt; 014import static java.lang.Math.tan; 015import static java.lang.Math.toRadians; 016import static org.openstreetmap.josm.tools.I18n.tr; 017 018import org.openstreetmap.josm.data.Bounds; 019import org.openstreetmap.josm.data.projection.Ellipsoid; 020import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 021 022/** 023 * Implementation of the stereographic double projection, 024 * also known as Oblique Stereographic and the Schreiber double stereographic projection. 025 * 026 * @author vholten 027 * 028 * Source: IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2, 029 * Sec. 1.3.7.1 Oblique and Equatorial Stereographic, http://www.epsg.org/GuidanceNotes 030 */ 031public class DoubleStereographic implements Proj { 032 033 private Ellipsoid ellps; 034 private double e; 035 private double n; 036 private double c; 037 private double chi0; 038 private double R; 039 040 private static final double EPSILON = 1e-12; 041 042 @Override 043 public String getName() { 044 return tr("Double Stereographic"); 045 } 046 047 @Override 048 public String getProj4Id() { 049 return "sterea"; 050 } 051 052 @Override 053 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 054 if (params.lat0 == null) 055 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 056 ellps = params.ellps; 057 e = ellps.e; 058 initialize(params.lat0); 059 } 060 061 private void initialize(double lat_0) { 062 double phi0 = toRadians(lat_0); 063 double e2 = ellps.e2; 064 R = sqrt(1-e2) / (1 - e2*pow(sin(phi0), 2)); 065 n = sqrt(1 + ellps.eb2 * pow(cos(phi0), 4)); 066 double S1 = (1 + sin(phi0)) / (1 - sin(phi0)); 067 double S2 = (1 - e * sin(phi0)) / (1 + e * sin(phi0)); 068 double w1 = pow(S1 * pow(S2, e), n); 069 double sinchi00 = (w1 - 1) / (w1 + 1); 070 c = (n + sin(phi0)) * (1 - sinchi00) / ((n - sin(phi0)) * (1 + sinchi00)); 071 double w2 = c * w1; 072 chi0 = asin((w2 - 1) / (w2 + 1)); 073 } 074 075 @Override 076 public double[] project(double phi, double lambda) { 077 double Lambda = n * lambda; 078 double Sa = (1 + sin(phi)) / (1 - sin(phi)); 079 double Sb = (1 - e * sin(phi)) / (1 + e * sin(phi)); 080 double w = c * pow(Sa * pow(Sb, e), n); 081 double chi = asin((w - 1) / (w + 1)); 082 double B = 1 + sin(chi) * sin(chi0) + cos(chi) * cos(chi0) * cos(Lambda); 083 double x = 2 * R * cos(chi) * sin(Lambda) / B; 084 double y = 2 * R * (sin(chi) * cos(chi0) - cos(chi) * sin(chi0) * cos(Lambda)) / B; 085 return new double[] {x, y}; 086 } 087 088 @Override 089 public double[] invproject(double x, double y) { 090 double e2 = ellps.e2; 091 double g = 2 * R * tan(PI/4 - chi0/2); 092 double h = 4 * R * tan(chi0) + g; 093 double i = atan(x/(h + y)); 094 double j = atan(x/(g - y)) - i; 095 double chi = chi0 + 2 * atan((y - x * tan(j/2)) / (2 * R)); 096 double Lambda = j + 2*i; 097 double lambda = Lambda / n; 098 double psi = 0.5 * log((1 + sin(chi)) / (c*(1 - sin(chi)))) / n; 099 double phiprev = -1000; 100 int iteration = 0; 101 double phi = 2 * atan(exp(psi)) - PI/2; 102 while (abs(phi - phiprev) > EPSILON) { 103 if (++iteration > 10) 104 throw new RuntimeException("Too many iterations"); 105 phiprev = phi; 106 double psii = log(tan(phi/2 + PI/4) * pow((1 - e * sin(phi)) / (1 + e * sin(phi)), e/2)); 107 phi = phi - (psii - psi) * cos(phi) * (1 - e2 * pow(sin(phi), 2)) / (1 - e2); 108 } 109 return new double[] {phi, lambda}; 110 } 111 112 @Override 113 public Bounds getAlgorithmBounds() { 114 return new Bounds(-89, -87, 89, 87, false); 115 } 116}