001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.data.projection.proj; 003 004import static org.openstreetmap.josm.tools.I18n.tr; 005 006import org.openstreetmap.josm.data.Bounds; 007import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 008import org.openstreetmap.josm.tools.Utils; 009 010/** 011 * Cassini-Soldner Projection (EPSG code 9806). 012 * The Cassini-Soldner Projection is the ellipsoidal version of the Cassini 013 * projection for the sphere. It is not conformal but as it is relatively simple 014 * to construct it was extensively used in the last century and is still useful 015 * for mapping areas with limited longitudinal extent. It has now largely 016 * been replaced by the conformal Transverse Mercator which it resembles. Like this, 017 * it has a straight central meridian along which the scale is true, all other 018 * meridians and parallels are curved, and the scale distortion increases 019 * rapidly with increasing distance from the central meridian. 020 * <p> 021 * 022 * This class has been derived from the implementation of the Geotools project; 023 * git 8cbf52d, org.geotools.referencing.operation.projection.CassiniSoldner 024 * at the time of migration. 025 */ 026public class CassiniSoldner extends AbstractProj { 027 028 /** 029 * Meridian distance at the {@code latitudeOfOrigin}. 030 * Used for calculations for the ellipsoid. 031 */ 032 private double ml0; 033 034 /** 035 * Latitude of origin. 036 */ 037 private double phi0; 038 039 /** 040 * Constants used for the forward and inverse transform for the elliptical 041 * case of the Cassini-Soldner. 042 */ 043 private static final double C1 = 1. / 6; 044 private static final double C2 = 1. / 120; 045 private static final double C3 = 1. / 24; 046 private static final double C4 = 1. / 3; 047 private static final double C5 = 1. / 15; 048 049 @Override 050 public String getName() { 051 return tr("Cassini-Soldner"); 052 } 053 054 @Override 055 public String getProj4Id() { 056 return "cass"; 057 } 058 059 @Override 060 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 061 super.initialize(params); 062 if (params.lat0 == null) 063 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 064 phi0 = Utils.toRadians(params.lat0); 065 ml0 = mlfn(phi0, Math.sin(phi0), Math.cos(phi0)); 066 } 067 068 @Override 069 public double[] project(double phi, double lam) { 070 if (spherical) { 071 double x = aasin(Math.cos(phi) * Math.sin(lam)); 072 double y = Math.atan2(Math.tan(phi), Math.cos(lam)); 073 return new double[] {x, y}; 074 } else { 075 double sinphi = Math.sin(phi); 076 double cosphi = Math.cos(phi); 077 078 double n = 1.0 / Math.sqrt(1.0 - e2 * sinphi * sinphi); 079 double tn = Math.tan(phi); 080 double t = tn * tn; 081 double a1 = lam * cosphi; 082 double c = cosphi * cosphi * e2 / (1 - e2); 083 double a2 = a1 * a1; 084 085 double x = n * a1 * (1.0 - a2 * t * (C1 - (8.0 - t + 8.0 * c) * a2 * C2)); 086 double y = mlfn(phi, sinphi, cosphi) - ml0 + n * tn * a2 * (0.5 + (5.0 - t + 6.0 * c) * a2 * C3); 087 return new double[] {x, y}; 088 } 089 } 090 091 @Override 092 public double[] invproject(double x, double y) { 093 if (spherical) { 094 double dd = y + phi0; 095 double phi = aasin(Math.sin(dd * Math.cos(x))); 096 double lam = Math.atan2(Math.tan(x), Math.cos(dd)); 097 return new double[] {phi, lam}; 098 } else { 099 double ph1 = invMlfn(ml0 + y); 100 double tn = Math.tan(ph1); 101 double t = tn * tn; 102 double n = Math.sin(ph1); 103 double r = 1.0 / (1.0 - e2 * n * n); 104 n = Math.sqrt(r); 105 r *= (1.0 - e2) * n; 106 double dd = x / n; 107 double d2 = dd * dd; 108 double phi = ph1 - (n * tn / r) * d2 * (0.5 - (1.0 + 3.0 * t) * d2 * C3); 109 double lam = dd * (1.0 + t * d2 * (-C4 + (1.0 + 3.0 * t) * d2 * C5)) / Math.cos(ph1); 110 return new double[] {phi, lam}; 111 } 112 } 113 114 @Override 115 public Bounds getAlgorithmBounds() { 116 return new Bounds(-89, -1.0, 89, 1.0, false); 117 } 118}