//SPCS83 - Converted from FORTRAN program by Edward E. Carlson
//Purpose:  To convert NAD 83 geodetic positions to NAD 83 state plane coordinates (meters)
//		Only handles a few states using Lambert Conformal Conic projection
// Latitude positive north, in radian measure.
// ER is equatorial radius of the ellipsoid (= major semiaxis).
// RF is reciprocal of flattening of the ellipsoid.
// FIS, FIN, FIB are respecitvely the latitudes of the south
//   standard parallel, the north standard parallel, and the
//   southernmost parallel.
// ESQ is the square of first eccentricity of the ellipsoid.
// ecc is first eccentricity.
// SINFO = SIN(FO), where FO is the central parallel.
// RB is mapping radius at the southernmost latitude.
// K is mapping radius at the equator.
// NB is false northing for the southernmost parallel.
// KO is scale factor at the central parallel.
// NO is northing of intersection of central meridian and parallel.
// G is a constant for computing chord-to-arc corrections.
// CONV is convergence.
// KP is point scale factor.
// CM is the central meridian of the projection zone.
// EO is false easting value at the central meridian.
// NB is false northing for the southernmost parallel of the projection, usually zero.

var deg2rad = Math.PI / 180;
var meter2feet = 1.0 / 0.3048; //3.28083989501;

var lccNorth = 0;	// required by calcLCCSP83 to return results
var lccEast = 0;	// required by calcLCCSP83 to return results

function splitdms(a) {
// Takes string in form "dd mm ss.s", returns decimal degrees
// No error checking
// Requires AllTrim function from "StringFunctions.js"
  a = AllTrim(a);
  var d = 0;
  var m = 0;
  var s = 0;
  var is_neg = false;
  if (a.charAt(0) == "-") {
    is_neg = true;
    a = AllTrim(a.substring(1,a.length));
  }
  var i = a.indexOf(" ");
  if (i > 0) {
    d = Number(a.substring(0,i+1));
    a = AllTrim(a.substring(i,a.length));
    var i = a.indexOf(" ");
    if (i > 0) {
      m = Number(a.substring(0,i+1));
      s = Number(AllTrim(a.substring(i,a.length)));
    } else {
      m = Number(a);
    }
  } else {
    d = Number(a);
  }
  if (is_neg == true) {
    return 0-dms2deg(d,m,s);  
  } else {
    return dms2deg(d,m,s);  
  }
}

function dms2deg(d,m,s) {
  return d + (m / 60.0) + (s / 3600.0);
}

function FormatNumber0(n,d,c) {
  s = FormatNumber(n,d,c);
  if (s.indexOf('.') == 0) s = "0" + s;
  return s;
}

function deg2dms(dd,sprec) {
// Requires FormatNumber function from "StringFunctions.js"
  if (!((sprec >= 0) && (sprec <= 8))) sprec = 2;	
  var d = Math.floor(dd);
  var dm = (dd-d) * 60;
  var m = Math.floor(dm);
  var s = FormatNumber0((dm-m) * 60, sprec , "");
  // handle rounding errors
  if (s == 60) {
	m = m + 1;
	s = FormatNumber0(0, sprec , "");
  }
  if (m == 60) {
	d = d + 1;
	m = "0";
  }
  if (m < 10) m = "0" + String(m);
  if (s < 10) s = "0" + String(s);
  return d + " " + m + " " + s;
}

function calcLCCq(Ecc,S) {
  var w = Ecc * S
  return (Math.log((1+S)/(1-S)) - Ecc * Math.log((1+w)/(1-w))) / 2.0;
}

function calcLccConsts(zone) {
// Precompute Lambert Conformal Conic values for given zone
	switch (zone) {
	case "CUS": // Conterminous US
		lccCentralMeridian = 96.0;
		lccFalseEasting = 0.0;
		lccFalseNorthing = 0.0;
		lccLatSouth = 33.0;
		lccLatNorth = 45.0;
		lccLatSmost = 39.0;
		break;
	case "EUS": // Same as Conterminous US except for CM
		lccCentralMeridian = 79.0;
		lccFalseEasting = 0.0;
		lccFalseNorthing = 0.0;
		lccLatSouth = 33.0;
		lccLatNorth = 45.0;
		lccLatSmost = 39.0;
		break;
	case "NC":
		lccCentralMeridian = 79.0;
		lccFalseEasting = 609601.22;
		lccFalseNorthing = 0.0;
		lccLatSouth = dms2deg(34,20,0);
		lccLatNorth = dms2deg(36,10,0);
		lccLatSmost = 33.75;
		break;
	case "PAN":	//PA NORTH
		lccCentralMeridian = dms2deg(77,45,0);
		lccFalseEasting = 600000.0;
		lccFalseNorthing = 0.0;
		lccLatSouth = dms2deg(40,53,0);
		lccLatNorth = dms2deg(41,57,0);
		lccLatSmost = dms2deg(40,10,0);
		break;
	case "PAS":	//PA SOUTH
		lccCentralMeridian = dms2deg(77,45,0);
		lccFalseEasting = 600000.0;
		lccFalseNorthing = 0.0;
		lccLatSouth = dms2deg(39,56,0);
		lccLatNorth = dms2deg(40,58,0);
		lccLatSmost = dms2deg(39,20,0);
		break;
	case "TN":
		lccCentralMeridian = 86.0;
		lccFalseEasting = 600000.0;
		lccFalseNorthing = 0.0;
		lccLatSouth = dms2deg(35,15,0);
		lccLatNorth = dms2deg(36,25,0);
		lccLatSmost = dms2deg(34,20,0);
		break;
	case "VAN":	//VA NORTH
		lccCentralMeridian = 78.5;
		lccFalseEasting = 3500000.0;
		lccFalseNorthing = 2000000.0;
		lccLatSouth = dms2deg(38,2,0);
		lccLatNorth = dms2deg(39,12,0);
		lccLatSmost = dms2deg(37,40,0);
		break;
	case "VAS":	//VA SOUTH
		lccCentralMeridian = 78.5;
		lccFalseEasting = 3500000.0;
		lccFalseNorthing = 1000000.0;
		lccLatSouth = dms2deg(36,46,0);
		lccLatNorth = dms2deg(37,58,0);
		lccLatSmost = dms2deg(36,20,0);
		break;
	case "VALAM":	//VA LAMBERT
		lccCentralMeridian = 79.5;
		lccFalseEasting = 0.0;
		lccFalseNorthing = 0.0;
		lccLatSouth = dms2deg(37,0,0);
		lccLatNorth = dms2deg(39,30,0);
		lccLatSmost = dms2deg(36,0,0);
		break;	
	case "UTM17":	//UTM ZONE 17
		lccCentralMeridian = 81.0;
		lccFalseEasting = 500000.0;
		lccFalseNorthing = 0.0;
		lccLatSouth = dms2deg(36,46,0);
		lccLatNorth = dms2deg(37,58,0);
		lccLatSmost = dms2deg(36,20,0);
		break;	
	default:
	} // end switch (zone)  

	var ER = 6378137.0;	//NAD83
	var RF = 298.257222101;	//NAD83

	var F = 1.0 / RF;
	var ESQ = F+F-(F*F);
	var ecc = Math.sqrt(ESQ);
	var FIS = lccLatSouth * deg2rad 
	var FIN = lccLatNorth * deg2rad 
	var FIB = lccLatSmost * deg2rad 
	var SINFS =Math.sin(FIS);
	var COSFS =Math.cos(FIS);
	var SINFN =Math.sin(FIN);
	var COSFN =Math.cos(FIN);
	var SINFB =Math.sin(FIB);

	var QS = calcLCCq(ecc,SINFS);
	var QN = calcLCCq(ecc,SINFN);
	var QB = calcLCCq(ecc,SINFB);
	var W1 = Math.sqrt(1.0 - (ESQ * (SINFS*SINFS)));
	var W2 = Math.sqrt(1.0 - (ESQ * (SINFN*SINFN)));
	var W4 = W2*COSFS/(W1*COSFN);
	var SINFO = Math.log(W2*COSFS/(W1*COSFN))/(QN-QS);
	var K = ER * COSFS * Math.exp(QS*SINFO)/(W1*SINFO);
	var RB = K / Math.exp(QB*SINFO);
	var QO = calcLCCq(ecc,SINFO);
	var RO = K / Math.exp(QO*SINFO);
	var SINFO2 = SINFO*SINFO;
	var COSFO = Math.sqrt(1.0 - SINFO2);
	var KO = Math.sqrt(1.0 - (ESQ * SINFO2)) * (SINFO/COSFO) * RO / ER;
	var NO = RB + lccFalseNorthing  - RO;
	var ERKO = ER*KO;
	var W3 = 1-(ESQ*(SINFO2));
	var G = (W3*W3) / (2 * (ERKO*ERKO)) * (1-ESQ);

	var results = new Array();
	results["cm"] = lccCentralMeridian * deg2rad;
	results["nb"] = lccFalseNorthing;
	results["eo"] = lccFalseEasting;
	results["er"] = ER;
	results["rf"] = RF;
	results["fis"] = FIS;
	results["fin"] = FIN;
	results["fib"] = FIB;
	results["esq"] = ESQ;
	results["ecc"] = ecc;
	results["sinfo"] = SINFO;
	results["rb"] = RB;
	results["k"] = K;
	results["ko"] = KO;
	results["no"] = NO;
	results["g"] = G;
	return results;
} // end function calcLccConsts(zone)

function calcLCCSP83(zone,lat,lng) {
// Given NAD83 LCC state plane code, and latitude and longitude in decimal degrees (longitude positive west)
// Returns northing and easting in meters
	var lc = calcLccConsts(zone);
	var FI = lat * deg2rad;
	var LAM = lng * deg2rad;
	var CM = lc["cm"];
	var NB = lc["nb"];
	var EO = lc["eo"];
	var ER = lc["er"];
	var ESQ = lc["esq"];
	var ecc = lc["ecc"];
	var SINFO = lc["sinfo"];
	var RB = lc["rb"];
	var K = lc["k"];

	var SINLAT= Math.sin(FI);
	var COSLAT= Math.cos(FI);
	var CONV = (CM - LAM) * SINFO;

	var Q = calcLCCq(ecc,SINLAT);
	var RPT = K / Math.exp(SINFO*Q);
	lccNorth = NB + RB - (RPT * Math.cos(CONV));
	lccEast = EO + (RPT * Math.sin(CONV));
	//var WP = Math.sqrt(1.0-ESQ*SINLAT*SINLAT);
	//var KP = WP * SINFO * RPT/(ER*COSLAT);
	var results = new Array();
	results["north"] = lccNorth ;
	results["east"] = lccEast ;
	//document.writeln(lat + ", " + lng + ", " + (1.0 / 0.3048 * lccNorth) + ", " + (1.0 / 0.3048 * lccEast));
	return results;
}

function calcLccLamr1(zone,north,east) {
// Given NAD83 LCC state plane code, and northing and easting in meters
// Returns latitude and longitude in decimal degrees (longitude positive west)
	var lc = calcLccConsts(zone);
	var CM = lc["cm"];
	var NB = lc["nb"];
	var EO = lc["eo"];
	var ER = lc["er"];
	var ESQ = lc["esq"];
	var ecc = lc["ecc"];
	var SINFO = lc["sinfo"];
	var RB = lc["rb"];
	var K = lc["k"];

	var NPR = RB - north + NB;
	var EPR = east - EO;
	var GAM = Math.atan(EPR / NPR);
	var LON = CM - (GAM / SINFO);
	var RPT = Math.sqrt((NPR * NPR) + (EPR * EPR));
	var Q = Math.log(K/RPT) / SINFO;
	var TEMP = Math.exp(Q + Q);
	var SINE = (TEMP - 1.0) / (TEMP + 1.0);
	for (i = 1; i <= 3; i++) {
		F1 = (Math.log((1.0 + SINE) / (1.0 - SINE)) - ecc * Math.log((1.0 + ecc*SINE) / (1.0 - ecc*SINE)) ) / 2.0 - Q;
		F2 = (1.0 / (1.0-SINE*SINE)) - ESQ /(1.0-ESQ*SINE*SINE);
		SINE = SINE - F1/F2;
	}
	var LAT = Math.asin(SINE);
	
	//var FI = LAT;
	//var LAM = LON;
	//var SINLAT = Math.sin(FI);
	//var COSLAT = Math.cos(FI);
	//var CONV = (CM-LAM) * SINFO;
	//var Q = (LOG((1+SINLAT)/(1-SINLAT)) - ecc * LOG((1+ecc*SINLAT)/(1-ecc*SINLAT))) / 2.0;
	//var RPT =K / Math.exp(SINFO*Q);
        //var WP = SQRT(1.0 - ESQ * SINLAT*SINLAT);
	//var KP = WP * SINFO * RPT / (ER*COSLAT);

	var results = new Array();
	results["lat"] = LAT / deg2rad;
	results["lng"] = LON / deg2rad;
	//document.writeln("<br>" + results["lat"] + ", " + results["lng"])
	return results
}