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 extends AbstractProj { 032 033 private Ellipsoid ellps; 034 private double n; 035 private double c; 036 private double chi0; 037 private double r; 038 039 private static final double EPSILON = 1e-12; 040 041 @Override 042 public String getName() { 043 return tr("Double Stereographic"); 044 } 045 046 @Override 047 public String getProj4Id() { 048 return "sterea"; 049 } 050 051 @Override 052 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 053 super.initialize(params); 054 if (params.lat0 == null) 055 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 056 ellps = params.ellps; 057 initialize(params.lat0); 058 } 059 060 private void initialize(double lat0) { 061 double phi0 = toRadians(lat0); 062 double e2 = ellps.e2; 063 r = sqrt(1-e2) / (1 - e2*pow(sin(phi0), 2)); 064 n = sqrt(1 + ellps.eb2 * pow(cos(phi0), 4)); 065 double s1 = (1 + sin(phi0)) / (1 - sin(phi0)); 066 double s2 = (1 - e * sin(phi0)) / (1 + e * sin(phi0)); 067 double w1 = pow(s1 * pow(s2, e), n); 068 double sinchi00 = (w1 - 1) / (w1 + 1); 069 c = (n + sin(phi0)) * (1 - sinchi00) / ((n - sin(phi0)) * (1 + sinchi00)); 070 double w2 = c * w1; 071 chi0 = asin((w2 - 1) / (w2 + 1)); 072 } 073 074 @Override 075 public double[] project(double phi, double lambda) { 076 double nLambda = n * lambda; 077 double sa = (1 + sin(phi)) / (1 - sin(phi)); 078 double sb = (1 - e * sin(phi)) / (1 + e * sin(phi)); 079 double w = c * pow(sa * pow(sb, e), n); 080 double chi = asin((w - 1) / (w + 1)); 081 double b = 1 + sin(chi) * sin(chi0) + cos(chi) * cos(chi0) * cos(nLambda); 082 double x = 2 * r * cos(chi) * sin(nLambda) / b; 083 double y = 2 * r * (sin(chi) * cos(chi0) - cos(chi) * sin(chi0) * cos(nLambda)) / b; 084 return new double[] {x, y}; 085 } 086 087 @Override 088 public double[] invproject(double x, double y) { 089 double e2 = ellps.e2; 090 double g = 2 * r * tan(PI/4 - chi0/2); 091 double h = 4 * r * tan(chi0) + g; 092 double i = atan(x/(h + y)); 093 double j = atan(x/(g - y)) - i; 094 double chi = chi0 + 2 * atan((y - x * tan(j/2)) / (2 * r)); 095 double lambda = (j + 2*i) / n; 096 double psi = 0.5 * log((1 + sin(chi)) / (c*(1 - sin(chi)))) / n; 097 double phiprev = -1000; 098 int iteration = 0; 099 double phi = 2 * atan(exp(psi)) - PI/2; 100 while (abs(phi - phiprev) > EPSILON) { 101 if (++iteration > 10) 102 throw new IllegalStateException("Too many iterations"); 103 phiprev = phi; 104 double psii = log(tan(phi/2 + PI/4) * pow((1 - e * sin(phi)) / (1 + e * sin(phi)), e/2)); 105 phi = phi - (psii - psi) * cos(phi) * (1 - e2 * pow(sin(phi), 2)) / (1 - e2); 106 } 107 return new double[] {phi, lambda}; 108 } 109 110 @Override 111 public Bounds getAlgorithmBounds() { 112 return new Bounds(-89, -87, 89, 87, false); 113 } 114}