| gvSIG desktop 1

package org.cresques.cts;

import junit.framework.Test;

import junit.framework.TestSuite;

public class AllTests {

	public static Test suite() {

		TestSuite suite = new TestSuite("Test for org.cresques.cts");

		//$JUnit-BEGIN$

		//$JUnit-END$

		return suite;

	}

}

#Translations for language [cs]

#Mon Nov 06 09:04:31 CET 2006

Aceptar=Budi\u017e

Aplicar=Pou\u017e\u00edt

Block_Size_=

Cancelar=Zru\u0161it

Compression_=

datum=Datum

Generate_Tfw_=

Interleave_=

Photometric_=

Progressive_=

projection=Projekce

reference_system=Prostorov\u00fd referen\u010dn\u00ed syst\u00e9m

zone=Z\u00f3na

#Translations for language [fr]

#Mon Oct 30 09:38:21 CET 2006

Aceptar=Accepter

Aplicar=Appliquer

Block_Size_=

Cancelar=Annuler

Compression_=

datum=

Generate_Tfw_=

Interleave_=

Photometric_=

Progressive_=

projection=

reference_system=

zone=

#Translations for language [de]

#Mon Oct 30 09:38:21 CET 2006

Aceptar=OK

Aplicar=Anwenden

Block_Size_=

Cancelar=Abbrechen

Compression_=

datum=

Generate_Tfw_=

Interleave_=

Photometric_=

Progressive_=

projection=

reference_system=

zone=

#Translations for language [eu]

#Mon Oct 30 09:38:21 CET 2006

Aceptar=Ados

Aplicar=Aplikatu

Block_Size_=

Cancelar=Utzi

Compression_=

datum=

Generate_Tfw_=

Interleave_=

Photometric_=

Progressive_=

projection=

reference_system=

zone=

#Translations for language [it]

#Tue Nov 07 12:30:01 CET 2006

Aceptar=Accetta

Aplicar=Applica

Block_Size_=Dimensioni blocco\:

Cancelar=Cancella

Compression_=Compressione

datum=Datum

Generate_Tfw_=Genera Tfw\:

Interleave_=Interlacciato\:

Photometric_=Fotometrica\:

Progressive_=Progressivo\:

projection=Proiezione

reference_system=Sistema di riferimento

zone=Fuso

#Translations for language [es]

#Mon Oct 30 09:38:21 CET 2006

Aceptar=Aceptar

Aplicar=Aplicar

Block_Size_=Tama\u00f1o de bloque\:

Cancelar=Cancelar

Compression_=Compresi\u00f3n\:

datum=Datum

Generate_Tfw_=Generar Tfw\:

Interleave_=Entrelazado\:

Photometric_=Fotom\u00e9trica\:

Progressive_=Progresivo\:

projection=Proyecci\u00f3n

reference_system=Sistema de Referencia

zone=Huso

#Translations for language [en]

#Mon Oct 30 09:38:21 CET 2006

Aceptar=Accept

Aplicar=Apply

Block_Size_=Block Size\:

Cancelar=Cancel

Compression_=Compression\:

datum=Datum

Generate_Tfw_=Generate Tfw\:

Interleave_=Interleave\:

Photometric_=Photometric\:

Progressive_=Progressive\:

projection=Projection

reference_system=Reference System

zone=Zone

#Translations for language [gl]

#Mon Oct 30 09:38:21 CET 2006

Aceptar=Aceptar

Aplicar=Aplicar

Block_Size_=

Cancelar=Cancelar

Compression_=

datum=

Generate_Tfw_=

Interleave_=

Photometric_=

Progressive_=

projection=

reference_system=

zone=

#Translations for language [ca]

#Mon Oct 30 09:38:21 CET 2006

Aceptar=Acceptar

Aplicar=Aplicar

Block_Size_=

Cancelar=Cancel\u00b7lar

Compression_=

datum=

Generate_Tfw_=

Interleave_=

Photometric_=

Progressive_=

projection=Projecci\u00f3

reference_system=Sistema de refer\u00e8ncia

zone=

#Translations for language [pt]

#Mon Oct 30 09:38:21 CET 2006

Aceptar=Aceitar

Aplicar=Aplicar

Block_Size_=

Cancelar=Cancelar

Compression_=

datum=

Generate_Tfw_=

Interleave_=

Photometric_=

Progressive_=

projection=

reference_system=

zone=

<?xml version="1.0" encoding="UTF-8" ?>

<!DOCTYPE log4j:configuration SYSTEM "log4j.dtd">

<log4j:configuration xmlns:log4j="http://jakarta.apache.org/log4j/">

	<appender name="CONSOLE" class="org.apache.log4j.ConsoleAppender">

		<layout class="org.apache.log4j.PatternLayout">

			<param name="ConversionPattern" value="%d{HH:mm:ss,SSS} %-5p [%c{2}.%M()]\n  %m%n" />

		</layout>

	</appender>

	<category name="org.gvsig.tools">

		<priority value="DEBUG" />

	</category>

	<category name="org.cresques">

		<priority value="DEBUG" /> 

	</category>

	<category name="org.gvsig.fmap.crs">

		<priority value="DEBUG" /> 

	</category>

	<root>

		<priority value="INFO" />

		<appender-ref ref="CONSOLE" />

	</root>

</log4j:configuration>

/* Cresques Mapping Suite. Graphic Library for constructing mapping applications.

* Copyright (C) 2006 IVER T.I. and Generalitat Valenciana.

*

* This program is free software; you can redistribute it and/or

* modify it under the terms of the GNU General Public License

* as published by the Free Software Foundation; either version 2

* of the License, or (at your option) any later version.

*

* This program is distributed in the hope that it will be useful,

* but WITHOUT ANY WARRANTY; without even the implied warranty of

* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

* GNU General Public License for more details.

*

* You should have received a copy of the GNU General Public License

* along with this program; if not, write to the Free Software

* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.

*

* For more information, contact:

*

*  Generalitat Valenciana

*   Conselleria d'Infraestructures i Transport

*   Av. Blasco Ib??ez, 50

*   46010 VALENCIA

*   SPAIN

*

*      +34 963862235

*   gvsig@gva.es

*      www.gvsig.gva.es

*

*    or

*

*   IVER T.I. S.A

*   Salamanca 50

*   46005 Valencia

*   Spain

*

*   +34 963163400

*   dac@iver.es

*/

package org.cresques;

import java.util.Locale;

/**

* Bridge class to provide internationalization services to the library.

* It uses the gvsig-i18n library as a backend, and includes its

* necessary initialization.

* 

* @author Cesar Martinez Izquierdo

*

*/

public class Messages {

	/**

	 * Whether the class has been initialized

	 */

	private static boolean isInitialized = false;

	/**

	 * The name of the Java package containing this class

	 */

//	private static final String packageName = Messages.class.getName() ;

	/**

	 * Loads the translations in the dictionary. It initializes the backend

	 * gvsig-i18n library

	 *

	 */

	private static void init() {

		if (!org.gvsig.i18n.Messages.hasLocales()) {

			org.gvsig.i18n.Messages.addLocale(Locale.getDefault());

		}

		org.gvsig.i18n.Messages.addResourceFamily("org.cresques.resources.i18n.text", Messages.class.getClassLoader(), Messages.class.getClass().getName());

	}

	/**

	 * Gets the translation associated with the provided translation key.

	 * 

	 * @param key The translation key which identifies the target text

	 * @return The translation associated with the provided translation key.

	 */

	public static String getText(String key) {

		if (isInitialized==false) {

			init();

			isInitialized = true;

		}

		return org.gvsig.i18n.Messages.getText(key, Messages.class.getName());

	}

}

/*

 * Cresques Mapping Suite. Graphic Library for constructing mapping applications.

 *

 * Copyright (C) 2004-5.

 *

 * This program is free software; you can redistribute it and/or

 * modify it under the terms of the GNU General Public License

 * as published by the Free Software Foundation; either version 2

 * of the License, or (at your option) any later version.

 *

 * This program is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.

 *

 * For more information, contact:

 *

 * cresques@gmail.com

 */

package org.cresques.px;

import java.awt.geom.Point2D;

import java.awt.geom.Rectangle2D;

import java.text.DecimalFormat;

/**

 *	Clase que getiona el extent de una imagen

 *	

 *  @author Luis W.Sevilla (sevilla_lui@gva.es)

 */

public class Extent {

    Point2D min = null;

    Point2D max = null;

    /**

     * Constructor sin par?metros

     */

    public Extent() {

        min = new Point2D.Double(999999999.0, 999999999.0);

        max = new Point2D.Double(-999999999.0, -999999999.0);

    }

    /**

     * Constructor 

     * @param pt1	punto que representa la esquina superior izquierda

     * @param pt2	punto que representa la esquina inferior derecha

     */

    public Extent(Point2D pt1, Point2D pt2) {

        newExtent(pt1.getX(), pt1.getY(), pt2.getX(), pt2.getY());

    }

    /**

     * Contructor

     * @param x1 punto que representa la coordenada X de la esquina superior izquierda

     * @param y1 punto que representa la coordenada Y de la esquina superior izquierda

     * @param x2 punto que representa la coordenada X de la esquina inferior derecha

     * @param y2 punto que representa la coordenada Y de la esquina inferior derecha

     */

    public Extent(double x1, double y1, double x2, double y2) {

        newExtent(x1, y1, x2, y2);

    }

    /**

     * Constructor

     * @param r	Rectangulo 2D

     */

    public Extent(Rectangle2D r) {

        newExtent(r.getX(), r.getY(), r.getX() + r.getWidth(),

                  r.getY() + r.getHeight());

    }

    /**

     * Constructor de copia

     * @param ext	Objeto Extent

     */

    public Extent(Extent ext) {

        newExtent(ext.minX(), ext.minY(), ext.maxX(), ext.maxY());

    }

    /**

     * Crea un objeto extent identico y lo retorna

     * @return Objeto extent

     */

    public Object clone() {

        Extent e = (Extent) clone();

        e.min = (Point2D) min.clone();

        e.max = (Point2D) max.clone();

        return e;

    }

    private void newExtent(double x1, double y1, double x2, double y2) {

        min = new Point2D.Double(Math.min(x1, x2), Math.min(y1, y2));

        max = new Point2D.Double(Math.max(x1, x2), Math.max(y1, y2));

    }

    /**

     * Obtiene la coordenada X m?nima

     * @return valor de la coordenada X m?nima

     */

    public double minX() {

        return min.getX();

    }

    /**

     * Obtiene la coordenada Y m?nima

     * @return valor de la coordenada X m?nima

     */

    public double minY() {

        return min.getY();

    }

    /**

     * Obtiene la coordenada X m?xima

     * @return valor de la coordenada X m?xima

     */

    public double maxX() {

        return max.getX();

    }

    /**

     * Obtiene la coordenada Y m?xima

     * @return valor de la coordenada Y m?xima

     */

    public double maxY() {

        return max.getY();

    }

    /**

     * Obtiene el punto m?nimo

     * @return m?nimo

     */

    public Point2D getMin() {

        return min;

    }

    /**

     * Obtiene el punto m?ximo

     * @return m?ximo

     */

    public Point2D getMax() {

        return max;

    }

    public boolean isAt(Point2D pt) {

        if (pt.getX() < minX()) {

            return false;

        }

        if (pt.getX() > maxX()) {

            return false;

        }

        if (pt.getY() < minY()) {

            return false;

        }

        if (pt.getY() > maxY()) {

            return false;

        }

        return true;

    }

    public double width() {

        return Math.abs(maxX() - minX());

    }

    public double height() {

        return Math.abs(maxY() - minY());

    }

    /**

     * Verifica un punto, y modifica el extent si no est? incluido

     */

    public void add(Point2D pt) {

        if (pt == null) {

            return;

        }

        min.setLocation(Math.min(pt.getX(), minX()), Math.min(pt.getY(), minY()));

        max.setLocation(Math.max(pt.getX(), maxX()), Math.max(pt.getY(), maxY()));

    }

    public void add(Extent ext) {

        if (ext == null) {

            return;

        }

        min.setLocation(Math.min(ext.minX(), minX()),

                        Math.min(ext.minY(), minY()));

        max.setLocation(Math.max(ext.maxX(), maxX()),

                        Math.max(ext.maxY(), maxY()));

    }

    /**

     * Obtiene la escala

     * @param width	Ancho

     * @param height	Alto

     * @return

     */

    public double[] getScale(int width, int height) {

        return getScale((double) width, (double) height);

    }

    public double[] getScale(double width, double height) {

        double[] scale = new double[2];

        scale[0] = ((float) width) / width();

        scale[1] = ((float) height) / height();

        return scale;

    }

    public Rectangle2D toRectangle2D() {

        return new Rectangle2D.Double(minX(), minY(), width(), height());

    }

    public String toString() {

        DecimalFormat format = new DecimalFormat("####.000");

        return "Extent: (" + format.format(minX()) + "," +

               format.format(minY()) + "), (" + format.format(maxX()) + "," +

               format.format(maxY()) + ")";

    }

    public interface Has {

        public Extent getExtent();

    }

}

*.dfPackage

*.wmf

<html>

	<body>Pixel: Clases para dibujar.

</body>

</html>

/*

 * Cresques Mapping Suite. Graphic Library for constructing mapping applications.

 *

 * Copyright (C) 2004-5.

 *

 * This program is free software; you can redistribute it and/or

 * modify it under the terms of the GNU General Public License

 * as published by the Free Software Foundation; either version 2

 * of the License, or (at your option) any later version.

 *

 * This program is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.

 *

 * For more information, contact:

 *

 * cresques@gmail.com

 */

package org.cresques.cts;

/**

 * @author "Luis W. Sevilla" (sevilla_lui@gva.es)

 */

public interface IDatum {

    /**

     * Semieje Mayor del elipsoide.

     * @return

     */

    public double getESemiMajorAxis();

    /**

     * Aplanamiento del elipsoide.

     * @return

     */

    public double getEIFlattening();

}

/*

 * Cresques Mapping Suite. Graphic Library for constructing mapping applications.

 *

 * Copyright (C) 2004-5.

 *

 * This program is free software; you can redistribute it and/or

 * modify it under the terms of the GNU General Public License

 * as published by the Free Software Foundation; either version 2

 * of the License, or (at your option) any later version.

 *

 * This program is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.

 *

 * For more information, contact:

 *

 * cresques@gmail.com

 */

package org.cresques.cts;

import org.cresques.geo.ViewPortData;

import java.awt.Color;

import java.awt.Graphics2D;

import java.awt.geom.Point2D;

import java.awt.geom.Rectangle2D;

import org.gvsig.tools.lang.Cloneable;

/**

 *

 * @author "Luis W. Sevilla" (sevilla_lui@gva.es)

 */

public interface IProjection extends Cloneable{

    public IDatum getDatum();

    public Point2D createPoint(double x, double y);

    // TODO Quitar si no son necesarias.

    public String getAbrev();

    /**

     * Devuelve getAbrev() mas los parametros de transformacion si los hay

     * ej.: (EPSG:23030:proj@+proj...@...)

     * 

     * @return getAbrev() o getAbrev()+parametros

     */

    public String getFullCode();

    public void drawGrid(Graphics2D g, ViewPortData vp);

    public void setGridColor(Color c);

    public Color getGridColor();

    /**

     * Crea un ICoordTrans para transformar coordenadas

     * desde el IProjection actual al dest.

     * @param dest

     * @return

     */

    public ICoordTrans getCT(IProjection dest);

    public Point2D toGeo(Point2D pt);

    public Point2D fromGeo(Point2D gPt, Point2D mPt);

    public boolean isProjected();

    public double getScale(double minX, double maxX, double width, double dpi);

    public Rectangle2D getExtent(Rectangle2D extent,double scale,double wImage,double hImage,double mapUnits,double distanceUnits,double dpi);

}

package org.cresques.cts;

public interface UTM extends IProjection {

}

/*

 * Cresques Mapping Suite. Graphic Library for constructing mapping applications.

 *

 * Copyright (C) 2004-5.

 *

 * This program is free software; you can redistribute it and/or

 * modify it under the terms of the GNU General Public License

 * as published by the Free Software Foundation; either version 2

 * of the License, or (at your option) any later version.

 *

 * This program is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.

 *

 * For more information, contact:

 *

 * cresques@gmail.com

 */

package org.cresques.cts;

import java.awt.geom.Point2D;

import java.awt.geom.Rectangle2D;

/**

 * @author "Luis W. Sevilla" (sevilla_lui@gva.es)

 */

public interface ICoordTrans {

    public IProjection getPOrig();

    public IProjection getPDest();

    public Point2D convert(Point2D ptOrig, Point2D ptDest);

    public Rectangle2D convert(Rectangle2D rect);

    public ICoordTrans getInverted();

}

package org.cresques.cts;

public interface ICRSFactory {

	public boolean doesRigurousTransformations();

    /**

     * Devuelve una proyeccion a partir de una cadena.

     * @param name abreviatura de la proyecccion (i.e. EPSG:23030)

     * @return Proyeccion si existe

     */

    public IProjection get(String name);

}

/*

 * Cresques Mapping Suite. Graphic Library for constructing mapping applications.

 *

 * Copyright (C) 2004-5.

 *

 * This program is free software; you can redistribute it and/or

 * modify it under the terms of the GNU General Public License

 * as published by the Free Software Foundation; either version 2

 * of the License, or (at your option) any later version.

 *

 * This program is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.

 *

 * For more information, contact:

 *

 * cresques@gmail.com

 */

package org.cresques.cts;

import java.awt.geom.Point2D;

/**

 * Operaciones relacionadas con las proyecciones y sistemas

 * de coordenadas.

 *

 * cmartinez: Esta clase no deber?a formar parte de una API, pero

 * se deja hasta que se aborde el refactoring de libProjection.

 *

 * @author Luis W. Sevilla (sevilla_lui@gva.es)

 */

public class GeoCalc {

    IProjection proj;

    /**

     *

     * @param proj

     */

    public GeoCalc(IProjection proj) {

        this.proj = proj;

    }

    /* (non-Javadoc)

	 * @see org.cresques.impl.cts.GeoCalc#distanceGeo(java.awt.geom.Point2D, java.awt.geom.Point2D)

	 */

    public double distanceGeo(Point2D pt1, Point2D pt2) {

        double R2 = Math.pow(proj.getDatum().getESemiMajorAxis(), 2);

        double dLat = Math.toRadians(pt2.getY() - pt1.getY());

        double dLong = Math.toRadians(pt2.getX() - pt1.getX());

        double alfa = Math.toRadians(pt1.getY());

        double alfa2 = Math.toRadians(pt2.getY());

        if (Math.abs(alfa2) < Math.abs(alfa)) {

            alfa = alfa2;

        }

        double ds2 = (R2 * dLat * dLat) +

                     (R2 * Math.cos(alfa) * Math.cos(alfa) * dLong * dLong);

        return Math.sqrt(ds2);

    }

    /* (non-Javadoc)

	 * @see org.cresques.impl.cts.GeoCalc#distanceEli(java.awt.geom.Point2D, java.awt.geom.Point2D)

	 */

    public double distanceEli(Point2D pt1, Point2D pt2) {

        double lat1 = Math.toRadians(pt1.getY());

        double lon1 = -Math.toRadians(pt1.getX());

        double lat2 = Math.toRadians(pt2.getY());

        double lon2 = -Math.toRadians(pt2.getX());

        double F = (lat1 + lat2) / 2D;

        double G = (lat1 - lat2) / 2D;

        double L = (lon1 - lon2) / 2D;

        double sing = Math.sin(G);

        double cosl = Math.cos(L);

        double cosf = Math.cos(F);

        double sinl = Math.sin(L);

        double sinf = Math.sin(F);

        double cosg = Math.cos(G);

        double flat = 1D / proj.getDatum().getEIFlattening();

        double S = (sing * sing * cosl * cosl) + (cosf * cosf * sinl * sinl);

        double C = (cosg * cosg * cosl * cosl) + (sinf * sinf * sinl * sinl);

        double W = Math.atan2(Math.sqrt(S), Math.sqrt(C));

        double R = Math.sqrt((S * C)) / W;

        double H1 = ((3D * R) - 1D) / (2D * C);

        double H2 = ((3D * R) + 1D) / (2D * S);

        double D = 2D * W * proj.getDatum().getESemiMajorAxis();

        return (D * ((1D + (flat * H1 * sinf * sinf * cosg * cosg)) -

               (flat * H2 * cosf * cosf * sing * sing)));

    }

    /**

     * Algrothims from Geocentric Datum of Australia Technical Manual

     *

     * http://www.anzlic.org.au/icsm/gdatum/chapter4.html

     *

     * This page last updated 11 May 1999

     *

     * Computations on the Ellipsoid

     *

     * There are a number of formulae that are available

     * to calculate accurate geodetic positions,

     * azimuths and distances on the ellipsoid.

     *

     * Vincenty's formulae (Vincenty, 1975) may be used

     * for lines ranging from a few cm to nearly 20,000 km,

     * with millimetre accuracy.

     * The formulae have been extensively tested

     * for the Australian region, by comparison with results

     * from other formulae (Rainsford, 1955 & Sodano, 1965).

     *

     * * Inverse problem: azimuth and distance from known

     *                 latitudes and longitudes

     * * Direct problem: Latitude and longitude from known

     *                 position, azimuth and distance.

     * * Sample data

     * * Excel spreadsheet

     *

     * Vincenty's Inverse formulae

     * Given: latitude and longitude of two points

     *                 (phi1, lembda1 and phi2, lembda2),

     * Calculate: the ellipsoidal distance (s) and

     * forward and reverse azimuths between the points (alpha12, alpha21).

     */

    /* (non-Javadoc)

	 * @see org.cresques.impl.cts.GeoCalc#distanceVincenty(java.awt.geom.Point2D, java.awt.geom.Point2D)

	 */

    public double distanceVincenty(Point2D pt1, Point2D pt2) {

        return distanceAzimutVincenty(pt1, pt2).dist;

    }

    /**

     * Returns the distance between two geographic points on the ellipsoid

     *        and the forward and reverse azimuths between these points.

     *       lats, longs and azimuths are in decimal degrees, distance in metres

     *        Returns ( s, alpha12,  alpha21 ) as a tuple

     * @param pt1

     * @param pt2

     * @return

     */

    protected GeoData distanceAzimutVincenty(Point2D pt1, Point2D pt2) {

        GeoData gd = new GeoData(0, 0);

        double f = 1D / proj.getDatum().getEIFlattening();

        double a = proj.getDatum().getESemiMajorAxis();

        double phi1 = pt1.getY();

        double lembda1 = pt1.getX();

        double phi2 = pt2.getY();

        double lembda2 = pt2.getX();

        if ((Math.abs(phi2 - phi1) < 1e-8) &&

                (Math.abs(lembda2 - lembda1) < 1e-8)) {

            return gd;

        }

        double piD4 = Math.atan(1.0);

        double two_pi = piD4 * 8.0;

        phi1 = (phi1 * piD4) / 45.0;

        lembda1 = (lembda1 * piD4) / 45.0; // unfortunately lambda is a key word!

        phi2 = (phi2 * piD4) / 45.0;

        lembda2 = (lembda2 * piD4) / 45.0;

        double b = a * (1.0 - f);

        double TanU1 = (1 - f) * Math.tan(phi1);

        double TanU2 = (1 - f) * Math.tan(phi2);

        double U1 = Math.atan(TanU1);

        double U2 = Math.atan(TanU2);

        double lembda = lembda2 - lembda1;

        double last_lembda = -4000000.0; // an impossibe value

        double omega = lembda;

        // Iterate the following equations,

        //  until there is no significant change in lembda

        double Sin_sigma = 0;

        // Iterate the following equations,

        //  until there is no significant change in lembda

        double Cos_sigma = 0;

        double Cos2sigma_m = 0;

        double alpha = 0;

        double sigma = 0;

        double sqr_sin_sigma = 0;

        while ((last_lembda < -3000000.0) ||

                   ((lembda != 0) &&

                   (Math.abs((last_lembda - lembda) / lembda) > 1.0e-9))) {

            sqr_sin_sigma = Math.pow(Math.cos(U2) * Math.sin(lembda), 2) +

                            Math.pow(((Math.cos(U1) * Math.sin(U2)) -

                                     (Math.sin(U1) * Math.cos(U2) * Math.cos(lembda))),

                                     2);

            Sin_sigma = Math.sqrt(sqr_sin_sigma);

            Cos_sigma = (Math.sin(U1) * Math.sin(U2)) +

                        (Math.cos(U1) * Math.cos(U2) * Math.cos(lembda));

            sigma = Math.atan2(Sin_sigma, Cos_sigma);

            double Sin_alpha = (Math.cos(U1) * Math.cos(U2) * Math.sin(lembda)) / Math.sin(sigma);

            alpha = Math.asin(Sin_alpha);

            Cos2sigma_m = Math.cos(sigma) -

                          ((2 * Math.sin(U1) * Math.sin(U2)) / Math.pow(Math.cos(alpha),

                                                                        2));

            double C = (f / 16) * Math.pow(Math.cos(alpha), 2) * (4 +

                       (f * (4 - (3 * Math.pow(Math.cos(alpha), 2)))));

            last_lembda = lembda;

            lembda = omega +

                     ((1 - C) * f * Math.sin(alpha) * (sigma +

                     (C * Math.sin(sigma) * (Cos2sigma_m +

                     (C * Math.cos(sigma) * (-1 +

                     (2 * Math.pow(Cos2sigma_m, 2))))))));

        }

        double u2 = (Math.pow(Math.cos(alpha), 2) * ((a * a) - (b * b))) / (b * b);

        double A = 1 +

                   ((u2 / 16384) * (4096 +

                   (u2 * (-768 + (u2 * (320 - (175 * u2)))))));

        double B = (u2 / 1024) * (256 +

                   (u2 * (-128 + (u2 * (74 - (47 * u2))))));

        double delta_sigma = B * Sin_sigma * (Cos2sigma_m +

                             ((B / 4) * ((Cos_sigma * (-1 +

                             (2 * Math.pow(Cos2sigma_m, 2)))) -

                             ((B / 6) * Cos2sigma_m * (-3 +

                             (4 * sqr_sin_sigma)) * (-3 +

                             (4 * Math.pow(Cos2sigma_m, 2)))))));

        double s = b * A * (sigma - delta_sigma);

        double alpha12 = Math.atan2((Math.cos(U2) * Math.sin(lembda)),

                                    ((Math.cos(U1) * Math.sin(U2)) -

                                    (Math.sin(U1) * Math.cos(U2) * Math.cos(lembda))));

        double alpha21 = Math.atan2((Math.cos(U1) * Math.sin(lembda)),

                                    ((-Math.sin(U1) * Math.cos(U2)) +

                                    (Math.cos(U1) * Math.sin(U2) * Math.cos(lembda))));

        if (alpha12 < 0.0) {

            alpha12 = alpha12 + two_pi;

        }

        if (alpha12 > two_pi) {

            alpha12 = alpha12 - two_pi;

        }

        alpha21 = alpha21 + (two_pi / 2.0);

        if (alpha21 < 0.0) {

            alpha21 = alpha21 + two_pi;

        }

        if (alpha21 > two_pi) {

            alpha21 = alpha21 - two_pi;

        }

        alpha12 = (alpha12 * 45.0) / piD4;

        alpha21 = (alpha21 * 45.0) / piD4;

        return new GeoData(0, 0, s, alpha12, alpha21);

    }

    /**

     * Vincenty's Direct formulae

     * Given: latitude and longitude of a point (phi1, lembda1) and

     * the geodetic azimuth (alpha12)

     * and ellipsoidal distance in metres (s) to a second point,

     *

     * Calculate: the latitude and longitude of the second point (phi2, lembda2)

     * and the reverse azimuth (alpha21).

     */

	/**

	 * Returns the lat and long of projected point and reverse azimuth

	 *        given a reference point and a distance and azimuth to project.

	 *  lats, longs and azimuths are passed in decimal degrees.

	 * Returns ( phi2,  lambda2,  alpha21 ) as a tuple

	 * @param pt

	 * @param azimut

	 * @param dist

	 * @return

	 */

    protected GeoData getPointVincenty(Point2D pt, double azimut, double dist) {

        GeoData ret = new GeoData(0, 0);

        double f = 1D / proj.getDatum().getEIFlattening();

        double a = proj.getDatum().getESemiMajorAxis();

        double phi1 = pt.getY();

        double lembda1 = pt.getX();

        double alpha12 = azimut;

        double s = dist;

        double piD4 = Math.atan(1.0);

        double two_pi = piD4 * 8.0;

        phi1 = (phi1 * piD4) / 45.0;

        lembda1 = (lembda1 * piD4) / 45.0;

        alpha12 = (alpha12 * piD4) / 45.0;

        if (alpha12 < 0.0) {

            alpha12 = alpha12 + two_pi;

        }

        if (alpha12 > two_pi) {

            alpha12 = alpha12 - two_pi;

        }

        double b = a * (1.0 - f);

        double TanU1 = (1 - f) * Math.tan(phi1);

        double U1 = Math.atan(TanU1);

        double sigma1 = Math.atan2(TanU1, Math.cos(alpha12));

        double Sinalpha = Math.cos(U1) * Math.sin(alpha12);

        double cosalpha_sq = 1.0 - (Sinalpha * Sinalpha);

        double u2 = (cosalpha_sq * ((a * a) - (b * b))) / (b * b);

        double A = 1.0 +

                   ((u2 / 16384) * (4096 +

                   (u2 * (-768 + (u2 * (320 - (175 * u2)))))));

        double B = (u2 / 1024) * (256 +

                   (u2 * (-128 + (u2 * (74 - (47 * u2))))));

        // Starting with the approximation

        double sigma = (s / (b * A));

        double last_sigma = (2.0 * sigma) + 2.0; // something impossible

        // Iterate the following three equations

        //  until there is no significant change in sigma

        // two_sigma_m , delta_sigma

        double two_sigma_m = 0;

        while (Math.abs((last_sigma - sigma) / sigma) > 1.0e-9) {

            two_sigma_m = (2 * sigma1) + sigma;

            double delta_sigma = B * Math.sin(sigma) * (Math.cos(two_sigma_m) +

                                 ((B / 4) * (Math.cos(sigma) * ((-1 +

                                 (2 * Math.pow(Math.cos(two_sigma_m), 2))) -

                                 ((B / 6) * Math.cos(two_sigma_m) * (-3 +

                                 (4 * Math.pow(Math.sin(sigma), 2))) * (-3 +

                                 (4 * Math.pow(Math.cos(two_sigma_m), 2))))))));

            last_sigma = sigma;

            sigma = (s / (b * A)) + delta_sigma;

        }

        double phi2 = Math.atan2(((Math.sin(U1) * Math.cos(sigma)) +

                                 (Math.cos(U1) * Math.sin(sigma) * Math.cos(alpha12))),

                                 ((1 - f) * Math.sqrt(Math.pow(Sinalpha, 2) +

                                                      Math.pow((Math.sin(U1) * Math.sin(sigma)) -

                                                               (Math.cos(U1) * Math.cos(sigma) * Math.cos(alpha12)),

                                                               2))));

        double lembda = Math.atan2((Math.sin(sigma) * Math.sin(alpha12)),

                                   ((Math.cos(U1) * Math.cos(sigma)) -

                                   (Math.sin(U1) * Math.sin(sigma) * Math.cos(alpha12))));

        double C = (f / 16) * cosalpha_sq * (4 + (f * (4 - (3 * cosalpha_sq))));

        double omega = lembda -

                       ((1 - C) * f * Sinalpha * (sigma +

                       (C * Math.sin(sigma) * (Math.cos(two_sigma_m) +

                       (C * Math.cos(sigma) * (-1 +

                       (2 * Math.pow(Math.cos(two_sigma_m), 2))))))));

        double lembda2 = lembda1 + omega;

        double alpha21 = Math.atan2(Sinalpha,

                                    ((-Math.sin(U1) * Math.sin(sigma)) +

                                    (Math.cos(U1) * Math.cos(sigma) * Math.cos(alpha12))));

        alpha21 = alpha21 + (two_pi / 2.0);

        if (alpha21 < 0.0) {

            alpha21 = alpha21 + two_pi;

        }

        if (alpha21 > two_pi) {

            alpha21 = alpha21 - two_pi;

        }

        phi2 = (phi2 * 45.0) / piD4;

        lembda2 = (lembda2 * 45.0) / piD4;

        alpha21 = (alpha21 * 45.0) / piD4;

        ret.pt = new Point2D.Double(lembda2, phi2);

        ret.azimut = alpha21;

        return ret;

    }

    /* (non-Javadoc)

	 * @see org.cresques.impl.cts.GeoCalc#surfaceSphere(java.awt.geom.Point2D, java.awt.geom.Point2D, java.awt.geom.Point2D)

	 */

    public double surfaceSphere(Point2D pt1, Point2D pt2, Point2D pt3) {

        double sup = -1;