package org.mbari.gis.util;

import java.awt.geom.*;

/**
 * <p>Static methods useful for GIS Applications</p><hr><br>
 *
 * <p>MBARI provides this documentation and code &quot;as is&quot;, with no
 * warranty, express or implied, of its quality or consistency. It is provided
 * without support and without obligation on the part of the Monterey Bay
 * Aquarium Research Institute to assist in its use, correction, modification,
 * or enhancement. This information should not be published or distributed to
 * third parties without specific written permission from MBARI.</p><br>
 *
 * Copyright 2002 MBARI.<br>
 * MBARI Proprietary Information. All rights reserved.<br><hr><br>
 *
 * $Log: GisUtil.java,v $
 * Revision 1.2  2002/01/30 00:23:56  brian
 * Added comments to all classes and methods. Tested the accuracy of AsciiRasterGridReader and AsciiRasterGridWriter. Also checked the the GridUtil was populating the spatial grids correctly and accurately. The data files generated matched up perfectly with the orginal text files. Also implemented an RovAltitudeFilter class to encapsulate different filtering methods.
 *<br><hr><br>
 *
 * @author  : $Author: brian $
 * @version : $Revision: 1.2 $
 */
public class GisUtil {

    private GisUtil() {
    }

    /**
     * Converts a set of Longitude and Latitude co-ordinates to UTM for the
     * WGS84 ellipsoid
     *
     * @param Lat A double value for the latitude to be converted.
     * @param Long A double value for the longitude to be converted.
     *
     * @return A java.awt.Point where getX() returns the Eastings and getY()
     * returns the Northings
     */
    public static final Point2D geoToUtm(double longitude, double latitude, int ZoneNumber){

        double a = 6378137.0d; // WGS84 ellipsoid radius
        double eccSquared = 0.00669438d; // WGS84 eccsq
        double k0 = 0.9996;

        double longOrigin;
        double eccPrimeSquared;
        double N, T, C, A, M;

        //Make sure the longitude is between -180.00 .. 179.9
        double longTemp = (longitude+180)-((int)((longitude+180)/360))*360-180; // -180.00 .. 179.9;

        double latRad = Math.toRadians(latitude);
        double longRad = Math.toRadians(longitude);
        double longOriginRad;
//	int  ZoneNumber = (int)((longTemp + 180)/6) + 1;
//
//	if ( Lat >= 56.0f && Lat < 64.0f &&
//	     longTemp >= 3.0f && longTemp < 12.0f ) {
//	    ZoneNumber = 32;
//	}
//
//	// Special zones for Svalbard
//	if( Lat >= 72.0f && Lat < 84.0f ) {
//	    if(longTemp >= 0.0f  && longTemp <  9.0f ) ZoneNumber = 31;
//	    else if( longTemp >= 9.0f  && longTemp < 21.0f ) ZoneNumber = 33;
//	    else if( longTemp >= 21.0f && longTemp < 33.0f ) ZoneNumber = 35;
//	    else if( longTemp >= 33.0f && longTemp < 42.0f ) ZoneNumber = 37;
//	}
        longOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
        longOriginRad = Math.toRadians(longOrigin);

        eccPrimeSquared = (eccSquared)/(1-eccSquared);

        N = a/Math.sqrt(1-eccSquared*Math.sin(latRad)*Math.sin(latRad));
        T = Math.tan(latRad)*Math.tan(latRad);
        C = eccPrimeSquared*Math.cos(latRad)*Math.cos(latRad);
        A = Math.cos(latRad)*(longRad-longOriginRad);

        M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 -
                5*eccSquared*eccSquared*eccSquared/256)*latRad -
               (3*eccSquared/8 + 3*eccSquared*eccSquared/32 +
                45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(2*latRad) +
               (15*eccSquared*eccSquared/256 +
                45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(4*latRad) -
               (35*eccSquared*eccSquared*eccSquared/3072)*Math.sin(6*latRad));

        double utmEasting =
                (double)(k0*N*(A+(1-T+C)*A*A*A/6.0d +
                  (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120.0d) + 500000.0d);

        double utmNorthing =
                (double)(k0*(M+N*Math.tan(latRad)*
                (A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24.0d
                 + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720.0d)));
        if(latitude < 0.0f){
            utmNorthing += 10000000.0f; //10000000 meter offset for southern hemisphere
        }

        return new Point2D.Double(utmEasting, utmNorthing);
    }
}