Changeset 4285 in josm

Timestamp:

2011-08-07T12:17:39+02:00 (13 years ago)

Author:

bastiK

Message:

major projection rework

More modular structure, inspired by Proj.4.

There are almost no semantic changes to the projection algorithms. Mostly factors of 'a' and 180/PI have been moved from one place to the other. In UTM_France_DOM, the ellipsoid conversion for the first 3 projections has been changed from hayford <-> GRS80 to hayford <-> WGS84.

Some redundant algorithms have been removed. In particular:

• UTM_France_DOM used to have its own Transverse Mercator implementation. It is different from the implementation in TransverseMercator.java as it has another series expansion. For EPSG::2975, there are numeric differences on centimeter scale. However, the new data fits better with Proj.4 output.
• Also removed are alternate implementations of LambertConformalConic. (They are all quite similar, though.)

Location:

trunk/src/org/openstreetmap/josm/data/projection

Files:

• AbstractProjection.java (added)
• Ellipsoid.java (modified) (1 diff)
• Epsg3008.java (modified) (4 diffs)
• Lambert.java (modified) (11 diffs)
• LambertCC9Zones.java (modified) (9 diffs)
• LambertEST.java (modified) (5 diffs)
• Mercator.java (modified) (4 diffs)
• Puwg.java (modified) (14 diffs)
• SwissGrid.java (modified) (3 diffs)
• TransverseMercatorLV.java (modified) (3 diffs)
• UTM.java (modified) (8 diffs)
• UTM_France_DOM.java (modified) (8 diffs)
• datum (added)
• datum/AbstractDatum.java (added)
• datum/CentricDatum.java (added)
• datum/Datum.java (added)
• datum/GRS80Datum.java (added)
• datum/SevenParameterDatum.java (added)
• datum/ThreeParameterDatum.java (added)
• datum/WGS84Datum.java (added)
• proj (added)
• proj/LambertConformalConic.java (added)
• proj/Mercator.java (added)
• proj/Proj.java (added)
• proj/SwissObliqueMercator.java (added)
• proj/TransverseMercator.java (moved) (moved from trunk/src/org/openstreetmap/josm/data/projection/TransverseMercator.java ) (7 diffs)

trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java

@@ -19 +19 @@
     public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
     /**
-     * Hayford's ellipsoid (ED50 system)
+     * Hayford's ellipsoid 1909 (ED50 system)
+     * Proj.4 code: intl
      */
     public static final Ellipsoid hayford =

trunk/src/org/openstreetmap/josm/data/projection/Epsg3008.java

@@ -6 +6 @@
 import org.openstreetmap.josm.data.Bounds;
 import org.openstreetmap.josm.data.coor.LatLon;
-import org.openstreetmap.josm.data.coor.EastNorth;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
+import org.openstreetmap.josm.data.projection.proj.TransverseMercator;
 
 /**
@@ -12 +13 @@
  * http://spatialreference.org/ref/epsg/3008/
  *
- * @author Hanno Hecker, based on the TransverseMercatorLV.java by Viesturs Zarins
+ * @author Hanno Hecker
  */
-public class Epsg3008 extends TransverseMercator {
+public class Epsg3008 extends AbstractProjection {
 
-    private final static double UTMScaleFactor = 1.0;
-    private double UTMCentralMeridianRad;
-    private double offsetEastMeters = 150000;
-    private double offsetNorthMeters = 0;
-
-    public Epsg3008()
-    {
-        UTMCentralMeridianRad = Math.toRadians(13.5);
+    public Epsg3008() {
+        ellps = Ellipsoid.GRS80;
+        proj = new TransverseMercator(ellps);
+        datum = GRS80Datum.INSTANCE;
+        lon_0 = 13.5;
+        x_0 = 150000;
     }
-
-    @Override public String toString() {
+    
+    @Override
+    public String toString() {
         return tr("SWEREF99 13 30 / EPSG:3008 (Sweden)");
     }
 
-    private int epsgCode() {
+    @Override
+    public Integer getEpsgCode() {
         return 3008;
-    }
-
-    @Override
-    public String toCode() {
-        return "EPSG:"+ epsgCode();
     }
 
@@ -44 +40 @@
     }
 
+    @Override
     public String getCacheDirectoryName() {
-        return "epsg"+ epsgCode();
+        return "epsg"+ getEpsgCode();
     }
 
@@ -55 +52 @@
     }
 
-    @Override
-    public EastNorth latlon2eastNorth(LatLon p) {
-        EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), UTMCentralMeridianRad);
-        return new EastNorth(a.east() * UTMScaleFactor + offsetEastMeters, a.north() * UTMScaleFactor + offsetNorthMeters);
-    }
-
-    @Override
-    public LatLon eastNorth2latlon(EastNorth p) {
-        return mapXYToLatLon((p.east() - offsetEastMeters)/UTMScaleFactor, (p.north() - offsetNorthMeters)/UTMScaleFactor, UTMCentralMeridianRad);
-    }
-
 }

trunk/src/org/openstreetmap/josm/data/projection/Lambert.java

@@ -6 +6 @@
 import java.awt.GridBagLayout;
 import java.awt.event.ActionListener;
-import java.io.IOException;
 import java.io.InputStream;
 import java.util.Collection;
@@ -17 +16 @@
 import org.openstreetmap.josm.Main;
 import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
 import org.openstreetmap.josm.tools.GBC;
 import org.openstreetmap.josm.tools.ImageProvider;
 
 /**
- * This class provides the two methods <code>latlon2eastNorth</code> and <code>eastNorth2latlon</code>
- * converting the JOSM LatLon coordinates in WGS84 system (GPS) to and from East North values in
- * the projection Lambert conic conform 4 zones using the French geodetic system NTF.
+ * Lambert conic conform 4 zones using the French geodetic system NTF.
  * This newer version uses the grid translation NTF<->RGF93 provided by IGN for a submillimetric accuracy.
  * (RGF93 is the French geodetic system similar to WGS84 but not mathematically equal)
  * @author Pieren
  */
-public class Lambert implements Projection, ProjectionSubPrefs {
+public class Lambert extends AbstractProjection implements ProjectionSubPrefs {
     /**
      * Lambert I, II, III, and IV projection exponents
      */
-    public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
+    private static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
 
     /**
      * Lambert I, II, III, and IV projection constants
      */
-    public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
+    private static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
 
     /**
      * Lambert I, II, III, and IV false east
      */
-    public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 };
+    private static final double x_0s[] = { 600000.0, 600000.0, 600000.0, 234.358 };
 
     /**
      * Lambert I, II, III, and IV false north
      */
-    public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
-
-    /**
-     * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
-     */
-    public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
-
-    /**
-     * precision in iterative schema
-     */
-
-    public static final double epsilon = 1e-11;
+    private static final double y_fs[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
 
     /**
@@ -68 +54 @@
     public static final double cMinLatZone1Radian = Math.toRadians(46.1 * 0.9);// lowest latitude of Zone 4 (South Corsica)
 
-    public static final double zoneLimitsDegree[][] = {
+    public static final double[][] zoneLimitsDegree = {
         {Math.toDegrees(cMaxLatZone1Radian), (53.5 * 0.9)}, // Zone 1  (reference values in grad *0.9)
         {(53.5 * 0.9), (50.5 * 0.9)}, // Zone 2
@@ -86 +72 @@
     public static final int DEFAULT_ZONE = 0;
 
-    private int layoutZone = DEFAULT_ZONE;
+    private int layoutZone;
 
     private static NTV2GridShiftFile ntf_rgf93Grid = null;
@@ -100 +86 @@
                 InputStream is = Main.class.getResourceAsStream("/data/"+gridFileName);
                 if (is == null) {
-                    System.err.println(tr("Warning: failed to open input stream for resource ''/data/{0}''. Cannot load NTF<->RGF93 grid", gridFileName));
-                    return;
+                    throw new RuntimeException(tr("Warning: failed to open input stream for resource ''/data/{0}''. Cannot load NTF<->RGF93 grid", gridFileName));
                 }
                 ntf_rgf93Grid = new NTV2GridShiftFile();
                 ntf_rgf93Grid.loadGridShiftFile(is, false);
-                //System.out.println("NTF<->RGF93 grid loaded.");
             } catch (Exception e) {
-                e.printStackTrace();
+                throw new RuntimeException(e);
             }
         }
-    }
-
-    /**
-     * @param p  WGS84 lat/lon (ellipsoid GRS80) (in degree)
-     * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
-     * @throws IOException
-     */
-    public EastNorth latlon2eastNorth(LatLon p) {
-        // translate ellipsoid GRS80 (WGS83) => Clark
-        LatLon geo = WGS84_to_NTF(p);
-        double lt = Math.toRadians(geo.lat()); // in radian
-        double lg = Math.toRadians(geo.lon());
-
-        // check if longitude and latitude are inside the French Lambert zones
-        if (lt >= cMinLatZone1Radian && lt <= cMaxLatZone1Radian && lg >= cMinLonZonesRadian && lg <= cMaxLonZonesRadian)
-            return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
-        return ConicProjection(lt, lg, Xs[0], Ys[0], c[0], n[0]);
-    }
-
-    public LatLon eastNorth2latlon(EastNorth p) {
-        LatLon geo;
-        geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
-        // translate geodetic system from NTF (ellipsoid Clark) to RGF93/WGS84 (ellipsoid GRS80)
-        return NTF_to_WGS84(geo);
-    }
-
-    @Override public String toString() {
+        updateParameters(DEFAULT_ZONE);
+    }
+    
+    private void updateParameters(int layoutZone) {
+        this.layoutZone = layoutZone;
+        ellps = Ellipsoid.clarke;
+        datum = null; // no datum needed, we have a shift file
+        nadgrids = ntf_rgf93Grid;
+        x_0 = x_0s[layoutZone];
+        lon_0 = 2.0 + 20.0 / 60 + 14.025 / 3600; // 0 grade Paris
+        if (proj == null) {
+            proj = new LambertConformalConic();
+        }
+        ((LambertConformalConic)proj).updateParametersDirect(
+                Ellipsoid.clarke, n[layoutZone], c[layoutZone] / ellps.a, y_fs[layoutZone] / ellps.a);
+    }
+
+    @Override 
+    public String toString() {
         return tr("Lambert 4 Zones (France)");
     }
 
-    public String toCode() {
-        return "EPSG:"+(27561+layoutZone);
+    @Override
+    public Integer getEpsgCode() {
+        return 27561+layoutZone;
     }
 
@@ -149 +126 @@
     }
 
+    @Override
     public String getCacheDirectoryName() {
         return "lambert";
     }
 
-    /**
-     * Initializes from geographic coordinates. Note that reference ellipsoid
-     * used by Lambert is the Clark ellipsoid.
-     *
-     * @param lat latitude in grad
-     * @param lon longitude in grad
-     * @param Xs  false east (coordinate system origin) in meters
-     * @param Ys  false north (coordinate system origin) in meters
-     * @param c   projection constant
-     * @param n   projection exponent
-     * @return EastNorth projected coordinates in meter
-     */
-    private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) {
-        double eslt = Ellipsoid.clarke.e * Math.sin(lat);
-        double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0))
-                * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0));
-        double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0));
-        double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0));
-        return new EastNorth(east, north);
-    }
-
-    /**
-     * Initializes from projected coordinates (conic projection). Note that
-     * reference ellipsoid used by Lambert is Clark
-     *
-     * @param eastNorth projected coordinates pair in meters
-     * @param Xs        false east (coordinate system origin) in meters
-     * @param Ys        false north (coordinate system origin) in meters
-     * @param c         projection constant
-     * @param n         projection exponent
-     * @return LatLon in degree
-     */
-    private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) {
-        double dx = eastNorth.east() - Xs;
-        double dy = Ys - eastNorth.north();
-        double R = Math.sqrt(dx * dx + dy * dy);
-        double gamma = Math.atan(dx / dy);
-        double l = -1.0 / n * Math.log(Math.abs(R / c));
-        l = Math.exp(l);
-        double lon = lg0 + gamma / n;
-        double lat = 2.0 * Math.atan(l) - Math.PI / 2.0;
-        double delta = 1.0;
-        while (delta > epsilon) {
-            double eslt = Ellipsoid.clarke.e * Math.sin(lat);
-            double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI
-            / 2.0;
-            delta = Math.abs(nlt - lat);
-            lat = nlt;
-        }
-        return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon)); // in radian
-    }
-
-    /**
-     * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert
-     * geographic, (ellipsoid Clark)
-     * @throws IOException
-     */
-    private LatLon WGS84_to_NTF(LatLon wgs) {
-        NTV2GridShift gs = new NTV2GridShift(wgs);
-        if (ntf_rgf93Grid != null) {
-            ntf_rgf93Grid.gridShiftReverse(gs);
-            return new LatLon(wgs.lat()+gs.getLatShiftDegrees(), wgs.lon()+gs.getLonShiftPositiveEastDegrees());
-        }
-        return new LatLon(0,0);
-    }
-
-    private LatLon NTF_to_WGS84(LatLon ntf) {
-        NTV2GridShift gs = new NTV2GridShift(ntf);
-        if (ntf_rgf93Grid != null) {
-            ntf_rgf93Grid.gridShiftForward(gs);
-            return new LatLon(ntf.lat()+gs.getLatShiftDegrees(), ntf.lon()+gs.getLonShiftPositiveEastDegrees());
-        }
-        return new LatLon(0,0);
-    }
-
+    @Override
     public Bounds getWorldBoundsLatLon()
     {
@@ -234 +138 @@
                 new LatLon(Math.min(zoneLimitsDegree[layoutZone][0] + cMaxOverlappingZonesDegree, Math.toDegrees(cMaxLatZone1Radian)), Math.toDegrees(cMaxLonZonesRadian)));
         return b;
-    }
-
-    /**
-     * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
-     */
-    public double getDefaultZoomInPPD() {
-        // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
-        return 10.0;
     }
 
@@ -273 +169 @@
     }
 
+    @Override
     public Collection<String> getPreferences(JPanel p) {
         Object prefcb = p.getComponent(2);
@@ -281 +178 @@
     }
 
+    @Override
     public void setPreferences(Collection<String> args) {
-        layoutZone = DEFAULT_ZONE;
+        int layoutZone = DEFAULT_ZONE;
         if (args != null) {
             try {
@@ -295 +193 @@
             } catch(NumberFormatException e) {}
         }
+        updateParameters(layoutZone);
     }
 
@@ -306 +205 @@
     }
 
+    @Override
     public Collection<String> getPreferencesFromCode(String code) {
         if (code.startsWith("EPSG:2756") && code.length() == 10) {

trunk/src/org/openstreetmap/josm/data/projection/LambertCC9Zones.java

@@ -14 +14 @@
 
 import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
+import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
 import org.openstreetmap.josm.tools.GBC;
 import org.openstreetmap.josm.tools.ImageProvider;
 
 /**
- * This class implements the Lambert Conic Conform 9 Zones projection as specified by the IGN
+ * Lambert Conic Conform 9 Zones projection as specified by the IGN
  * in this document http://professionnels.ign.fr/DISPLAY/000/526/700/5267002/transformation.pdf
  * @author Pieren
  *
  */
-public class LambertCC9Zones implements Projection, ProjectionSubPrefs {
-
-    /**
-     * Lambert 9 zones projection exponents
-     */
-    public static final double n[] = { 0.6691500006885269, 0.682018118346418, 0.6946784863203991, 0.7071272481559119,
-        0.7193606118567315, 0.7313748510399917, 0.7431663060711892, 0.7547313851789208, 0.7660665655489937};
-
-    /**
-     * Lambert 9 zones projection constants
-     */
-    public static final double c[] = { 1.215363305807804E7, 1.2050261119223533E7, 1.195716926884592E7, 1.18737533925172E7,
-        1.1799460698022118E7, 1.17337838820243E7, 1.16762559948139E7, 1.1626445901183508E7, 1.1583954251630554E7};
-
-    /**
-     * Lambert 9 zones false east
-     */
-    public static final double Xs = 1700000;
-
-    /**
-     * Lambert 9 zones false north
-     */
-    public static final double Ys[] = { 8293467.503439436, 9049604.665107645, 9814691.693461388, 1.0588107871787189E7,
-        1.1369285637569271E7, 1.2157704903382052E7, 1.2952888086405803E7, 1.3754395745267643E7, 1.4561822739114787E7};
-
-    /**
-     * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
-     */
-    public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
-
-    /**
-     * precision in iterative schema
-     */
-
-    public static final double epsilon = 1e-12;
+public class LambertCC9Zones extends AbstractProjection implements ProjectionSubPrefs {
 
     /**
@@ -67 +34 @@
 
     public static final double cMinLatZonesDegree = 41.0;
-    public static final double cMinLatZonesRadian = Math.toRadians(cMinLatZonesDegree);
-
-    public static final double cMinLonZonesRadian = Math.toRadians(-5.0);
-
-    public static final double cMaxLonZonesRadian = Math.toRadians(10.2);
-
-    public static final double lambda0 = Math.toRadians(3);
-    public static final double e = Ellipsoid.GRS80.e; // but in doc=0.08181919112
-    public static final double e2 =Ellipsoid.GRS80.e2;
-    public static final double a = Ellipsoid.GRS80.a;
 
     public static final double cMaxOverlappingZones = 1.5;
@@ -82 +39 @@
     public static final int DEFAULT_ZONE = 0;
 
-    private static int layoutZone = DEFAULT_ZONE;
+    private int layoutZone = DEFAULT_ZONE;
 
-    private double L(double phi, double e) {
-        double sinphi = Math.sin(phi);
-        return (0.5*Math.log((1+sinphi)/(1-sinphi))) - e/2*Math.log((1+e*sinphi)/(1-e*sinphi));
+    public LambertCC9Zones() {
+        this(DEFAULT_ZONE);
     }
 
-    /**
-     * @param p  WGS84 lat/lon (ellipsoid GRS80) (in degree)
-     * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
-     */
-    public EastNorth latlon2eastNorth(LatLon p) {
-        double lt = Math.toRadians(p.lat());
-        double lg = Math.toRadians(p.lon());
-        if (lt >= cMinLatZonesRadian && lt <= cMaxLatZonesRadian && lg >= cMinLonZonesRadian && lg <= cMaxLonZonesRadian)
-            return ConicProjection(lt, lg, layoutZone);
-        return ConicProjection(lt, lg, 0);
+    public LambertCC9Zones(int layoutZone) {
+        updateParameters(layoutZone);
     }
 
-    /**
-     *
-     * @param lat latitude in grad
-     * @param lon longitude in grad
-     * @param nz Lambert CC zone number (from 1 to 9) - 1 !
-     * @return EastNorth projected coordinates in meter
-     */
-    private EastNorth ConicProjection(double lat, double lon, int nz) {
-        double R = c[nz]*Math.exp(-n[nz]*L(lat,e));
-        double gamma = n[nz]*(lon-lambda0);
-        double X = Xs + R*Math.sin(gamma);
-        double Y = Ys[nz] + -R*Math.cos(gamma);
-        return new EastNorth(X, Y);
+    public void updateParameters(int layoutZone) {
+        ellps = Ellipsoid.GRS80;
+        datum = GRS80Datum.INSTANCE;
+        this.layoutZone = layoutZone;
+        x_0 = 1700000;
+        y_0 = (layoutZone+1) * 1000000 + 200000;
+        lon_0 = 3;
+        double lat_0 = 42 + layoutZone;
+        double lat_1 = 41.25 + layoutZone;
+        double lat_2 = 42.75 + layoutZone;
+        if (proj == null) {
+            proj = new LambertConformalConic();
+        }
+        ((LambertConformalConic)proj).updateParameters2SP(ellps, lat_0, lat_1, lat_2);
     }
 
-    public LatLon eastNorth2latlon(EastNorth p) {
-        return Geographic(p, layoutZone);
-    }
-
-    private LatLon Geographic(EastNorth ea, int nz) {
-        double R = Math.sqrt(Math.pow(ea.getX()-Xs,2)+Math.pow(ea.getY()-Ys[nz], 2));
-        double gamma = Math.atan((ea.getX()-Xs)/(Ys[nz]-ea.getY()));
-        double lon = lambda0+gamma/n[nz];
-        double latIso = (-1/n[nz])*Math.log(Math.abs(R/c[nz]));
-        double lat = Ellipsoid.GRS80.latitude(latIso, e, epsilon);
-        return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon));
-    }
-
-    @Override public String toString() {
+    @Override 
+    public String toString() {
         return tr("Lambert CC9 Zone (France)");
     }
@@ -140 +77 @@
     }
 
-    public String toCode() {
-        return "EPSG:"+(3942+layoutZone); //CC42 is EPSG:3942 (up to EPSG:3950 for CC50)
+    @Override
+    public Integer getEpsgCode() {
+        return 3942+layoutZone; //CC42 is EPSG:3942 (up to EPSG:3950 for CC50)
     }
 
@@ -149 +87 @@
     }
 
+    @Override
     public String getCacheDirectoryName() {
         return "lambert";
     }
 
-    /**
-     * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
-     */
-    public double getDefaultZoomInPPD() {
-        // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
-        return 10.0;
-    }
-
+    @Override
     public Bounds getWorldBoundsLatLon()
     {
@@ -203 +135 @@
     }
 
+    @Override
     public Collection<String> getPreferences(JPanel p) {
         Object prefcb = p.getComponent(2);
         if(!(prefcb instanceof JComboBox))
             return null;
-        layoutZone = ((JComboBox)prefcb).getSelectedIndex();
+        int layoutZone = ((JComboBox)prefcb).getSelectedIndex();
         return Collections.singleton(Integer.toString(layoutZone+1));
     }
 
+    @Override
     public void setPreferences(Collection<String> args)
     {
-        layoutZone = DEFAULT_ZONE;
+        int layoutZone = DEFAULT_ZONE;
         if (args != null) {
             try {
@@ -226 +160 @@
             } catch(NumberFormatException e) {}
         }
+        updateParameters(layoutZone);
     }
 
@@ -237 +172 @@
     }
 
+    @Override
     public Collection<String> getPreferencesFromCode(String code)
     {
@@ -246 +182 @@
                 if(zoneval >= 0 && zoneval <= 8)
                     return Collections.singleton(String.valueOf(zoneval+1));
-            } catch(NumberFormatException e) {}
+            } catch(NumberFormatException ex) {}
         }
         return null;

trunk/src/org/openstreetmap/josm/data/projection/LambertEST.java

@@ -1 @@
-//License: GPL. For details, see LICENSE file.
-//Thanks to Johan Montagnat and its geoconv java converter application
-//(http://www.i3s.unice.fr/~johan/gps/ , published under GPL license)
-//from which some code and constants have been reused here.
+// License: GPL. For details, see LICENSE file.
 package org.openstreetmap.josm.data.projection;
 
@@ -8 +5 @@
 
 import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
+import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
 
-public class LambertEST implements Projection {
+/**
+ * Estonian Coordinate System of 1997.
+ * 
+ * Thanks to Johan Montagnat and its geoconv java converter application
+ * (http://www.i3s.unice.fr/~johan/gps/ , published under GPL license)
+ * from which some code and constants have been reused here.
+ */
+public class LambertEST extends AbstractProjection {
 
-    public static final double ef = 500000; //Easting at false origin            = 5000000 m
-    public static final double nf = 6375000; //Northing at false origin           = 6375000 m
-    public static final double lat1 = Math.toRadians(59 + 1.0/3.0); //Latitude of 1st standard parallel  = 59o20`0`` N
-    public static final double lat2 = Math.toRadians( 58);//Latitude of 2nd standard parallel  = 58o0`0`` N
-    public static final double latf = Math.toRadians(57.517553930555555555555555555556);//'Latitude of false origin = 57o31`3,19415`` N
-    public static final double lonf = Math.toRadians( 24.0);
-    public static final double a = 6378137;
-    public static final double ee = 0.081819191042815792;
-    public static final double m1 = Math.cos(lat1) / (Math.sqrt(1 - ee *ee * Math.pow(Math.sin(lat1), 2)));
-    public static final double m2 = Math.cos(lat2) / (Math.sqrt(1 - ee *ee * Math.pow(Math.sin(lat2), 2)));
-    public static final double t1 = Math.tan(Math.PI / 4.0 - lat1 / 2.0)
-    / Math.pow(( (1.0 - ee * Math.sin(lat1)) / (1.0 + ee * Math.sin(lat1))) ,(ee / 2.0));
-    public static final double t2 = Math.tan(Math.PI / 4.0 - lat2 / 2.0)
-    / Math.pow(( (1.0 - ee * Math.sin(lat2)) / (1.0 + ee * Math.sin(lat2))) ,(ee / 2.0));
-    public static final double tf = Math.tan(Math.PI / 4.0 - latf / 2.0)
-    / Math.pow(( (1.0 - ee * Math.sin(latf)) / (1.0 + ee * Math.sin(latf))) ,(ee / 2.0));
-    public static final double n  = (Math.log(m1) - Math.log(m2))
-    / (Math.log(t1) - Math.log(t2));
-    public static final double f  = m1 / (n * Math.pow(t1, n));
-    public static final double rf  = a * f * Math.pow(tf, n);
-
-    /**
-     * precision in iterative schema
-     */
-    public static final double epsilon = 1e-11;
-
-    /**
-     * @param p  WGS84 lat/lon (ellipsoid GRS80) (in degree)
-     * @return eastnorth projection in Lambert Zone (ellipsoid GRS80)
-     */
-    public EastNorth latlon2eastNorth(LatLon p)
-    {
-
-        double t = Math.tan(Math.PI / 4.0 - Math.toRadians(p.lat()) / 2.0)
-        / Math.pow(( (1.0 - ee * Math.sin(Math.toRadians(p.lat()))) / (1.0
-                + ee * Math.sin(Math.toRadians(p.lat())))) ,(ee / 2.0));
-        double r = a * f * Math.pow(t, n);
-        double theta = n * (Math.toRadians(p.lon()) - lonf);
-
-        double x = ef + r * Math.sin(theta);     //587446.7
-        double y = nf + rf - r * Math.cos(theta); //6485401.6
-
-        return new EastNorth(x,y);
-    }
-
-    public static double iterateAngle(double e, double t)
-    {
-        double a1 = 0.0;
-        double a2 = 3.1415926535897931;
-        double a = 1.5707963267948966;
-        double b = 1.5707963267948966 - (2.0 * Math.atan(t * Math.pow((1.0
-                - (e * Math.sin(a))) / (1.0 + (e * Math.sin(a))), e / 2.0)));
-        while (Math.abs(a-b) > epsilon)
-        {
-            a = a1 + ((a2 - a1) / 2.0);
-            b = 1.5707963267948966 - (2.0 * Math.atan(t * Math.pow((1.0
-                    - (e * Math.sin(a))) / (1.0 + (e * Math.sin(a))), e / 2.0)));
-            if (a1 == a2)
-                return 0.0;
-            if (b > a) {
-                a1 = a;
-            } else {
-                a2 = a;
-            }
-        }
-        return b;
-    }
-
-    public LatLon eastNorth2latlon(EastNorth p)
-    {
-        double r = Math.sqrt(Math.pow((p.getX() - ef), 2.0) + Math.pow((rf
-                - p.getY() + nf), 2.0) ) * Math.signum(n);
-        double T = Math.pow((r / (a * f)), (1.0/ n)) ;
-        double theta = Math.atan((p.getX() - ef) / (rf - p.getY() + nf));
-        double y = (theta / n + lonf) ;
-        double x = iterateAngle(ee, T);
-        return new LatLon(Math.toDegrees(x),Math.toDegrees(y));
+    public LambertEST() {
+        ellps = Ellipsoid.GRS80;
+        datum = GRS80Datum.INSTANCE;
+        lon_0 = 24;
+        double lat_0 = 57.517553930555555555555555555556;
+        double lat_1 = 59 + 1.0/3.0;
+        double lat_2 = 58;
+        x_0 = 500000;
+        y_0 = 6375000;
+        proj = new LambertConformalConic();
+        ((LambertConformalConic) proj).updateParameters2SP(ellps, lat_0, lat_1, lat_2);
     }
 
@@ -97 +36 @@
     }
 
-    public String toCode() {
-        return "EPSG:3301";
+    @Override
+    public Integer getEpsgCode() {
+        return 3301;
     }
 
@@ -106 +46 @@
     }
 
+    @Override
     public String getCacheDirectoryName() {
         return "lambertest";
     }
 
+    @Override
     public Bounds getWorldBoundsLatLon()
     {
@@ -117 +59 @@
     }
 
-    public double getDefaultZoomInPPD() {
-        return 10;
-    }
 }

trunk/src/org/openstreetmap/josm/data/projection/Mercator.java

@@ -5 +5 @@
 
 import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.WGS84Datum;
 
 /**
- * Implement Mercator Projection code, coded after documentation
- * from wikipedia.
+ * Mercator Projection
  *
  * The center of the mercator projection is always the 0 grad
@@ -20 +19 @@
  * @author imi
  */
-public class Mercator implements Projection {
+public class Mercator extends AbstractProjection {
 
-    final double radius = 6378137.0;
-
-    public EastNorth latlon2eastNorth(LatLon p) {
-        return new EastNorth(
-                p.lon()*Math.PI/180*radius,
-                Math.log(Math.tan(Math.PI/4+p.lat()*Math.PI/360))*radius);
+    public Mercator() {
+        ellps = Ellipsoid.WGS84;
+        datum = WGS84Datum.INSTANCE;
+        proj = new org.openstreetmap.josm.data.projection.proj.Mercator();
     }
 
-    public LatLon eastNorth2latlon(EastNorth p) {
-        return new LatLon(
-                Math.atan(Math.sinh(p.north()/radius))*180/Math.PI,
-                p.east()/radius*180/Math.PI);
-    }
-
-    @Override public String toString() {
+    @Override 
+    public String toString() {
         return tr("Mercator");
     }
 
-    public String toCode() {
-        return "EPSG:3857"; /* initially they used 3785 but that has been superseded, see http://www.epsg-registry.org/ */
+    @Override
+    public Integer getEpsgCode() {
+        /* initially they used 3785 but that has been superseded, 
+         * see http://www.epsg-registry.org/ */
+        return 3857;
     }
 
@@ -49 +44 @@
     }
 
+    @Override
     public String getCacheDirectoryName() {
         return "mercator";
     }
 
+    @Override
     public Bounds getWorldBoundsLatLon()
     {
@@ -60 +57 @@
     }
 
-    public double getDefaultZoomInPPD() {
-        // This will set the scale bar to about 100 km
-        return 1000.0;/*0.000158*/
-    }
 }

trunk/src/org/openstreetmap/josm/data/projection/Puwg.java

@@ -7 +7 @@
 import java.awt.GridBagLayout;
 import java.awt.event.ActionListener;
-import java.text.DecimalFormat;
 import java.util.Collection;
 import java.util.Collections;
@@ -16 +15 @@
 
 import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
 import org.openstreetmap.josm.tools.GBC;
 
@@ -26 +25 @@
  * @author steelman
  */
-public class Puwg extends UTM implements Projection,ProjectionSubPrefs {
+public class Puwg extends AbstractProjection implements ProjectionSubPrefs {
+    
     public static final int DEFAULT_ZONE = 0;
-    private int zone = DEFAULT_ZONE;
-
-    static PuwgData[] Zones = new PuwgData[]{
+    
+    private int zone;
+
+    static PuwgData[] Zones = new PuwgData[] {
         new Epsg2180(),
         new Epsg2176(),
@@ -38 +39 @@
     };
 
-    private static DecimalFormat decFormatter = new DecimalFormat("###0");
-
-    @Override
-    public EastNorth latlon2eastNorth(LatLon p) {
+    public Puwg() {
+        this(DEFAULT_ZONE);
+    }
+
+    public Puwg(int zone) {
+        ellps = Ellipsoid.GRS80;
+        proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
+        datum = GRS80Datum.INSTANCE;
+        updateParameters(zone);
+    }
+    
+    public void updateParameters(int zone) {
+        this.zone = zone;
         PuwgData z = Zones[zone];
-        double easting = z.getPuwgFalseEasting();
-        double northing = z.getPuwgFalseNorthing();
-        double scale = z.getPuwgScaleFactor();
-        double center = z.getPuwgCentralMeridian(); /* in radians */
-        EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), center);
-        return new EastNorth(a.east() * scale + easting, a.north() * scale + northing);
-    }
-
-    @Override
-    public LatLon eastNorth2latlon(EastNorth p) {
-        PuwgData z = Zones[zone];
-        double easting = z.getPuwgFalseEasting();
-        double northing = z.getPuwgFalseNorthing();
-        double scale = z.getPuwgScaleFactor();
-        double center = z.getPuwgCentralMeridian(); /* in radians */
-        return mapXYToLatLon((p.east() - easting)/scale, (p.north() - northing)/scale, center);
-    }
-
-    @Override public String toString() {
+        x_0 = z.getPuwgFalseEasting();
+        y_0 = z.getPuwgFalseNorthing();
+        lon_0 = z.getPuwgCentralMeridianDeg();
+        k_0 = z.getPuwgScaleFactor();
+    }
+
+    @Override 
+    public String toString() {
         return tr("PUWG (Poland)");
     }
 
     @Override
-    public String toCode() {
-        return Zones[zone].toCode();
+    public Integer getEpsgCode() {
+        return Zones[zone].getEpsgCode();
     }
 
@@ -83 +82 @@
     public Bounds getWorldBoundsLatLon() {
         return Zones[zone].getWorldBoundsLatLon();
-    }
-
-    @Override
-    public double getDefaultZoomInPPD() {
-        // This will set the scale bar to about 100 km
-        return 0.009;
-    }
-
-    public String eastToString(EastNorth p) {
-        return decFormatter.format(p.east());
-    }
-
-    public String northToString(EastNorth p) {
-        return decFormatter.format(p.north());
     }
 
@@ -135 +120 @@
 
     @Override
-    public Collection<String> getPreferencesFromCode(String code)
-    {
-        for (PuwgData p : Puwg.Zones)
-        {
-            if(code.equals(p.toCode()))
+    public Collection<String> getPreferencesFromCode(String code) {
+        for (PuwgData p : Puwg.Zones) {
+            if (code.equals(p.toCode()))
                 return Collections.singleton(code);
         }
@@ -146 +129 @@
 
     @Override
-    public void setPreferences(Collection<String> args)
-    {
-        zone = DEFAULT_ZONE;
-        if(args != null)
-        {
+    public void setPreferences(Collection<String> args) {
+        int z = DEFAULT_ZONE;
+        if (args != null) {
             try {
-                for(String s : args)
-                {
-                    for (int i=0; i < Puwg.Zones.length; ++i)
-                        if(s.equals(Zones[i].toCode())) {
-                            zone = i;
+                for (String s : args) {
+                    for (int i=0; i < Zones.length; ++i)
+                        if (s.equals(Zones[i].toCode())) {
+                            z = i;
+                            break;
                         }
                     break;
@@ -162 +143 @@
             } catch (NullPointerException e) {}
         }
+        updateParameters(z);
     }
 }
@@ -173 +155 @@
 
     // Projection methods
+    Integer getEpsgCode();
     String toCode();
     String getCacheDirectoryName();
@@ -189 +172 @@
     }
 
+    @Override
+    public Integer getEpsgCode() {
+        return 2180;
+    }
+    
+    @Override
     public String toCode() {
-        return "EPSG:2180";
-    }
-
+        return "EPSG:" + getEpsgCode();
+    }
+
+    @Override
     public String getCacheDirectoryName() {
         return "epsg2180";
     }
 
+    @Override
     public Bounds getWorldBoundsLatLon()
     {
@@ -204 +195 @@
     }
 
-    public double getPuwgCentralMeridianDeg() { return Epsg2180CentralMeridian; }
-    public double getPuwgCentralMeridian() { return Math.toRadians(Epsg2180CentralMeridian); }
-    public double getPuwgFalseEasting() { return Epsg2180FalseEasting; }
-    public double getPuwgFalseNorthing() { return Epsg2180FalseNorthing; }
-    public double getPuwgScaleFactor() { return Epsg2180ScaleFactor; }
+    @Override public double getPuwgCentralMeridianDeg() { return Epsg2180CentralMeridian; }
+    @Override public double getPuwgCentralMeridian() { return Math.toRadians(Epsg2180CentralMeridian); }
+    @Override public double getPuwgFalseEasting() { return Epsg2180FalseEasting; }
+    @Override public double getPuwgFalseNorthing() { return Epsg2180FalseNorthing; }
+    @Override public double getPuwgScaleFactor() { return Epsg2180ScaleFactor; }
 }
 
@@ -217 +208 @@
     private static final double PuwgScaleFactor = 0.999923;
     //final private double[] Puwg2000CentralMeridian = {15.0, 18.0, 21.0, 24.0};
-    final private String[] Puwg2000Code = { "EPSG:2176",  "EPSG:2177", "EPSG:2178", "EPSG:2179"};
-    final private String[] Puwg2000CDName = { "epsg2176",  "epsg2177", "epsg2178", "epsg2179"};
+    final private Integer[] Puwg2000Code = { 2176,  2177, 2178, 2179 };
+    final private String[] Puwg2000CDName = { "epsg2176",  "epsg2177", "epsg2178", "epsg2179" };
 
     @Override public String toString() {
@@ -224 +215 @@
     }
 
+    @Override
+    public Integer getEpsgCode() {
+        return Puwg2000Code[getZoneIndex()];
+    }
+
+    @Override
     public String toCode() {
-        return Puwg2000Code[getZoneIndex()];
-    }
-
+        return "EPSG:" + getEpsgCode();
+    }
+
+    @Override
     public String getCacheDirectoryName() {
         return Puwg2000CDName[getZoneIndex()];
     }
 
+    @Override
     public Bounds getWorldBoundsLatLon()
     {
@@ -239 +238 @@
     }
 
-    public double getPuwgCentralMeridianDeg() { return getZone() * 3.0; }
-    public double getPuwgCentralMeridian() { return Math.toRadians(getZone() * 3.0); }
-    public double getPuwgFalseNorthing() { return PuwgFalseNorthing;}
-    public double getPuwgFalseEasting() { return 1e6 * getZone() + PuwgFalseEasting; }
-    public double getPuwgScaleFactor() { return PuwgScaleFactor; }
+    @Override public double getPuwgCentralMeridianDeg() { return getZone() * 3.0; }
+    @Override public double getPuwgCentralMeridian() { return Math.toRadians(getZone() * 3.0); }
+    @Override public double getPuwgFalseNorthing() { return PuwgFalseNorthing;}
+    @Override public double getPuwgFalseEasting() { return 1e6 * getZone() + PuwgFalseEasting; }
+    @Override public double getPuwgScaleFactor() { return PuwgScaleFactor; }
     public abstract int getZone();
 

trunk/src/org/openstreetmap/josm/data/projection/SwissGrid.java

@@ -12 +12 @@
 
 import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.ThreeParameterDatum;
+import org.openstreetmap.josm.data.projection.proj.SwissObliqueMercator;
 import org.openstreetmap.josm.gui.widgets.HtmlPanel;
 import org.openstreetmap.josm.tools.GBC;
 
 /**
- * Projection for the SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid.
+ * SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid.
+ * 
+ * Actually, what we have here, is CH1903+ (EPSG:2056), but without
+ * the additional false easting of 2000km and false northing 1000 km.
  *
- * Calculations are based on formula from
- * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf
- *
- * August 2010 update to this formula (rigorous formulas)
- * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf
+ * To get to CH1903, a shift file is required. So currently, there are errors
+ * up to 1.6m (depending on the location).
+ * 
+ * This projection does not have any parameters, it only implements
+ * ProjectionSubPrefs to show a warning that the grid file correction is not done.
  */
-public class SwissGrid implements Projection, ProjectionSubPrefs {
+public class SwissGrid extends AbstractProjection implements ProjectionSubPrefs {
 
-    private static final double dX = 674.374;
-    private static final double dY = 15.056;
-    private static final double dZ = 405.346;
-
-    private static final double phi0 = Math.toRadians(46.0 + 57.0 / 60 + 8.66 / 3600);
-    private static final double lambda0 = Math.toRadians(7.0 + 26.0 / 60 + 22.50 / 3600);
-    private static final double R = Ellipsoid.Bessel1841.a * Math.sqrt(1 - Ellipsoid.Bessel1841.e2) / (1 - (Ellipsoid.Bessel1841.e2 * Math.pow(Math.sin(phi0), 2)));
-    private static final double alpha = Math.sqrt(1 + (Ellipsoid.Bessel1841.eb2 * Math.pow(Math.cos(phi0), 4)));
-    private static final double b0 = Math.asin(Math.sin(phi0) / alpha);
-    private static final double K = Math.log(Math.tan(Math.PI / 4 + b0 / 2)) - alpha
-    * Math.log(Math.tan(Math.PI / 4 + phi0 / 2)) + alpha * Ellipsoid.Bessel1841.e / 2
-    * Math.log((1 + Ellipsoid.Bessel1841.e * Math.sin(phi0)) / (1 - Ellipsoid.Bessel1841.e * Math.sin(phi0)));
-
-    private static final double xTrans = 200000;
-    private static final double yTrans = 600000;
-
-    private static final double DELTA_PHI = 1e-11;
-
-    private LatLon correctEllipoideGSR80toBressel1841(LatLon coord) {
-        double[] XYZ = Ellipsoid.WGS84.latLon2Cart(coord);
-        XYZ[0] -= dX;
-        XYZ[1] -= dY;
-        XYZ[2] -= dZ;
-        return Ellipsoid.Bessel1841.cart2LatLon(XYZ);
-    }
-
-    private LatLon correctEllipoideBressel1841toGRS80(LatLon coord) {
-        double[] XYZ = Ellipsoid.Bessel1841.latLon2Cart(coord);
-        XYZ[0] += dX;
-        XYZ[1] += dY;
-        XYZ[2] += dZ;
-        return Ellipsoid.WGS84.cart2LatLon(XYZ);
-    }
-
-    /**
-     * @param wgs  WGS84 lat/lon (ellipsoid GRS80) (in degree)
-     * @return eastnorth projection in Swiss National Grid (ellipsoid Bessel)
-     */
-    @Override
-    public EastNorth latlon2eastNorth(LatLon wgs) {
-        LatLon coord = correctEllipoideGSR80toBressel1841(wgs);
-        double phi = Math.toRadians(coord.lat());
-        double lambda = Math.toRadians(coord.lon());
-
-        double S = alpha * Math.log(Math.tan(Math.PI / 4 + phi / 2)) - alpha * Ellipsoid.Bessel1841.e / 2
-        * Math.log((1 + Ellipsoid.Bessel1841.e * Math.sin(phi)) / (1 - Ellipsoid.Bessel1841.e * Math.sin(phi))) + K;
-        double b = 2 * (Math.atan(Math.exp(S)) - Math.PI / 4);
-        double l = alpha * (lambda - lambda0);
-
-        double lb = Math.atan2(Math.sin(l), Math.sin(b0) * Math.tan(b) + Math.cos(b0) * Math.cos(l));
-        double bb = Math.asin(Math.cos(b0) * Math.sin(b) - Math.sin(b0) * Math.cos(b) * Math.cos(l));
-
-        double y = R * lb;
-        double x = R / 2 * Math.log((1 + Math.sin(bb)) / (1 - Math.sin(bb)));
-
-        return new EastNorth(y + yTrans, x + xTrans);
-    }
-
-    /**
-     * @param xy SwissGrid east/north (in meters)
-     * @return LatLon WGS84 (in degree)
-     */
-    @Override
-    public LatLon eastNorth2latlon(EastNorth xy) {
-        double x = xy.north() - xTrans;
-        double y = xy.east() - yTrans;
-
-        double lb = y / R;
-        double bb = 2 * (Math.atan(Math.exp(x / R)) - Math.PI / 4);
-
-        double b = Math.asin(Math.cos(b0) * Math.sin(bb) + Math.sin(b0) * Math.cos(bb) * Math.cos(lb));
-        double l = Math.atan2(Math.sin(lb), Math.cos(b0) * Math.cos(lb) - Math.sin(b0) * Math.tan(bb));
-
-        double lambda = lambda0 + l / alpha;
-        double phi = b;
-        double S = 0;
-
-        double prevPhi = -1000;
-        int iteration = 0;
-        // iteration to finds S and phi
-        while (Math.abs(phi - prevPhi) > DELTA_PHI) {
-            if (++iteration > 30)
-                throw new RuntimeException("Two many iterations");
-            prevPhi = phi;
-            S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + Ellipsoid.Bessel1841.e
-            * Math.log(Math.tan(Math.PI / 4 + Math.asin(Ellipsoid.Bessel1841.e * Math.sin(phi)) / 2));
-            phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2;
-        }
-
-        LatLon coord = correctEllipoideBressel1841toGRS80(new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda)));
-        return coord;
+    public SwissGrid() {
+        ellps = Ellipsoid.Bessel1841;
+        datum = new ThreeParameterDatum("SwissGrid Datum", null, ellps, 674.374, 15.056, 405.346);
+        x_0 = 600000;
+        y_0 = 200000;
+        lon_0 = 7.0 + 26.0/60 + 22.50/3600;
+        double lat_0 = 46.0 + 57.0/60 + 8.66/3600;
+        proj = new SwissObliqueMercator(ellps, lat_0);
     }
 
@@ -127 +48 @@
 
     @Override
-    public String toCode() {
-        return "EPSG:21781";
+    public Integer getEpsgCode() {
+        return 21781;
     }
 
@@ -144 +65 @@
     public Bounds getWorldBoundsLatLon() {
         return new Bounds(new LatLon(45.7, 5.7), new LatLon(47.9, 10.6));
-    }
-
-    @Override
-    public double getDefaultZoomInPPD() {
-        // This will set the scale bar to about 100 m
-        return 1.01;
     }
 

trunk/src/org/openstreetmap/josm/data/projection/TransverseMercatorLV.java

@@ -6 +6 @@
 import org.openstreetmap.josm.data.Bounds;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
 
 /**
@@ -13 +14 @@
  * @author Viesturs Zarins
  */
-public class TransverseMercatorLV extends TransverseMercator {
+public class TransverseMercatorLV extends AbstractProjection {
 
-    public TransverseMercatorLV()
-    {
-        setProjectionParameters(24, 500000, -6000000);
+    public TransverseMercatorLV() {
+        ellps = Ellipsoid.GRS80;
+        proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
+        datum = GRS80Datum.INSTANCE;
+        lon_0 = 24;
+        x_0 = 500000;
+        y_0 = -6000000;
+        k_0 = 0.9996;
     }
-
-    @Override public String toString() {
+    
+    @Override 
+    public String toString() {
         return tr("LKS-92 (Latvia TM)");
     }
 
-    private int epsgCode() {
+    @Override
+    public Integer getEpsgCode() {
         return 3059;
-    }
-
-    @Override
-    public String toCode() {
-        return "EPSG:"+ epsgCode();
     }
 
@@ -38 +41 @@
     }
 
+    @Override
     public String getCacheDirectoryName() {
-        return "epsg"+ epsgCode();
+        return "epsg"+ getEpsgCode();
     }
 

trunk/src/org/openstreetmap/josm/data/projection/UTM.java

@@ -19 +19 @@
 import org.openstreetmap.josm.data.Bounds;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
 import org.openstreetmap.josm.tools.GBC;
 
@@ -25 +26 @@
  * @author Dirk Stöcker
  * code based on JavaScript from Chuck Taylor
+ *
  */
-public class UTM extends TransverseMercator implements ProjectionSubPrefs {
+public class UTM extends AbstractProjection implements ProjectionSubPrefs {
 
     private static final int DEFAULT_ZONE = 30;
-    private int zone = DEFAULT_ZONE;
-
-    public enum Hemisphere {North, South}
+    private int zone;
+
+    public enum Hemisphere { North, South }
     private static final Hemisphere DEFAULT_HEMISPHERE = Hemisphere.North;
-    private Hemisphere hemisphere = DEFAULT_HEMISPHERE;
-
-    private boolean offset = false;
-
-    public UTM()
-    {
-        updateParameters();
+    private Hemisphere hemisphere;
+
+    /**
+     * Applies an additional false easting of 3000000 m if true.
+     */
+    private boolean offset;
+
+    public UTM() {
+        this(DEFAULT_ZONE, DEFAULT_HEMISPHERE, false);
+    }
+
+    public UTM(int zone, Hemisphere hemisphere, boolean offset) {
+        ellps = Ellipsoid.GRS80;
+        proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
+        datum = GRS80Datum.INSTANCE;
+        updateParameters(zone, hemisphere, offset);
+    }
+
+    public void updateParameters(int zone, Hemisphere hemisphere, boolean offset) {
+        this.zone = zone;
+        this.hemisphere = hemisphere;
+        this.offset = offset;
+        x_0 = 500000 + (offset ? 3000000 : 0);
+        y_0 = hemisphere == Hemisphere.North ? 0 : 10000000;
+        lon_0 = getUtmCentralMeridianDeg(zone);
+        k_0 = 0.9996;
     }
 
@@ -56 +77 @@
      *
      */
-    private double UTMCentralMeridianDeg(int zone)
+    private double getUtmCentralMeridianDeg(int zone)
     {
         return -183.0 + (zone * 6.0);
     }
 
-    @Override public String toString() {
+    public int getzone() {
+        return zone;
+    }
+
+    @Override
+    public String toString() {
         return tr("UTM");
     }
 
-    private void updateParameters() {
-        setProjectionParameters(this.UTMCentralMeridianDeg(getzone()), getEastOffset(), getNorthOffset());
-    }
-
-    public int getzone()
-    {
-        return zone;
-    }
-
-    private double getEastOffset() {
-        return 500000 + (offset?3000000:0);
-    }
-
-    private double getNorthOffset() {
-        if (hemisphere == Hemisphere.North)
-            return 0;
-        else
-            return 10000000;
-    }
-
-    private int epsgCode() {
+    @Override
+    public Integer getEpsgCode() {
         return ((offset?325800:32600) + getzone() + (hemisphere == Hemisphere.South?100:0));
-    }
-
-    public String toCode() {
-        return "EPSG:"+ epsgCode();
     }
 
@@ -98 +101 @@
     }
 
+    @Override
     public String getCacheDirectoryName() {
-        return "epsg"+ epsgCode();
-    }
-
+        return "epsg"+ getEpsgCode();
+    }
+
+    @Override
     public Bounds getWorldBoundsLatLon()
     {
         if (hemisphere == Hemisphere.North)
             return new Bounds(
-                    new LatLon(-5.0, UTMCentralMeridianDeg(getzone())-5.0),
-                    new LatLon(85.0, UTMCentralMeridianDeg(getzone())+5.0));
+                    new LatLon(-5.0, getUtmCentralMeridianDeg(getzone())-5.0),
+                    new LatLon(85.0, getUtmCentralMeridianDeg(getzone())+5.0));
         else
             return new Bounds(
-                    new LatLon(-85.0, UTMCentralMeridianDeg(getzone())-5.0),
-                    new LatLon(5.0, UTMCentralMeridianDeg(getzone())+5.0));
+                    new LatLon(-85.0, getUtmCentralMeridianDeg(getzone())-5.0),
+                    new LatLon(5.0, getUtmCentralMeridianDeg(getzone())+5.0));
     }
 
@@ -173 +178 @@
     }
 
+    @Override
     public Collection<String> getPreferences(JPanel p) {
         int zone = DEFAULT_ZONE;
@@ -199 +205 @@
     }
 
-    public void setPreferences(Collection<String> args)
-    {
-        zone = DEFAULT_ZONE;
-        hemisphere = DEFAULT_HEMISPHERE;
-        offset = false;
+    @Override
+    public void setPreferences(Collection<String> args) {
+        int zone = DEFAULT_ZONE;
+        Hemisphere hemisphere = DEFAULT_HEMISPHERE;
+        boolean offset = false;
 
         if(args != null)
@@ -223 +229 @@
             }
         }
-        updateParameters();
-    }
-
+        updateParameters(zone, hemisphere, offset);
+    }
+
+    @Override
     public String[] allCodes() {
         ArrayList<String> projections = new ArrayList<String>(60*4);
@@ -236 +243 @@
         }
         return projections.toArray(new String[0]);
-
-    }
-
-    public Collection<String> getPreferencesFromCode(String code)
-    {
+    }
+
+    @Override
+    public Collection<String> getPreferencesFromCode(String code) {
+
         boolean offset = code.startsWith("EPSG:3258") || code.startsWith("EPSG:3259");
 

trunk/src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java

@@ -2 +2 @@
 package org.openstreetmap.josm.data.projection;
 
-/**
- * This class implements all projections for French departements in the Caribbean Sea and
- * Indian Ocean using the UTM transvers Mercator projection and specific geodesic settings (7 parameters transformation algorithm).
- */
 import static org.openstreetmap.josm.tools.I18n.tr;
 
@@ -18 +14 @@
 
 import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
 import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.Datum;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
+import org.openstreetmap.josm.data.projection.datum.SevenParameterDatum;
+import org.openstreetmap.josm.data.projection.datum.ThreeParameterDatum;
 import org.openstreetmap.josm.tools.GBC;
 
-public class UTM_France_DOM implements Projection, ProjectionSubPrefs {
+/**
+ * This class implements all projections for French departements in the Caribbean Sea and
+ * Indian Ocean using the UTM transvers Mercator projection and specific geodesic settings (7 parameters transformation algorithm).
+ * 
+ */
+public class UTM_France_DOM extends AbstractProjection implements ProjectionSubPrefs {
 
-    private static String FortMarigotName = tr("Guadeloupe Fort-Marigot 1949");
-    private static String SainteAnneName = tr("Guadeloupe Ste-Anne 1948");
-    private static String MartiniqueName = tr("Martinique Fort Desaix 1952");
-    private static String Reunion92Name = tr("Reunion RGR92");
-    private static String Guyane92Name = tr("Guyane RGFG95");
-    public static String[] utmGeodesicsNames = { FortMarigotName, SainteAnneName, MartiniqueName, Reunion92Name, Guyane92Name};
+    private final static String FortMarigotName = tr("Guadeloupe Fort-Marigot 1949");
+    private final static String SainteAnneName = tr("Guadeloupe Ste-Anne 1948");
+    private final static String MartiniqueName = tr("Martinique Fort Desaix 1952");
+    private final static String Reunion92Name = tr("Reunion RGR92");
+    private final static String Guyane92Name = tr("Guyane RGFG95");
+    private final static String[] utmGeodesicsNames = { FortMarigotName, SainteAnneName, MartiniqueName, Reunion92Name, Guyane92Name};
 
-    private Bounds FortMarigotBounds = new Bounds( new LatLon(17.6,-63.25), new LatLon(18.5,-62.5));
-    private Bounds SainteAnneBounds = new Bounds( new LatLon(15.8,-61.9), new LatLon(16.6,-60.9));
-    private Bounds MartiniqueBounds = new Bounds( new LatLon(14.25,-61.25), new LatLon(15.025,-60.725));
-    private Bounds ReunionBounds = new Bounds( new LatLon(-25.92,37.58), new LatLon(-10.6, 58.27));
-    private Bounds GuyaneBounds = new Bounds( new LatLon(2.16 , -54.0), new LatLon(9.06 , -49.62));
-    private Bounds[] utmBounds = { FortMarigotBounds, SainteAnneBounds, MartiniqueBounds, ReunionBounds, GuyaneBounds};
+    private final static Bounds FortMarigotBounds = new Bounds( new LatLon(17.6,-63.25), new LatLon(18.5,-62.5));
+    private final static Bounds SainteAnneBounds = new Bounds( new LatLon(15.8,-61.9), new LatLon(16.6,-60.9));
+    private final static Bounds MartiniqueBounds = new Bounds( new LatLon(14.25,-61.25), new LatLon(15.025,-60.725));
+    private final static Bounds ReunionBounds = new Bounds( new LatLon(-25.92,37.58), new LatLon(-10.6, 58.27));
+    private final static Bounds GuyaneBounds = new Bounds( new LatLon(2.16 , -54.0), new LatLon(9.06 , -49.62));
+    private final static Bounds[] utmBounds = { FortMarigotBounds, SainteAnneBounds, MartiniqueBounds, ReunionBounds, GuyaneBounds };
 
-    private String FortMarigotEPSG = "EPSG::2969";
-    private String SainteAnneEPSG = "EPSG::2970";
-    private String MartiniqueEPSG = "EPSG::2973";
-    private String ReunionEPSG = "EPSG::2975";
-    private String GuyaneEPSG = "EPSG::2972";
-    private String[] utmEPSGs = { FortMarigotEPSG, SainteAnneEPSG, MartiniqueEPSG, ReunionEPSG, GuyaneEPSG};
+    private final static Integer FortMarigotEPSG = 2969;
+    private final static Integer SainteAnneEPSG = 2970;
+    private final static Integer MartiniqueEPSG = 2973;
+    private final static Integer ReunionEPSG = 2975;
+    private final static Integer GuyaneEPSG = 2972;
+    private final static Integer[] utmEPSGs = { FortMarigotEPSG, SainteAnneEPSG, MartiniqueEPSG, ReunionEPSG, GuyaneEPSG };
 
-    /**
-     * false east in meters (constant)
-     */
-    private static final double Xs = 500000.0;
-    /**
-     * false north in meters (0 in northern hemisphere, 10000000 in southern
-     * hemisphere)
-     */
-    private static double Ys = 0;
-    /**
-     * origin meridian longitude
-     */
-    protected double lg0;
+    private final static Datum FortMarigotDatum = new ThreeParameterDatum("FortMarigot Datum", null, Ellipsoid.hayford, 136.596, 248.148, -429.789);
+    private final static Datum SainteAnneDatum = new SevenParameterDatum("SainteAnne Datum", null, Ellipsoid.hayford, -472.29, -5.63, -304.12, 0.4362, -0.8374, 0.2563, 1.8984);
+    private final static Datum MartiniqueDatum = new SevenParameterDatum("Martinique Datum", null, Ellipsoid.hayford, 126.926, 547.939, 130.409, -2.78670, 5.16124, -0.85844, 13.82265);
+    private final static Datum ReunionDatum = GRS80Datum.INSTANCE;
+    private final static Datum GuyaneDatum = GRS80Datum.INSTANCE;
+    private final static Datum[] utmDatums = { FortMarigotDatum, SainteAnneDatum, MartiniqueDatum, ReunionDatum, GuyaneDatum };
+
+    private final static int[] utmZones = { 20, 20, 20, 40, 22 };
+
     /**
      * UTM zone (from 1 to 60)
@@ -69 +69 @@
     public static final int DEFAULT_GEODESIC = 0;
 
-    public static int currentGeodesic = DEFAULT_GEODESIC;
+    public int currentGeodesic;
 
-    /**
-     * 7 parameters transformation
-     */
-    private static double tx = 0.0;
-    private static double ty = 0.0;
-    private static double tz = 0.0;
-    private static double rx = 0;
-    private static double ry = 0;
-    private static double rz = 0;
-    private static double scaleDiff = 0;
-    /**
-     * precision in iterative schema
-     */
-    public static final double epsilon = 1e-11;
 
-    private void refresh7ParametersTranslation() {
-        if (currentGeodesic == 0) { // UTM_20N_Guadeloupe_Fort_Marigot
-            set7ParametersTranslation(new double[]{136.596, 248.148, -429.789},
-                    new double[]{0, 0, 0},
-                    0,
-                    true, 20);
-        } else if (currentGeodesic == 1) { // UTM_20N_Guadeloupe_Ste_Anne
-            set7ParametersTranslation(new double[]{-472.29, -5.63, -304.12},
-                    new double[]{0.4362, -0.8374, 0.2563},
-                    1.8984E-6,
-                    true, 20);
-        } else if (currentGeodesic == 2) { // UTM_20N_Martinique_Fort_Desaix
-            set7ParametersTranslation(new double[]{126.926, 547.939, 130.409},
-                    new double[]{-2.78670, 5.16124,  -0.85844},
-                    13.82265E-6
-                    , true, 20);
-        } else if (currentGeodesic == 3) { // UTM_40S_Reunion_RGR92 (translation only required for re-projections from Gauss-Laborde)
-            set7ParametersTranslation(new double[]{789.524, -626.486, -89.904},
-                    new double[]{0.6006, 76.7946,  -10.5788},
-                    -32.3241E-6
-                    , false, 40);
-        } else if (currentGeodesic == 4) { // UTM_22N_Guyane_RGFG95 (translation only required for re-projections from CSG67)
-            set7ParametersTranslation(new double[]{-193.066 , 236.993, 105.447},
-                    new double[]{0.4814, -0.8074,  0.1276},
-                    1.5649E-6
-                    , true, 22);
-        }
+    public UTM_France_DOM() {
+        updateParameters(DEFAULT_GEODESIC);
+    }
+    
+    public void updateParameters(int currentGeodesic) {
+        this.currentGeodesic = currentGeodesic;
+        datum = utmDatums[currentGeodesic];
+        ellps = datum.getEllipsoid();
+        proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
+        isNorth = currentGeodesic != 3;
+        zone = utmZones[currentGeodesic];
+        x_0 = 500000;
+        y_0 = isNorth ? 0.0 : 10000000.0;
+        lon_0 = 6 * zone - 183;
+        k_0 = 0.9996;
+    }
+    
+    public int getCurrentGeodesic() {
+        return currentGeodesic;
     }
 
-    private void set7ParametersTranslation(double[] translation, double[] rotation, double scalediff, boolean north, int utmZone) {
-        tx = translation[0];
-        ty = translation[1];
-        tz = translation[2];
-        rx = rotation[0]/206264.806247096355; // seconds to radian
-        ry = rotation[1]/206264.806247096355;
-        rz = rotation[2]/206264.806247096355;
-        scaleDiff = scalediff;
-        isNorth = north;
-        Ys = isNorth ? 0.0 : 10000000.0;
-        zone = utmZone;
+    @Override 
+    public String toString() {
+        return tr("UTM France (DOM)");
     }
 
-    public EastNorth latlon2eastNorth(LatLon p) {
-        if (currentGeodesic < 3 ) {
-            // translate ellipsoid GRS80 (WGS83) => reference ellipsoid geographic
-            LatLon geo = GRS802Hayford(p);
-            // reference ellipsoid geographic => UTM projection
-            return MTProjection(geo, Ellipsoid.hayford.a, Ellipsoid.hayford.e);
-        } else { // UTM_40S_Reunion_RGR92 or UTM_22N_Guyane_RGFG95
-            LatLon geo = new LatLon(Math.toRadians(p.lat()), Math.toRadians(p.lon()));
-            return MTProjection(geo, Ellipsoid.GRS80.a, Ellipsoid.GRS80.e);
-        }
-    }
-
-    /**
-     * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to UTM
-     * geographic, (ellipsoid Hayford)
-     */
-    private LatLon GRS802Hayford(LatLon wgs) {
-        double lat = Math.toRadians(wgs.lat()); // degree to radian
-        double lon = Math.toRadians(wgs.lon());
-        // WGS84 geographic => WGS84 cartesian
-        double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
-        double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
-        double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
-        double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
-        // translation
-        double coord[] = invSevenParametersTransformation(X, Y, Z);
-        // UTM cartesian => UTM geographic
-        return Geographic(coord[0], coord[1], coord[2], Ellipsoid.hayford);
-    }
-
-    /**
-     * initializes from cartesian coordinates
-     *
-     * @param X
-     *            1st coordinate in meters
-     * @param Y
-     *            2nd coordinate in meters
-     * @param Z
-     *            3rd coordinate in meters
-     * @param ell
-     *            reference ellipsoid
-     */
-    private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
-        double norm = Math.sqrt(X * X + Y * Y);
-        double lg = 2.0 * Math.atan(Y / (X + norm));
-        double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
-        double delta = 1.0;
-        while (delta > epsilon) {
-            double s2 = Math.sin(lt);
-            s2 *= s2;
-            double l = Math.atan((Z / norm)
-                    / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
-            delta = Math.abs(l - lt);
-            lt = l;
-        }
-        double s2 = Math.sin(lt);
-        s2 *= s2;
-        // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
-        return new LatLon(lt, lg);
-    }
-
-    /**
-     * initalizes from geographic coordinates
-     *
-     * @param coord geographic coordinates triplet
-     * @param a reference ellipsoid long axis
-     * @param e reference ellipsoid excentricity
-     */
-    private EastNorth MTProjection(LatLon coord, double a, double e) {
-        double n = 0.9996 * a;
-        Ys = (coord.lat() >= 0.0) ? 0.0 : 10000000.0;
-        double r6d = Math.PI / 30.0;
-        //zone = (int) Math.floor((coord.lon() + Math.PI) / r6d) + 1;
-        lg0 = r6d * (zone - 0.5) - Math.PI;
-        double e2 = e * e;
-        double e4 = e2 * e2;
-        double e6 = e4 * e2;
-        double e8 = e4 * e4;
-        double C[] = {
-                1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
-                e2/8.0 - e4/96.0 - 9.0*e6/1024.0 - 901.0*e8/184320.0,
-                13.0*e4/768.0 + 17.0*e6/5120.0 - 311.0*e8/737280.0,
-                61.0*e6/15360.0 + 899.0*e8/430080.0,
-                49561.0*e8/41287680.0
-        };
-        double s = e * Math.sin(coord.lat());
-        double l = Math.log(Math.tan(Math.PI/4.0 + coord.lat()/2.0) *
-                Math.pow((1.0 - s) / (1.0 + s), e/2.0));
-        double phi = Math.asin(Math.sin(coord.lon() - lg0) /
-                ((Math.exp(l) + Math.exp(-l)) / 2.0));
-        double ls = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
-        double lambda = Math.atan(((Math.exp(l) - Math.exp(-l)) / 2.0) /
-                Math.cos(coord.lon() - lg0));
-
-        double north = C[0] * lambda;
-        double east = C[0] * ls;
-        for(int k = 1; k < 5; k++) {
-            double r = 2.0 * k * lambda;
-            double m = 2.0 * k * ls;
-            double em = Math.exp(m);
-            double en = Math.exp(-m);
-            double sr = Math.sin(r)/2.0 * (em + en);
-            double sm = Math.cos(r)/2.0 * (em - en);
-            north += C[k] * sr;
-            east += C[k] * sm;
-        }
-        east *= n;
-        east += Xs;
-        north *= n;
-        north += Ys;
-        return new EastNorth(east, north);
-    }
-
-    public LatLon eastNorth2latlon(EastNorth p) {
-        if (currentGeodesic < 3) {
-            MTProjection(p.east(), p.north(), zone, isNorth);
-            LatLon geo = Geographic(p, Ellipsoid.hayford.a, Ellipsoid.hayford.e, 0.0 /* z */);
-
-            // reference ellipsoid geographic => reference ellipsoid cartesian
-            double N = Ellipsoid.hayford.a / (Math.sqrt(1.0 - Ellipsoid.hayford.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
-            double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
-            double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
-            double Z = (N * (1.0-Ellipsoid.hayford.e2) /*+ h*/) * Math.sin(geo.lat());
-            // translation
-            double coord[] = sevenParametersTransformation(X, Y, Z);
-            // WGS84 cartesian => WGS84 geographic
-            LatLon wgs = cart2LatLon(coord[0], coord[1], coord[2], Ellipsoid.GRS80);
-            return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
-        } else {
-            // UTM_40S_Reunion_RGR92 or UTM_22N_Guyane_RGFG95
-            LatLon geo = Geographic(p, Ellipsoid.GRS80.a, Ellipsoid.GRS80.e, 0.0 /* z */);
-            double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
-            double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
-            double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
-            double Z = (N * (1.0-Ellipsoid.GRS80.e2) /*+ h*/) * Math.sin(geo.lat());
-            LatLon wgs = cart2LatLon(X, Y, Z, Ellipsoid.GRS80);
-            return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
-        }
-    }
-
-    /**
-     * initializes new projection coordinates (in north hemisphere)
-     *
-     * @param east east from origin in meters
-     * @param north north from origin in meters
-     * @param zone zone number (from 1 to 60)
-     * @param isNorth true in north hemisphere, false in south hemisphere
-     */
-    private void MTProjection(double east, double north, int zone, boolean isNorth) {
-        Ys = isNorth ? 0.0 : 10000000.0;
-        double r6d = Math.PI / 30.0;
-        lg0 = r6d * (zone - 0.5) - Math.PI;
-    }
-
-    public double scaleFactor() {
-        return 1/Math.PI/2;
-    }
-
-    /**
-     * initalizes from projected coordinates (Mercator transverse projection)
-     *
-     * @param coord projected coordinates pair
-     * @param e reference ellipsoid excentricity
-     * @param a reference ellipsoid long axis
-     * @param z altitude in meters
-     */
-    private LatLon Geographic(EastNorth coord, double a, double e, double z) {
-        double n = 0.9996 * a;
-        double e2 = e * e;
-        double e4 = e2 * e2;
-        double e6 = e4 * e2;
-        double e8 = e4 * e4;
-        double C[] = {
-                1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
-                e2/8.0 + e4/48.0 + 7.0*e6/2048.0 + e8/61440.0,
-                e4/768.0 + 3.0*e6/1280.0 + 559.0*e8/368640.0,
-                17.0*e6/30720.0 + 283.0*e8/430080.0,
-                4397.0*e8/41287680.0
-        };
-        double l = (coord.north() - Ys) / (n * C[0]);
-        double ls = (coord.east() - Xs) / (n * C[0]);
-        double l0 = l;
-        double ls0 = ls;
-        for(int k = 1; k < 5; k++) {
-            double r = 2.0 * k * l0;
-            double m = 2.0 * k * ls0;
-            double em = Math.exp(m);
-            double en = Math.exp(-m);
-            double sr = Math.sin(r)/2.0 * (em + en);
-            double sm = Math.cos(r)/2.0 * (em - en);
-            l -= C[k] * sr;
-            ls -= C[k] * sm;
-        }
-        double lon = lg0 + Math.atan(((Math.exp(ls) - Math.exp(-ls)) / 2.0) /
-                Math.cos(l));
-        double phi = Math.asin(Math.sin(l) /
-                ((Math.exp(ls) + Math.exp(-ls)) / 2.0));
-        l = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
-        double lat = 2.0 * Math.atan(Math.exp(l)) - Math.PI / 2.0;
-        double lt0;
-        do {
-            lt0 = lat;
-            double s = e * Math.sin(lat);
-            lat = 2.0 * Math.atan(Math.pow((1.0 + s) / (1.0 - s), e/2.0) *
-                    Math.exp(l)) - Math.PI / 2.0;
-        }
-        while(Math.abs(lat-lt0) >= epsilon);
-        //h = z;
-
-        return new LatLon(lat, lon);
-    }
-
-    /**
-     * initializes from cartesian coordinates
-     *
-     * @param X 1st coordinate in meters
-     * @param Y 2nd coordinate in meters
-     * @param Z 3rd coordinate in meters
-     * @param ell reference ellipsoid
-     */
-    private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) {
-        double[] XYZ = {X, Y, Z};
-        LatLon coord = ell.cart2LatLon(XYZ, epsilon);
-        return new LatLon(Math.toRadians(coord.lat()), Math.toRadians(coord.lon()));
-    }
-
-    /**
-     * 7 parameters transformation
-     * @param coord X, Y, Z in array
-     * @return transformed X, Y, Z in array
-     */
-    private double[] sevenParametersTransformation(double Xa, double Ya, double Za){
-        double Xb = tx + Xa*(1+scaleDiff) + Za*ry - Ya*rz;
-        double Yb = ty + Ya*(1+scaleDiff) + Xa*rz - Za*rx;
-        double Zb = tz + Za*(1+scaleDiff) + Ya*rx - Xa*ry;
-        return new double[]{Xb, Yb, Zb};
-    }
-
-    /**
-     * 7 parameters inverse transformation
-     * @param coord X, Y, Z in array
-     * @return transformed X, Y, Z in array
-     */
-    private double[] invSevenParametersTransformation(double Xa, double Ya, double Za){
-        double Xb = (1-scaleDiff)*(-tx + Xa + ((-tz+Za)*(-ry) - (-ty+Ya)*(-rz)));
-        double Yb = (1-scaleDiff)*(-ty + Ya + ((-tx+Xa)*(-rz) - (-tz+Za)*(-rx)));
-        double Zb = (1-scaleDiff)*(-tz + Za + ((-ty+Ya)*(-rx) - (-tx+Xa)*(-ry)));
-        return new double[]{Xb, Yb, Zb};
-    }
-
+    @Override
     public String getCacheDirectoryName() {
         return this.toString();
     }
 
-    /**
-     * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
-     */
-    public double getDefaultZoomInPPD() {
-        // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
-        return 10.0;
-    }
-
+    @Override
     public Bounds getWorldBoundsLatLon() {
         return utmBounds[currentGeodesic];
     }
 
-    public String toCode() {
+    @Override
+    public Integer getEpsgCode() {
         return utmEPSGs[currentGeodesic];
     }
@@ -401 +116 @@
     public int hashCode() {
         return getClass().getName().hashCode()+currentGeodesic; // our only real variable
-    }
-
-    @Override public String toString() {
-        return (tr("UTM France (DOM)"));
-    }
-
-    public int getCurrentGeodesic() {
-        return currentGeodesic;
     }
 
@@ -426 +133 @@
     }
 
+    @Override
     public Collection<String> getPreferences(JPanel p) {
         Object prefcb = p.getComponent(2);
@@ -431 +139 @@
             return null;
         currentGeodesic = ((JComboBox)prefcb).getSelectedIndex();
-        refresh7ParametersTranslation();
         return Collections.singleton(Integer.toString(currentGeodesic+1));
     }
@@ -437 +144 @@
     @Override
     public String[] allCodes() {
-        return utmEPSGs;
+        String[] res = new String[utmEPSGs.length];
+        for (int i=0; i<utmEPSGs.length; ++i) {
+            res[i] = "EPSG:"+utmEPSGs[i];
+        }
+        return res;
     }
 
+    @Override
     public Collection<String> getPreferencesFromCode(String code) {
         for (int i=0; i < utmEPSGs.length; i++ )
-            if (utmEPSGs[i].endsWith(code))
+            if (("EPSG:"+utmEPSGs[i]).equals(code))
                 return Collections.singleton(Integer.toString(i+1));
         return null;
     }
 
+    @Override
     public void setPreferences(Collection<String> args) {
-        currentGeodesic = DEFAULT_GEODESIC;
+        int currentGeodesic = DEFAULT_GEODESIC;
         if (args != null) {
             try {
@@ -461 +174 @@
             } catch(NumberFormatException e) {}
         }
-        refresh7ParametersTranslation();
+        updateParameters(currentGeodesic);
     }
-
 }

trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java

@@ -1 +1 @@
 // License: GPL. For details, see LICENSE file.
-package org.openstreetmap.josm.data.projection;
-
-import org.openstreetmap.josm.data.coor.EastNorth;
-import org.openstreetmap.josm.data.coor.LatLon;
+package org.openstreetmap.josm.data.projection.proj;
+
+import static java.lang.Math.*;
+
+import static org.openstreetmap.josm.tools.I18n.tr;
+
+import org.openstreetmap.josm.data.projection.Ellipsoid;
 
 /**
- * This is a base class to do projections based on Transverse Mercator projection.
+ * Transverse Mercator projection.
  *
  * @author Dirk Stöcker
  * code based on JavaScript from Chuck Taylor
  * 
- * NOTE: Uses polygon approximation to translate to WGS84.
  */
-public abstract class TransverseMercator implements Projection {
-
-    private final static double UTMScaleFactor = 0.9996;
-
-    private double UTMCentralMeridianRad = 0;
-    private double offsetEastMeters = 500000;
-    private double offsetNorthMeters = 0;
-
-
-    protected void setProjectionParameters(double centralMeridianDegress, double offsetEast, double offsetNorth)
-    {
-        UTMCentralMeridianRad = Math.toRadians(centralMeridianDegress);
-        offsetEastMeters = offsetEast;
-        offsetNorthMeters = offsetNorth;
-    }
-
-    /*
-     * ArcLengthOfMeridian
-     *
-     * Computes the ellipsoidal distance from the equator to a point at a
-     * given latitude.
-     *
-     * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
-     * GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
-     *
-     * Inputs:
-     *     phi - Latitude of the point, in radians.
-     *
-     * Globals:
-     *     Ellipsoid.GRS80.a - Ellipsoid model major axis.
-     *     Ellipsoid.GRS80.b - Ellipsoid model minor axis.
-     *
-     * Returns:
-     *     The ellipsoidal distance of the point from the equator, in meters.
-     *
-     */
-    private double ArcLengthOfMeridian(double phi)
-    {
-        /* Precalculate n */
-        double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
-
-        /* Precalculate alpha */
-        double alpha = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
-        * (1.0 + (Math.pow (n, 2.0) / 4.0) + (Math.pow (n, 4.0) / 64.0));
-
-        /* Precalculate beta */
-        double beta = (-3.0 * n / 2.0) + (9.0 * Math.pow (n, 3.0) / 16.0)
-        + (-3.0 * Math.pow (n, 5.0) / 32.0);
-
-        /* Precalculate gamma */
-        double gamma = (15.0 * Math.pow (n, 2.0) / 16.0)
-        + (-15.0 * Math.pow (n, 4.0) / 32.0);
-
-        /* Precalculate delta */
-        double delta = (-35.0 * Math.pow (n, 3.0) / 48.0)
-        + (105.0 * Math.pow (n, 5.0) / 256.0);
-
-        /* Precalculate epsilon */
-        double epsilon = (315.0 * Math.pow (n, 4.0) / 512.0);
-
-        /* Now calculate the sum of the series and return */
-        return alpha
-        * (phi + (beta * Math.sin (2.0 * phi))
-                + (gamma * Math.sin (4.0 * phi))
-                + (delta * Math.sin (6.0 * phi))
-                + (epsilon * Math.sin (8.0 * phi)));
-    }
-
-    /*
-     * FootpointLatitude
-     *
-     * Computes the footpoint latitude for use in converting transverse
-     * Mercator coordinates to ellipsoidal coordinates.
-     *
-     * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
-     *   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
-     *
-     * Inputs:
-     *   y - The UTM northing coordinate, in meters.
-     *
-     * Returns:
-     *   The footpoint latitude, in radians.
-     *
-     */
-    private double FootpointLatitude(double y)
-    {
-        /* Precalculate n (Eq. 10.18) */
-        double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
-
-        /* Precalculate alpha_ (Eq. 10.22) */
-        /* (Same as alpha in Eq. 10.17) */
-        double alpha_ = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
-        * (1 + (Math.pow (n, 2.0) / 4) + (Math.pow (n, 4.0) / 64));
-
-        /* Precalculate y_ (Eq. 10.23) */
-        double y_ = y / alpha_;
-
-        /* Precalculate beta_ (Eq. 10.22) */
-        double beta_ = (3.0 * n / 2.0) + (-27.0 * Math.pow (n, 3.0) / 32.0)
-        + (269.0 * Math.pow (n, 5.0) / 512.0);
-
-        /* Precalculate gamma_ (Eq. 10.22) */
-        double gamma_ = (21.0 * Math.pow (n, 2.0) / 16.0)
-        + (-55.0 * Math.pow (n, 4.0) / 32.0);
-
-        /* Precalculate delta_ (Eq. 10.22) */
-        double delta_ = (151.0 * Math.pow (n, 3.0) / 96.0)
-        + (-417.0 * Math.pow (n, 5.0) / 128.0);
-
-        /* Precalculate epsilon_ (Eq. 10.22) */
-        double epsilon_ = (1097.0 * Math.pow (n, 4.0) / 512.0);
-
-        /* Now calculate the sum of the series (Eq. 10.21) */
-        return y_ + (beta_ * Math.sin (2.0 * y_))
-        + (gamma_ * Math.sin (4.0 * y_))
-        + (delta_ * Math.sin (6.0 * y_))
-        + (epsilon_ * Math.sin (8.0 * y_));
-    }
-
-    /*
-     * MapLatLonToXY
-     *
+public class TransverseMercator implements Proj {
+
+    protected double a, b;
+
+    public TransverseMercator(Ellipsoid ellps) {
+        this.a = ellps.a;
+        this.b = ellps.b;
+    }
+    
+    @Override
+    public String getName() {
+        return tr("Transverse Mercator");
+    }
+
+    @Override
+    public String getProj4Id() {
+        return "tmerc";
+    }
+
+    /**
      * Converts a latitude/longitude pair to x and y coordinates in the
      * Transverse Mercator projection.  Note that Transverse Mercator is not
@@ -141 +41 @@
      * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
      * GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
-     *
-     * Inputs:
-     *    phi - Latitude of the point, in radians.
-     *    lambda - Longitude of the point, in radians.
-     *    lambda0 - Longitude of the central meridian to be used, in radians.
-     *
-     * Outputs:
-     *    xy - A 2-element array containing the x and y coordinates
-     *         of the computed point.
-     *
-     * Returns:
-     *    The function does not return a value.
-     *
-     */
-    public EastNorth mapLatLonToXY(double phi, double lambda, double lambda0)
-    {
+     * 
+     * @param phi Latitude of the point, in radians
+     * @param lambda Longitude of the point, in radians
+     * @return A 2-element array containing the x and y coordinates
+     *         of the computed point
+     */
+    @Override
+    public double[] project(double phi, double lambda) {
+        
         /* Precalculate ep2 */
-        double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0)) / Math.pow (Ellipsoid.GRS80.b, 2.0);
+        double ep2 = (pow(a, 2.0) - pow(b, 2.0)) / pow(b, 2.0);
 
         /* Precalculate nu2 */
-        double nu2 = ep2 * Math.pow (Math.cos (phi), 2.0);
-
-        /* Precalculate N */
-        double N = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nu2));
+        double nu2 = ep2 * pow(cos(phi), 2.0);
+
+        /* Precalculate N / a */
+        double N_a = a / (b * sqrt(1 + nu2));
 
         /* Precalculate t */
-        double t = Math.tan (phi);
+        double t = tan(phi);
         double t2 = t * t;
 
         /* Precalculate l */
-        double l = lambda - lambda0;
+        double l = lambda;
 
         /* Precalculate coefficients for l**n in the equations below
@@ -191 +84 @@
         double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
 
-        return new EastNorth(
+        return new double[] {
                 /* Calculate easting (x) */
-                N * Math.cos (phi) * l
-                + (N / 6.0 * Math.pow (Math.cos (phi), 3.0) * l3coef * Math.pow (l, 3.0))
-                + (N / 120.0 * Math.pow (Math.cos (phi), 5.0) * l5coef * Math.pow (l, 5.0))
-                + (N / 5040.0 * Math.pow (Math.cos (phi), 7.0) * l7coef * Math.pow (l, 7.0)),
+                N_a * cos(phi) * l
+                + (N_a / 6.0 * pow(cos(phi), 3.0) * l3coef * pow(l, 3.0))
+                + (N_a / 120.0 * pow(cos(phi), 5.0) * l5coef * pow(l, 5.0))
+                + (N_a / 5040.0 * pow(cos(phi), 7.0) * l7coef * pow(l, 7.0)),
                 /* Calculate northing (y) */
-                ArcLengthOfMeridian (phi)
-                + (t / 2.0 * N * Math.pow (Math.cos (phi), 2.0) * Math.pow (l, 2.0))
-                + (t / 24.0 * N * Math.pow (Math.cos (phi), 4.0) * l4coef * Math.pow (l, 4.0))
-                + (t / 720.0 * N * Math.pow (Math.cos (phi), 6.0) * l6coef * Math.pow (l, 6.0))
-                + (t / 40320.0 * N * Math.pow (Math.cos (phi), 8.0) * l8coef * Math.pow (l, 8.0)));
-    }
-
-    /*
-     * MapXYToLatLon
-     *
+                ArcLengthOfMeridian (phi) / a
+                + (t / 2.0 * N_a * pow(cos(phi), 2.0) * pow(l, 2.0))
+                + (t / 24.0 * N_a * pow(cos(phi), 4.0) * l4coef * pow(l, 4.0))
+                + (t / 720.0 * N_a * pow(cos(phi), 6.0) * l6coef * pow(l, 6.0))
+                + (t / 40320.0 * N_a * pow(cos(phi), 8.0) * l8coef * pow(l, 8.0)) };
+    }
+
+    /**
      * Converts x and y coordinates in the Transverse Mercator projection to
      * a latitude/longitude pair.  Note that Transverse Mercator is not
@@ -214 +105 @@
      * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
      *   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
-     *
-     * Inputs:
-     *   x - The easting of the point, in meters.
-     *   y - The northing of the point, in meters.
-     *   lambda0 - Longitude of the central meridian to be used, in radians.
-     *
-     * Outputs:
-     *   philambda - A 2-element containing the latitude and longitude
-     *               in radians.
-     *
-     * Returns:
-     *   The function does not return a value.
-     *
+     * 
      * Remarks:
      *   The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
@@ -234 +113 @@
      *   x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
      *   to optimize computations.
-     *
-     */
-    public LatLon mapXYToLatLon(double x, double y, double lambda0)
-    {
+     * 
+     * @param x The easting of the point, in meters, divided by the semi major axis of the ellipsoid
+     * @param y The northing of the point, in meters, divided by the semi major axis of the ellipsoid
+     * @return A 2-element containing the latitude and longitude
+     *               in radians
+     */
+    @Override
+    public double[] invproject(double x, double y) {
         /* Get the value of phif, the footpoint latitude. */
-        double phif = FootpointLatitude (y);
+        double phif = footpointLatitude(y);
 
         /* Precalculate ep2 */
-        double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0))
-        / Math.pow (Ellipsoid.GRS80.b, 2.0);
-
-        /* Precalculate cos (phif) */
-        double cf = Math.cos (phif);
+        double ep2 = (a*a - b*b)
+        / (b*b);
+
+        /* Precalculate cos(phif) */
+        double cf = cos(phif);
 
         /* Precalculate nuf2 */
-        double nuf2 = ep2 * Math.pow (cf, 2.0);
-
-        /* Precalculate Nf and initialize Nfpow */
-        double Nf = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nuf2));
-        double Nfpow = Nf;
+        double nuf2 = ep2 * pow(cf, 2.0);
+
+        /* Precalculate Nf / a and initialize Nfpow */
+        double Nf_a = a / (b * sqrt(1 + nuf2));
+        double Nfpow = Nf_a;
 
         /* Precalculate tf */
-        double tf = Math.tan (phif);
+        double tf = tan(phif);
         double tf2 = tf * tf;
         double tf4 = tf2 * tf2;
@@ -264 +147 @@
         double x1frac = 1.0 / (Nfpow * cf);
 
-        Nfpow *= Nf;   /* now equals Nf**2) */
+        Nfpow *= Nf_a;   /* now equals Nf**2) */
         double x2frac = tf / (2.0 * Nfpow);
 
-        Nfpow *= Nf;   /* now equals Nf**3) */
+        Nfpow *= Nf_a;   /* now equals Nf**3) */
         double x3frac = 1.0 / (6.0 * Nfpow * cf);
 
-        Nfpow *= Nf;   /* now equals Nf**4) */
+        Nfpow *= Nf_a;   /* now equals Nf**4) */
         double x4frac = tf / (24.0 * Nfpow);
 
-        Nfpow *= Nf;   /* now equals Nf**5) */
+        Nfpow *= Nf_a;   /* now equals Nf**5) */
         double x5frac = 1.0 / (120.0 * Nfpow * cf);
 
-        Nfpow *= Nf;   /* now equals Nf**6) */
+        Nfpow *= Nf_a;   /* now equals Nf**6) */
         double x6frac = tf / (720.0 * Nfpow);
 
-        Nfpow *= Nf;   /* now equals Nf**7) */
+        Nfpow *= Nf_a;   /* now equals Nf**7) */
         double x7frac = 1.0 / (5040.0 * Nfpow * cf);
 
-        Nfpow *= Nf;   /* now equals Nf**8) */
+        Nfpow *= Nf_a;   /* now equals Nf**8) */
         double x8frac = tf / (40320.0 * Nfpow);
 
@@ -295 +178 @@
         double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2);
 
-        return new LatLon(
+        return new double[] {
                 /* Calculate latitude */
-                Math.toDegrees(
                         phif + x2frac * x2poly * (x * x)
-                        + x4frac * x4poly * Math.pow (x, 4.0)
-                        + x6frac * x6poly * Math.pow (x, 6.0)
-                        + x8frac * x8poly * Math.pow (x, 8.0)),
-                        Math.toDegrees(
-                                /* Calculate longitude */
-                                lambda0 + x1frac * x
-                                + x3frac * x3poly * Math.pow (x, 3.0)
-                                + x5frac * x5poly * Math.pow (x, 5.0)
-                                + x7frac * x7poly * Math.pow (x, 7.0)));
-    }
-
-    @Override
-    public EastNorth latlon2eastNorth(LatLon p) {
-        EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), UTMCentralMeridianRad);
-        return new EastNorth(a.east() * UTMScaleFactor + offsetEastMeters, a.north() * UTMScaleFactor + offsetNorthMeters);
-    }
-
-    @Override
-    public LatLon eastNorth2latlon(EastNorth p) {
-        return mapXYToLatLon((p.east() - offsetEastMeters)/UTMScaleFactor, (p.north() - offsetNorthMeters)/UTMScaleFactor, UTMCentralMeridianRad);
-    }
-
-    @Override
-    public double getDefaultZoomInPPD() {
-        // this will set the map scaler to about 1000 m
-        return 10;
-    }
+                        + x4frac * x4poly * pow(x, 4.0)
+                        + x6frac * x6poly * pow(x, 6.0)
+                        + x8frac * x8poly * pow(x, 8.0),
+                        /* Calculate longitude */
+                        x1frac * x
+                        + x3frac * x3poly * pow(x, 3.0)
+                        + x5frac * x5poly * pow(x, 5.0)
+                        + x7frac * x7poly * pow(x, 7.0) };
+    }
+    
+    /**
+     * ArcLengthOfMeridian
+     *
+     * Computes the ellipsoidal distance from the equator to a point at a
+     * given latitude.
+     *
+     * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
+     * GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
+     * 
+     * @param phi Latitude of the point, in radians
+     * @return The ellipsoidal distance of the point from the equator
+     *         (in meters, divided by the semi major axis of the ellipsoid)
+     */
+    private double ArcLengthOfMeridian(double phi) {
+        /* Precalculate n */
+        double n = (a - b) / (a + b);
+
+        /* Precalculate alpha */
+        double alpha = ((a + b) / 2.0)
+            * (1.0 + (pow(n, 2.0) / 4.0) + (pow(n, 4.0) / 64.0));
+
+        /* Precalculate beta */
+        double beta = (-3.0 * n / 2.0) + (9.0 * pow(n, 3.0) / 16.0)
+            + (-3.0 * pow(n, 5.0) / 32.0);
+
+        /* Precalculate gamma */
+        double gamma = (15.0 * pow(n, 2.0) / 16.0)
+            + (-15.0 * pow(n, 4.0) / 32.0);
+
+        /* Precalculate delta */
+        double delta = (-35.0 * pow(n, 3.0) / 48.0)
+            + (105.0 * pow(n, 5.0) / 256.0);
+
+        /* Precalculate epsilon */
+        double epsilon = (315.0 * pow(n, 4.0) / 512.0);
+
+        /* Now calculate the sum of the series and return */
+        return alpha
+            * (phi + (beta * sin(2.0 * phi))
+                    + (gamma * sin(4.0 * phi))
+                    + (delta * sin(6.0 * phi))
+                    + (epsilon * sin(8.0 * phi)));
+    }
+        
+    /**
+     * FootpointLatitude
+     *
+     * Computes the footpoint latitude for use in converting transverse
+     * Mercator coordinates to ellipsoidal coordinates.
+     *
+     * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
+     *   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
+     * 
+     * @param y northing coordinate, in meters, divided by the semi major axis of the ellipsoid
+     * @return The footpoint latitude, in radians
+     */
+    private double footpointLatitude(double y) {
+        /* Precalculate n (Eq. 10.18) */
+        double n = (a - b) / (a + b);
+
+        /* Precalculate alpha_ (Eq. 10.22) */
+        /* (Same as alpha in Eq. 10.17) */
+        double alpha_ = ((a + b) / 2.0)
+            * (1 + (pow(n, 2.0) / 4) + (pow(n, 4.0) / 64));
+
+        /* Precalculate y_ (Eq. 10.23) */
+        double y_ = y / alpha_ * a;
+
+        /* Precalculate beta_ (Eq. 10.22) */
+        double beta_ = (3.0 * n / 2.0) + (-27.0 * pow(n, 3.0) / 32.0)
+            + (269.0 * pow(n, 5.0) / 512.0);
+
+        /* Precalculate gamma_ (Eq. 10.22) */
+        double gamma_ = (21.0 * pow(n, 2.0) / 16.0)
+            + (-55.0 * pow(n, 4.0) / 32.0);
+
+        /* Precalculate delta_ (Eq. 10.22) */
+        double delta_ = (151.0 * pow(n, 3.0) / 96.0)
+            + (-417.0 * pow(n, 5.0) / 128.0);
+
+        /* Precalculate epsilon_ (Eq. 10.22) */
+        double epsilon_ = (1097.0 * pow(n, 4.0) / 512.0);
+
+        /* Now calculate the sum of the series (Eq. 10.21) */
+        return y_ + (beta_ * sin(2.0 * y_))
+            + (gamma_ * sin(4.0 * y_))
+            + (delta_ * sin(6.0 * y_))
+            + (epsilon_ * sin(8.0 * y_));
+    }
+    
 }