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