class Geocoder::EastingNorthing

def to_WGS84(latlng)
  latitude = latlng[0]
  longitude = latlng[1]
  a = 6_377_563.396
  b = 6_356_256.909
  eSquared = ((a * a) - (b * b)) / (a * a)
  phi = deg_to_rad(latitude)
  lambda = deg_to_rad(longitude)
  v = a / Math.sqrt(1 - eSquared * sin_pow_2(phi))
  h = 0
  x = (v + h) * Math.cos(phi) * Math.cos(lambda)
  y = (v + h) * Math.cos(phi) * Math.sin(lambda)
  z = ((1 - eSquared) * v + h) * Math.sin(phi)
  tx = 446.448
  ty = -124.157
  tz = 542.060
  s  = -0.0000204894
  rx = deg_to_rad(0.00004172222)
  ry = deg_to_rad(0.00006861111)
  rz = deg_to_rad(0.00023391666)
  xB = tx + (x * (1 + s)) + (-rx * y) + (ry * z)
  yB = ty + (rz * x) + (y * (1 + s)) + (-rx * z)
  zB = tz + (-ry * x) + (rx * y) + (z * (1 + s))
  a = 6_378_137.000
  b = 6_356_752.3141
  eSquared = ((a * a) - (b * b)) / (a * a)
  lambdaB = rad_to_deg(Math.atan(yB / xB))
  p = Math.sqrt((xB * xB) + (yB * yB))
  phiN = Math.atan(zB / (p * (1 - eSquared)))
  (1..10).each do |_i|
    v = a / Math.sqrt(1 - eSquared * sin_pow_2(phiN))
    phiN1 = Math.atan((zB + (eSquared * v * Math.sin(phiN))) / p)
    phiN = phiN1
  end
  phiB = rad_to_deg(phiN)
  [phiB, lambdaB]
end