var half  =  0.5
   var one   =  1.0
   var two   =  2.0
   var three =  3.0
   var four  =  4.0
   var five  =  5.0
   var six   =  6.0
   var eight =  8.0
   var nine  =  9.0
   var ten   = 10.0

   var pie  = 3.1415926536

   var one_hundred_eighty   = 180.0
   var one_hundred_eighty_3 = 183.0 // combines constants in computing UTM easting

   var five_hundred_thousand = 500000.0 // UTM easting offset in meters 
   var ten_million = 10000000.0 // UTM northing offset in meters for southern hemisphere

   var deg_to_rad = pie / one_hundred_eighty
   var rad_to_deg = one_hundred_eighty / pie

   var WGS84_a  = 6378137.0
   var WGS84_b  = 6356752.3142 // polar radius (meters)
   var WGS84_f  = one / 298.257223563
   var WGS84_e2 = 0.00669437999013 // eccentricity squared 
   var WGS84_e  = Math.sqrt(WGS84_e2)

   var WGS84_central_meridian_scale_factor = 0.9996 // for UTM coordinates
    // for the NAD27 datum use this ellipsoid -
   var NAD27_a  =  WGS84_a + 69.4
   var NAD27_f  =  WGS84_f + 0.000037264639
   var NAD27_b  = (one - NAD27_f)*NAD27_a
   var NAD27_e2 =  NAD27_f*(two-NAD27_f)
   var NAD27_e  =  Math.sqrt(NAD27_e2)

   var wgs84 = 0
   var nad27 = 1

   var geodetic = 0
   var utm = 1

   var coord_type = geodetic
   var initial_datum = wgs84

   var initial_utm_zone = 0
   var south = 0 // default is northern Hemisphere for UTM coordinates

   var WGS84_lat, WGS84_lon, WGS84_easting, WGS84_northing, WGS84_utm_zone, WGS84_X, WGS84_Y, WGS84_Z
   var NAD27_lat, NAD27_lon, NAD27_easting, NAD27_northing, NAD27_utm_zone, NAD27_X, NAD27_Y, NAD27_Z

   var loop_index = 0 // sadism-only

function setCoordinateFrameType(initial_frame) {

   if (initial_frame == 'geodetic')
       coord_type = geodetic
   else 
       coord_type = utm
}   

function setDatum(intal_datum) {

   if (intal_datum == 'wgs84')
       initial_datum = wgs84
   else 
       initial_datum = nad27
}

function getandsetTextInput() {

   if (coord_type == geodetic) {
       northlike = window.document.convert_coordinates.lat_input.value
        eastlike = window.document.convert_coordinates.lon_input.value
   }

   if (coord_type == utm) {

       northlike = window.document.convert_coordinates.northing_input.value
        eastlike = window.document.convert_coordinates.easting_input.value

       initial_utm_zone = window.document.convert_coordinates.utm_zone_input.value

       if (initial_datum == wgs84)
           WGS84_utm_zone = initial_utm_zone
       else 
           NAD27_utm_zone = initial_utm_zone
   }
}

function set_SoHsphere() {
   south = window.document.convert_coordinates.SoHsphere_input.checked
}

function checkValues() {

   if (coord_type == geodetic) {

       if ((northlike < -90.0) || (northlike > 90.0)) {
            alert("Latitude out of bounds - reenter value.")
            return false
       }

       if ((northlike == -90.0) || (northlike == 90.0)) {
            
            if (loop_index == 0) 
                alert("These coordinate systems have a discontinuity at each pole. Hence specifying 90 N or 90 S is somewhat meaningless if you expect conversion to a proper set of UTM grid coordinates.")
            else if (loop_index == 1)
                alert("I see that you insist on obtaining an answer despite these considerations. You should really abort this calculation and input a reasonable latitude.")
            else if (loop_index == 2)
                alert("I must say you are nothing if persistent. Unfortunately your strategy will fail to bear fruit. Please provide a latitude less than 90 degrees magnitude.")
            else if (loop_index == 3)
                alert("Do you REALLY want to continue with the current input values?")
            else if (loop_index == 4)
                alert("Are you absolutely certain?")
            else if (loop_index == 5)
                alert("Are you 100% totally insistent on being a pig-headed fool who does not know when to quit?")
            else if (loop_index == 6)
                alert("You should get help. I know a friend of a friend, and it's 50% covered by most health care policies.")
            else if (loop_index == 7)
                alert("This is your final warning. Quit this stupidity or suffer the consequences.")
            else 
                alert("HA HA! HA HA! HA HA! HA HA! HA HA!")

            ++loop_index;
            return checkValues()
       }

       if ((eastlike < -180.0) || (eastlike > 180.0)) {
            alert("Longitude out of bounds - reenter value.")
            return false
       }
   }

   if (coord_type == utm) {

       if ((initial_utm_zone < 0) || (initial_utm_zone > 60)) {
            alert("UTM out of bounds [1 to 60] - reenter value.")
            return false
       }

       if ((northlike < 0) || (northlike > ten_million)) {
            alert("UTM northing out of bounds - reenter value.")
            return false
       }

       if (south == 0) { // N hemisphere
           if (initial_datum == wgs84) {
               if (northlike > 9997965) { // 9997965 is max WGS84 northing in N hemisphere
                   alert("UTM northing out of bounds - reenter value.")
                   return false
               }
           }
           if (initial_datum == nad27) {
               if (northlike > 9997887) { // 9997887 is max NAD27 northing in N hemisphere
                   alert("UTM northing out of bounds - reenter value.")
                   return false
               }
           }
       } 
       else { // S hemisphere
          if (initial_datum == wgs84) {
               if (northlike < 2035) { // 2035 is min WGS84 northing in S hemisphere
                   alert("UTM northing out of bounds - reenter value.")
                   return false
               }
           }
           if (initial_datum == nad27) {
               if (northlike < 2113) { // 2113 is min NAD27 northing in S hemisphere
                   alert("UTM northing out of bounds - reenter value.")
                   return false
               }
           }
       }
            // Based on +/- 3 deg at the Equator -
       if ((eastlike < 166000) || (eastlike > 834000)) {
            alert("UTM easting out of bounds - reenter value.")
            return false
       }
   }

   return true
}

function DisplayTextOutput() {
   
   window.document.convert_coordinates.wgs84_deg.value = geodeticOutput(WGS84_lat, WGS84_lon)
   window.document.convert_coordinates.nad27_deg.value = geodeticOutput(NAD27_lat, NAD27_lon)
 
   window.document.convert_coordinates.wgs84_utm.value = utmOutput(WGS84_northing, WGS84_easting, WGS84_utm_zone)
   window.document.convert_coordinates.nad27_utm.value = utmOutput(NAD27_northing, NAD27_easting, NAD27_utm_zone)
}

function geodeticOutput(lat, lon) {

   var rounded  = Math.round(100000.0*lat)
       rounded /= 100000.0

   var txtMsg = "("
   
   if (rounded < 0.0) {
       txtMsg += -rounded
       txtMsg += " S, "
   }
   else {
       txtMsg +=  rounded
       txtMsg += " N, "
   }

   rounded = Math.round(100000.0*lon)
   rounded /= 100000.0

   if (rounded < 0.0) {
       txtMsg += -rounded
       txtMsg += " W)"
   }
   else {
       txtMsg +=  rounded
       txtMsg += " E)"
   }
   return txtMsg
}

function utmOutput(northing, easting, zone) {

   var txtMsg = "zone "
   txtMsg += zone
   txtMsg += " "
   
   var rounded  = Math.round(northing)

   txtMsg += rounded
   txtMsg += " N "
  
   rounded = Math.round(easting)

   txtMsg += rounded
   txtMsg += " E "

   return txtMsg
}

function ResetValues() {

   window.document.convert_coordinates.coord_frame[0].checked = false
   window.document.convert_coordinates.coord_frame[1].checked = false

   window.document.convert_coordinates.datum[0].checked = false
   window.document.convert_coordinates.datum[1].checked = false

   window.document.convert_coordinates.lat_input.value = ""
   window.document.convert_coordinates.lon_input.value = ""

   window.document.convert_coordinates.northing_input.value = ""
   window.document.convert_coordinates.easting_input.value = ""

   window.document.convert_coordinates.utm_zone_input.value = ""
   window.document.convert_coordinates.SoHsphere_input.checked = false
}

function get_Set_Transform_displayCoordinates() {

   getandsetTextInput()

   if (!checkValues(loop_index))
        return

   TransformCoordinates() // calculation engine 

   DisplayTextOutput()
}

function formDefaultAcmeURL(latitude, longitude) {

   var querystring = "http://mapper.acme.com/?ll="
 
   querystring += latitude 
   querystring += ","
   querystring += longitude
  
   querystring += "&z=15&t=T" 

   window.location.href = querystring
}

function getMap() {

   get_Set_Transform_displayCoordinates()

   formDefaultAcmeURL(WGS84_lat, WGS84_lon)
}

  // This function is called from individual trip reports wherein the 
  // option of linking to anywhere other than Acme.com does not exist.

function formDefaultURL(lat, lon, coord_typ, intal_datum, intal_zone, sud) { 

   var querystring = "http://mapper.acme.com/?ll="

   northlike = lat
    eastlike = lon

   coord_type = coord_typ
   initial_datum = intal_datum

   initial_utm_zone = intal_zone
   south = sud

   TransformCoordinates() // Obtain WGS84 values

   querystring += WGS84_lat
   querystring += ","
   querystring += WGS84_lon

   querystring += "&z=15&t=T"

   window.location.href = querystring
}

  function selectChartSize(site) {

    if (site == 'A') {
        window.document.map.site[0].checked = true
        window.document.map.site[1].checked = false
    }
    else if (site == 'F') {
        window.document.map.site[0].checked = false
        window.document.map.site[1].checked = true
    }
  
  }

  function selectChartSize(size) {

    if (size == 's') {
        window.document.map.size[1].checked = false
        window.document.map.size[2].checked = false
    }
    else if (size == 'm') {
        window.document.map.size[0].checked = false
        window.document.map.size[2].checked = false
    }
    else if (size == 'l') {
        window.document.map.size[0].checked = false
        window.document.map.size[1].checked = false
    }
    
  }
  
  function selectZoomParameter(zoom) {

    if (zoom == 'z15') {
        window.document.map.zoom[1].checked = false

    }
    else if (zoom == 'z12') {
        window.document.map.zoom[0].checked = false
    }
  
  }

    // This function is called from state trip report index pages.
    // Input is in the NAD27 CONUS datum - thence transformed to the WGS84 datum
    // prior to constructing URLs to Acme.com and FlashEarth.com .

    // Exceptionally, the Tennessee index pages use WGS84 input.
    
function formURL(lat, lon, coord_typ, intal_datum, intal_zone, sud) { 

   var querystring

   northlike = lat
    eastlike = lon

   coord_type = coord_typ
   initial_datum = intal_datum

   initial_utm_zone = intal_zone
   south = sud

   TransformCoordinates() // Obtain WGS84 values

   if (window.document.map.site[0].checked == true) {

       querystring = "http://mapper.acme.com/?ll="
 
       querystring += WGS84_lat
       querystring += ","
       querystring += WGS84_lon
   }

   else if (window.document.map.site[1].checked == true) {

       querystring = "http://mapper.acme.com/?ll="
 
       querystring += WGS84_lat
       querystring += ","
       querystring += WGS84_lon
   }

   else if (window.document.map.site[2].checked == true) {

       querystring = "http://www.flashearth.com/?lat="

       querystring += WGS84_lat
       querystring += "&lon="
       querystring += WGS84_lon
   }

   if (window.document.map.site[0].checked == true) {

      if (window.document.map.zoom[0].checked == true) 
          querystring += "&z=10&t=M" 
      else if (window.document.map.zoom[1].checked == true) 
          querystring += "&z=7&t=M"
   }

   else if (window.document.map.site[1].checked == true) {

      if (window.document.map.zoom[0].checked == true) 
          querystring += "&z=15&t=T" 
      else if (window.document.map.zoom[1].checked == true) 
          querystring += "&z=12&t=T"
   }

   else if (window.document.map.site[2].checked == true) {

      if (window.document.map.zoom[0].checked == true) 
          querystring += "&z=15&r=0&src=msa" 
      else if (window.document.map.zoom[1].checked == true) 
          querystring += "&z=12&r=0&src=msa"
   }

   window.location.href = querystring
}

  // Main function. In Javascript code, there are six input parameters -
  // northlike, eastlike, coord_type, initial_datum, utm_zone, south.

function TransformCoordinates() {

   if (initial_datum == wgs84) {

       if (coord_type == geodetic) {

           WGS84_lat = deg_to_rad*northlike
           WGS84_lon = deg_to_rad*eastlike
     
           geodetic_to_UTM(WGS84_lat, WGS84_lon, WGS84_a, wgs84)
       }
       else {
   
           WGS84_northing = northlike
           WGS84_easting  = eastlike

           UTM_to_geodetic(WGS84_northing, WGS84_easting, WGS84_a, wgs84)
       }
    
       geodetic_to_XYZ(WGS84_lat, WGS84_lon, WGS84_a, WGS84_e, wgs84)
       XYZ_to_Geodetic(NAD27_lat, NAD27_lon, NAD27_a, NAD27_b, NAD27_e, NAD27_f, nad27)

       geodetic_to_UTM(NAD27_lat, NAD27_lon, NAD27_a, nad27)
   }

   else {

       if (coord_type == geodetic) {

	   NAD27_lat = deg_to_rad*northlike 
           NAD27_lon = deg_to_rad*eastlike

           geodetic_to_UTM(NAD27_lat, NAD27_lon, NAD27_a, nad27)
       }
       else {

           NAD27_northing = northlike
           NAD27_easting  = eastlike

           UTM_to_geodetic(NAD27_northing, NAD27_easting, NAD27_a, nad27)
       }

       geodetic_to_XYZ(NAD27_lat, NAD27_lon, NAD27_a, NAD27_e, nad27)
       XYZ_to_Geodetic(WGS84_lat, WGS84_lon, WGS84_a, WGS84_b, WGS84_e, WGS84_f, wgs84)

       geodetic_to_UTM(WGS84_lat, WGS84_lon, WGS84_a, wgs84)
   }

   WGS84_lat *= rad_to_deg 
   WGS84_lon *= rad_to_deg

   NAD27_lat *= rad_to_deg 
   NAD27_lon *= rad_to_deg
}
  
function geodetic_to_UTM(lat, lon, a, datum) {

   var northing, easting, utm_zone

   var vsc // v * sin phi cos phi
   var voc // v * omega   cos phi

   var om2 // omega^2

   var  c2 // cos^2 phi
   var  t2 // tan^2 phi

   var term1, term2, term3, term4

   var e2, e4, e6

   var k0 = WGS84_central_meridian_scale_factor
   var s, c, t

   if (datum == wgs84) 
       e2 = WGS84_e2 
   else 
       e2 = NAD27_e2

   e4 = e2 * e2 // e^4
   e6 = e2 * e4 // e^6

   var a0 = one - 0.25 * e2 - 0.046875*e4 - ( 5.0/256.0)*e6
   var a2 =       0.375*(e2 + 0.25    *e4 + (15.0/128.0)*e6)

   var a4 =              (15.0/256.0)*(e4 +   0.75      *e6)
   var a6 =                                (35.0/3072.0)*e6

   s = Math.sin(lat) 
   c = Math.cos(lat) 
   t = Math.tan(lat)

   var v   = a / Math.sqrt(one - e2*s*s) 
   var psi = one + (e2/(one-e2))*c*c
     // longitude is in radians so must first convert to degrees prior to computing the UTM zone number -
   utm_zone = Math.floor((rad_to_deg*lon + one_hundred_eighty)/six) + 1

   if (datum == wgs84) 
       WGS84_utm_zone = utm_zone
   else 
       NAD27_utm_zone = utm_zone

     // omega is in radians  
   var omega  = lon + deg_to_rad*(one_hundred_eighty_3 - six*utm_zone) 

   vsc = v * s * c
   om2 = omega * omega

   c2 = c * c
   t2 = t * t
               // zeroth order northing 
   var m = a * (a0*lat - a2*Math.sin(two*lat) + a4*Math.sin(four*lat) - a6*Math.sin(six*lat))

   term1 = vsc *   om2         / two 
   term2 = vsc * ((om2*om2)    / 24.0)*c2   *(four*psi*psi + psi - t2) 
  
   term3 = vsc * ((om2*om2*om2)/720.0)*c2*c2*(eight*psi*psi*psi*psi*(11.0 - 24.0*t2)
                                                 - 28.0*psi*psi*psi*(one  -  six*t2)
		                                         +  psi*psi*(one  - 32.0*t2)
                                                         -  two*psi*t2 +      t2*t2)
  
   term4 = vsc * ((om2*om2*om2*om2)/40320.0)*c2*c2*c2*(1385.0 - 3111.0*t2 
                                                              +  453.0*t2*t2 - t2*t2*t2)
  
   northing = k0 * (m + term1 + term2 + term3 + term4) // the northing

   if (s < 0.0) 
       northing += ten_million // add offset if in southern hemisphere

   if (datum == wgs84) 
       WGS84_northing = northing
   else 
       NAD27_northing = northing

   voc = v * omega * c // zeroth order easting 

   term1 = voc 
   term2 = voc * (om2/      six)*c2*   (     psi -                    t2) 
  
   term3 = voc * (om2*om2/120.0)*c2*c2*(four*psi*psi*psi*(one -   six*t2)
                                      +          psi*psi*(one + eight*t2)
                                      -          two*psi             *t2
		                      +                            t2*t2)

   term4 = voc * (om2*om2*om2/5040.0)*c2*c2*c2*(61.0 - 479.0*t2
                                                     + 179.0*t2*t2 - t2*t2*t2)

   easting  = k0 * (term1 + term2 + term3 + term4) 
   easting += five_hundred_thousand; // add offset in meters

   if (datum == wgs84) 
       WGS84_easting = easting
   else 
       NAD27_easting = easting
}

function UTM_to_geodetic(northing, easting, a, datum) {

   var phi
   var temp1, temp2, temp3 
   var x, y

   var k0 = WGS84_central_meridian_scale_factor
   var omega

   var e2, e4, e6
   var e2prime, e1

   if (datum == wgs84) 
       e2 = WGS84_e2 
   else 
       e2 = NAD27_e2

   e4 = e2 * e2 // e^4
   e6 = e2 * e4 // e^6

   e2prime = e2/(one - e2) // e^2 / (1 - e^2) 
   e1 = (one - Math.sqrt(one - e2)) / (one + Math.sqrt(one - e2))
  
         // remove easting offset in meters
   x = easting - five_hundred_thousand 
         
   y = northing
   if (south == 1) // remove 10,000,000 meter offset if in southern hemisphere
       y -= ten_million

   y /= (k0*a) // normalize false northing 
           // obtain what I a perceive as a rectifying latitude - 
   phi = y / (one - 0.25*e2 - (3.0/64.0)*e4 - (5.0/256.0)*e6)
           // add terms recognized as used in converting from rectifying to geodetic latitude -
   phi += (1.5*e1       - 0.84375*e1*e1*e1   ) * Math.sin(two *phi)
   phi += (1.3125*e1*e1 - 1.71875*e1*e1*e1*e1) * Math.sin(four*phi)

   phi += ((151.0 / 96.0)*e1*e1*e1   ) * Math.sin(six  *phi)
   phi += ((1097.0/512.0)*e1*e1*e1*e1) * Math.sin(eight*phi) // the author calls this "phi1Rad"

     // ^ For internal consistency this fourth order term should not be included. 
     //   However the term's magnitude is far less than the inherent error owing to 
     //   imprecision of input constants such as a, f - even for WGS84 reference ellipsoid.

   var s, c, t
   var s2prime, c2prime
   var v, t2, curvature

   s  = Math.sin(phi)
   s2prime = Math.sqrt(one - e2*s*s) // sqrt[1 - e^2 sin^2 (phi)]

   v  = a / s2prime // N1 of code
  
   t  = Math.tan(phi)
   t2 = t * t // t2 = T1 of code

   c       = Math.cos(phi)
   c2prime = e2prime * c * c // C1 of code

   curvature = a * (one - e2) / Math.pow(s2prime, 3) 

   x /= (k0*v) // zeroth order offset from central meridian

   temp1  = half*x*x 

   temp2  = five + three*t2 + ten*c2prime - four*c2prime*c2prime - nine*e2prime
   temp2 *= x*x*x*x / 24.0  

   temp3  = 61.0 + 90.0*t2 +    298.0*c2prime + 45.0*t2*t2 
              - 252.0*e2prime - three*c2prime*c2prime
   temp3 *= x*x*x*x*x*x / 720.0  
  
   phi -= (v * t / curvature)*(temp1 - temp2 + temp3) // geodetic latitude in radians
  
   temp1  = x // compute offset in radians from UTM central meridian 

   temp2  = one + two*t2 + c2prime 
   temp2 *= x*x*x/six

   temp3  = five - two*c2prime + 28.0*t2 - three*c2prime*c2prime
               + eight*e2prime + 24.0*t2*t2
   temp3 *= x*x*x*x*x / 120.0  

   omega  = temp1 - temp2 + temp3
   omega /= c  
	
   var zone_center = six*initial_utm_zone - one_hundred_eighty_3 // in degrees
   var lon = deg_to_rad * zone_center + omega

   if (datum == wgs84) {
       WGS84_lat = phi
       WGS84_lon = lon
   }
   else {
       NAD27_lat = phi 
       NAD27_lon = lon
   }

}

function geodetic_to_XYZ(lat, lon, a, e, datum) {

   var r_curv, s_lat, temp
   var x, y, z

   s_lat  = Math.sin(lat) // sin phi

   temp   = e * s_lat
   r_curv = a / Math.sqrt(one - temp*temp) // radius of curvature in prime meridian ("N")

   temp = r_curv*Math.cos(lat) // (N + h) cos phi 

   x = temp * Math.cos(lon) 
   y = temp * Math.sin(lon)

   z = (r_curv*(one - e*e))*s_lat  

   if (datum == wgs84) {

       WGS84_X = x               
       WGS84_Y = y                 
       WGS84_Z = z

       NAD27_X = x +   8.0 
       NAD27_Y = y - 160.0 
       NAD27_Z = z - 176.0
   }

   else {
	   
       NAD27_X = x               
       NAD27_Y = y                 
       NAD27_Z = z

       WGS84_X = x -   8.0 
       WGS84_Y = y + 160.0 
       WGS84_Z = z + 176.0
   }

}

function XYZ_to_Geodetic(lat, lon, a, b, e, f, datum) {

   var e_p2, p, temp
   var theta, s, c 
   var x, y, z

   if (datum == wgs84) {
       x = WGS84_X 
       y = WGS84_Y 
       z = WGS84_Z
   }
   else {
       x = NAD27_X 
       y = NAD27_Y 
       z = NAD27_Z
   }

   lon = Math.atan2(y, x) // the geodetic longitude -

      // intermediate variables for the geodetic latitude -
   temp = one/ (one - f) 
   e_p2 = temp*temp - one // (a^2 - b^2) / b^2 = (1 - f)^{-2}  - 1

   p = Math.sqrt(x*x + y*y) 

   theta = Math.atan(z*a / (p*b))

   temp = Math.sin(theta) 
   s = temp*temp*temp
   temp = Math.cos(theta) 
   c = temp*temp*temp

             // the geodetic latitude -
   lat = Math.atan((z + e_p2*b*s) / (p - e*e*a*c)) 

   if (datum == wgs84) {
       WGS84_lat = lat 
       WGS84_lon = lon
   }
   else {
       NAD27_lat = lat 
       NAD27_lon = lon
   }

}