/*
 *    GeoTools - OpenSource mapping toolkit
 *    http://geotools.org
 *
 *   (C) 2003, 2004 Geotools Project Managment Committee (PMC)
 *   (C) 2001, Institut de Recherche pour le Développement
 *   (C) 2000, Frank Warmerdam
 *   (C) 1999, Fisheries and Oceans Canada
 *
 *    This library is free software; you can redistribute it and/or
 *    modify it under the terms of the GNU Lesser General Public
 *    License as published by the Free Software Foundation; either
 *    version 2.1 of the License, or (at your option) any later version.
 *
 *    This library 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
 *    Lesser General Public License for more details.
 *
 *    This package contains formulas from the PROJ package of USGS.
 *    USGS's work is fully acknowledged here. This derived work has
 *    been relicensed under LGPL with Frank Warmerdam's permission.
 */
package org.geotools.referencing.operation.projection;

// J2SE dependencies
import java.awt.geom.Point2D;
import java.util.Collection;

// OpenGIS dependencies
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.parameter.ParameterNotFoundException;

// Geotools dependencies
import org.geotools.measure.Latitude;
import org.geotools.resources.i18n.Errors;
import org.geotools.resources.i18n.ErrorKeys;


/**
 * Lambert Conical Conformal Projection.  Areas and shapes are deformed
 * as one moves away from standard parallels.  The angles are true in
 * a limited area.  This projection is used for the charts of North America.
 * <p>
 *
 * This implementation provides transforms for three cases of the lambert conic 
 * conformal projection:
 * <ul>
 *   <li>{@code Lambert_Conformal_Conic_1SP} (EPSG code 9801)</li>
 *   <li>{@code Lambert_Conformal_Conic_2SP} (EPSG code 9802)</li>
 *   <li>{@code Lambert_Conic_Conformal_2SP_Belgium} (EPSG code 9803)</li>
 *   <li>{@code Lambert_Conformal_Conic} - An alias for the ESRI 2SP case
 *       that includes a scale_factor parameter</li>
 * </ul>
 *
 * For the 1SP case the latitude of origin is used as the standard parallel (SP). 
 * To use 1SP with a latitude of origin different from the SP, use the 2SP
 * and set the SP1 to the single SP. The "standard_parallel_2" 
 * parameter is optional and will be given the same value as "standard_parallel_1" 
 * if not set (creating a 1 standard parallel projection). 
 * <p>
 *
 * <strong>References:</strong><ul>
 *   <li>John P. Snyder (Map Projections - A Working Manual,<br>
 *       U.S. Geological Survey Professional Paper 1395, 1987)</li>
 *   <li>"Coordinate Conversions and Transformations including Formulas",<br>
 *       EPSG Guidence Note Number 7, Version 19.</li>
 * </ul>
 *
 * @see <A HREF="http://mathworld.wolfram.com/LambertConformalConicProjection.html">Lambert conformal conic projection on MathWorld</A>
 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html">lambert_conic_conformal_1sp</A>
 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html">lambert_conic_conformal_2sp</A>
 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp_belgium.html">lambert_conic_conformal_2sp_belgium</A>
 *
 * @since 2.1
 * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/LambertConformal.java $
 * @version $Id: LambertConformal.java 25485 2007-05-11 19:12:35Z desruisseaux $
 * @author André Gosselin
 * @author Martin Desruisseaux
 * @author Rueben Schulz
 */
public abstract class LambertConformal extends MapProjection {
    /**
     * Maximum difference allowed when comparing real numbers.
     */
    private static final double EPSILON = 1E-6;
    
    /** 
     * Constant for the belgium 2SP case. This is 29.2985 seconds, given here in radians.
     */
    private static final double BELGE_A = 0.00014204313635987700;

    /**
     * Standards parallel 1 in radians, for {@link #getParameterValues} implementation.
     */
    private final double phi1;
    
    /**
     * Standards parallel 2 in radians, for {@link #getParameterValues} implementation.
     */
    final private double phi2;
    
    /**
     * Internal variables for computation.
     */
    private final double n,F,rho0;
    
    /**
     * {@code true} for Belgium 2SP.
     */
    private final boolean belgium;
    
    /**
     * Constructs a new map projection from the supplied parameters.
     *
     * @param  parameters The parameter values in standard units.
     * @throws ParameterNotFoundException if a mandatory parameter is missing.
     */
    protected LambertConformal(final ParameterValueGroup parameters)
            throws ParameterNotFoundException
    {
        this(parameters, false);
    }
    
    /**
     * Constructs a new map projection from the supplied parameters.
     *
     * @param  parameters The parameter values in standard units.
     * @param  belgium {@code true} for the Belgium 2SP case.
     * @throws ParameterNotFoundException if a mandatory parameter is missing.
     */
    LambertConformal(final ParameterValueGroup parameters, final boolean belgium) 
            throws ParameterNotFoundException
    {
        //Fetch parameters 
        super(parameters);
        final Collection expected = getParameterDescriptors().descriptors();
        final boolean sp2 = expected.contains(AbstractProvider.STANDARD_PARALLEL_2);
        this.belgium = belgium;
        if (sp2) {
            double phi2;
            phi1 = doubleValue(expected, AbstractProvider.STANDARD_PARALLEL_1, parameters);
            ensureLatitudeInRange(       AbstractProvider.STANDARD_PARALLEL_1, phi1, true);
            phi2 = doubleValue(expected, AbstractProvider.STANDARD_PARALLEL_2, parameters);
            if (Double.isNaN(phi2)) {
                phi2 = phi1;
            }
            this.phi2 = phi2;
            ensureLatitudeInRange(AbstractProvider.STANDARD_PARALLEL_2, phi2, true);
        } else {
            if (belgium) {
                throw new IllegalArgumentException();
            }
            // EPSG says the 1SP case uses the latitude of origin as the SP
            phi1 = phi2 = latitudeOfOrigin;
        }
        // Compute constants
        if (Math.abs(phi1 + phi2) < EPSILON) {
            throw new IllegalArgumentException(Errors.format(ErrorKeys.ANTIPODE_LATITUDES_$2,
                                               new Latitude(Math.toDegrees(phi1)),
                                               new Latitude(Math.toDegrees(phi2))));
        }
        final double  cosphi1 = Math.cos(phi1);
        final double  sinphi1 = Math.sin(phi1);
        final boolean  secant = Math.abs(phi1-phi2) > EPSILON; // Should be 'true' for 2SP case.
        if (isSpherical) {
            if (secant) {
                n = Math.log(cosphi1 / Math.cos(phi2)) /
                    Math.log(Math.tan((Math.PI/4) + 0.5*phi2) /
                             Math.tan((Math.PI/4) + 0.5*phi1));
            } else {
                n = sinphi1;
            }
            F = cosphi1 * Math.pow(Math.tan((Math.PI/4) + 0.5*phi1), n) / n;
            if (Math.abs(Math.abs(latitudeOfOrigin) - (Math.PI/2)) >= EPSILON) {
                rho0 = F * Math.pow(Math.tan((Math.PI/4) + 0.5*latitudeOfOrigin), -n);
            } else {
                rho0 = 0.0;
            }
        } else {
            final double m1 = msfn(sinphi1, cosphi1);
            final double t1 = tsfn(phi1, sinphi1);
            if (secant) {
                final double sinphi2 = Math.sin(phi2);
                final double m2 = msfn(sinphi2, Math.cos(phi2));
                final double t2 = tsfn(phi2, sinphi2);
                n = Math.log(m1/m2) / Math.log(t1/t2);
            } else {
                n = sinphi1;
            }
            F = m1 * Math.pow(t1, -n) / n;
            if (Math.abs(Math.abs(latitudeOfOrigin) - (Math.PI/2)) >= EPSILON) {
                rho0 = F * Math.pow(tsfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin)), n);
            } else {
                rho0 = 0.0;
            }
        } 
    }
    
    /**
     * {@inheritDoc}
     */
    public ParameterValueGroup getParameterValues() {
        final ParameterValueGroup values = super.getParameterValues();
        final Collection expected = getParameterDescriptors().descriptors();
        set(expected, AbstractProvider.STANDARD_PARALLEL_1, values, phi1);
        set(expected, AbstractProvider.STANDARD_PARALLEL_2, values, phi2);
        return values;
    }
    
    /**
     * Transforms the specified (<var>&lambda;</var>,<var>&phi;</var>) coordinates
     * (units in radians) and stores the result in {@code ptDst} (linear distance
     * on a unit sphere).
     */
    protected Point2D transformNormalized(double x, double y, Point2D ptDst) 
            throws ProjectionException 
    {
        double rho;
        //Snyder p. 108
        if (Math.abs(Math.abs(y) - (Math.PI/2)) < EPSILON) {
            if (y*n <= 0) {
                throw new ProjectionException(Errors.format(ErrorKeys.POLE_PROJECTION_$1,
                                              new Latitude(Math.toDegrees(y))));
            } else {
                rho = 0;
            }
        } else if (isSpherical) {
            rho = F * Math.pow(Math.tan((Math.PI/4) + 0.5*y), -n);
        } else {
            rho = F * Math.pow(tsfn(y, Math.sin(y)), n);
        }
        
        x *= n;
        if (belgium) {
            x -= BELGE_A;
        }
        y = rho0 - rho * Math.cos(x);
        x =        rho * Math.sin(x);
        
        if (ptDst != null) {
            ptDst.setLocation(x,y);
            return ptDst;
        }
        return new Point2D.Double(x,y);
    }
     
    /**
     * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
     * and stores the result in {@code ptDst}.
     */
    protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst) 
            throws ProjectionException 
    {
        double theta;
        y = rho0 - y;
        double rho = Math.sqrt(x*x + y*y);  // Zero when the latitude is 90 degrees.
        if (rho > EPSILON) {
            if (n < 0) {
                rho = -rho;
                x = -x;
                y = -y;
            }
            theta = Math.atan2(x, y);
            if (belgium) {
                theta += BELGE_A;
            }
            x = theta/n;
            if (isSpherical) {
                y = 2.0 * Math.atan(Math.pow(F/rho, 1.0/n)) - (Math.PI/2);
            } else {
                y = cphi2(Math.pow(rho/F, 1.0/n));
            }
        } else {
            x = 0.0;
            y = n < 0 ? -(Math.PI/2) : (Math.PI/2);
        }
        if (ptDst != null) {
            ptDst.setLocation(x,y);
            return ptDst;
        }
        return new Point2D.Double(x,y);
    }
    
    /**
     * Returns a hash value for this projection.
     */
    public int hashCode() {
        /*
         * This code should be computed fast. Consequently, we do not use all fields
         * in this object.  Two {@code LambertConformal} objects with different
         * {@link #phi1} and {@link #phi2} should compute a F value different enough.
         */
        final long code = Double.doubleToLongBits(F);
        return ((int)code ^ (int)(code >>> 32)) + 37*super.hashCode();
    }
    
    /**
     * Compares the specified object with this map projection for equality.
     */
    public boolean equals(final Object object) {
        if (object == this) {
            // Slight optimization
            return true;
        }
        if (super.equals(object)) {
            final LambertConformal that = (LambertConformal) object;
            return (this.belgium == that.belgium) &&
                   equals(this.n,      that.n)    &&
                   equals(this.F,      that.F)    &&
                   equals(this.rho0,   that.rho0) &&
                   equals(this.phi1,   that.phi1) &&
                   equals(this.phi2,   that.phi2);
        }
        return false;
    }
}
