/*
 * Copyright (C) 2007 Mark Thomas. All rights reserved.
 *
 * $Id: StatePlaneCoordinates.java 194 2007-07-21 04:27:34Z mark $
 */

package net.spatialguru.util;

import java.util.Map;
import java.util.WeakHashMap;

/**
 * This class is a port of the SPCS83 program that converts state plane
 * coordinates to latitude/longitude and vice versa; however this class
 * currently only supports converting from state plane coordinates to
 * latitude/longitude.<br/>
 * <p>The original code may be obtained
 * <a href="http://www.ngs.noaa.gov/PC_PROD/SPCS83/">here</a>
 * </p>
 */
public class StatePlaneCoordinates {

    private static final double RAD = (180d / Math.PI);
    private static final double ER  = 6378137d;
    private static final double RF  = 298.257222101d;
    private static final double F   = (1d / RF);
    private static final double ESQ = (F + F - F * F);
    private static final double E   = StrictMath.sqrt(ESQ);

    private static Map<Integer, double[]> ZONES =
      new WeakHashMap<Integer, double[]>();

    static {
        // These values are only for zone 0403 (California III)
        ZONES.put(403, new double[] {
          120.5d, 2000000.0001016d, 500000.0001016001d, 37.06666666666667d,
          38.43333333333333d, 36.5d
        });
    }

    public static double[] toLatLng(final double easting,
      final double northing, final int zone) {
        final double[] SPCC = ZONES.get(zone);
        if (SPCC == null)
            return new double[] { 0, 0 };
        final double CM = SPCC[0] / RAD;
        final double EO = SPCC[1];
        final double NB = SPCC[2];
        final double FIS = SPCC[3] / RAD;
        final double FIN = SPCC[4] / RAD;
        final double FIB = SPCC[5] / RAD;
        final double SINFS = StrictMath.sin(FIS);
        final double COSFS = StrictMath.cos(FIS);
        final double SINFN = StrictMath.sin(FIN);
        final double COSFN = StrictMath.cos(FIN);
        final double SINFB = StrictMath.sin(FIB);
        final double QS = q(E, SINFS);
        final double QN = q(E, SINFN);
        final double QB = q(E, SINFB);
        final double W1 = StrictMath.sqrt(1d - ESQ * SINFS * SINFS);
        final double W2 = StrictMath.sqrt(1d - ESQ * SINFN * SINFN);
        final double SINFO =
          (StrictMath.log(W2 * COSFS / (W1 * COSFN)) / (QN - QS));
        final double K =
          (ER * COSFS * StrictMath.exp(QS * SINFO) / (W1 * SINFO));
        final double RB = (K / StrictMath.exp(QB * SINFO));
        double NPR = RB - northing + NB;
        double EPR = easting - EO;
        double GAM = StrictMath.atan(EPR / NPR);
        double LON = CM - (GAM / SINFO);
        double RPT = StrictMath.sqrt(NPR * NPR + EPR * EPR);
        double Q = StrictMath.log(K / RPT) / SINFO;
        double TEMP = StrictMath.exp(Q + Q);
        double SINE = (TEMP - 1d) / (TEMP + 1d);
        double F1, F2;
        for (int i = 0; i < 2; i++) {
            F1 = ((StrictMath.log((1d + SINE) / (1d - SINE)) - E *
              StrictMath.log((1d + E * SINE) / (1d - E * SINE))) / 2d) - Q;
            F2 = (1d / (1d - SINE * SINE) - ESQ / (1d - ESQ * SINE * SINE));
            SINE -= (F1 / F2);
        }
        return new double[] { StrictMath.toDegrees(LON) * - 1,
          StrictMath.toDegrees(StrictMath.asin(SINE)) };
    }

    private static double q(double e, double s) {
        return ((StrictMath.log((1 + s) / (1 - s)) - e *
          StrictMath.log((1 + e * s) / (1 - e * s))) / 2);
    }

    /**
     * As an example...output should be:<br/>
     * <code>-121.2284594353128<br/>
     * 37.9075149470656</code>
     */
    public static void main(String[] a) {
        final double feet2Meters = 0.3048006096012192d;
        for (Double d : toLatLng(feet2Meters * 6351504, feet2Meters * 2153727, 403))
            System.out.println(d);
    }
}