Convert Albers projection values to WGS84 projection values (ruby)

Albers to WGS84 - Ruby

  # CONSUS Albers consts (EPSG: 5070)
  RE_NAD83 = 6378137.0
  E_NAD83 = 0.0818187034 #Eccentricity
  D2R = 0.01745329251 #Pi/180

  def calcPhi(i)
    sinPhi = Math.sin(i * D2R)
    return (1.0 - pow(E_NAD83, 2.0)) * ((sinPhi/(1.0 - pow((E_NAD83 * sinPhi), 2.0))) - 1.0/(2.0 * E_NAD83) * Math.log((1.0 - E_NAD83 * sinPhi)/(1.0 + E_NAD83 * sinPhi)))
  end
  
  def pow(a, b)
    a ** b
  end
  
  def albers_to_wgs84(northing, easting)

    standardParallel1 = 29.5
    standardParallel2 = 45.5
    centralMeridian = -96.0
    originLat = 23.0
    originLon = -96.0

    m1 = Math.cos(standardParallel1 * D2R)/Math.sqrt(1.0 - pow((E_NAD83 * Math.sin(standardParallel1 * D2R)), 2.0))
    m2 = Math.cos(standardParallel2 * D2R)/Math.sqrt(1.0 - pow((E_NAD83 * Math.sin(standardParallel2 * D2R)), 2.0))

    q0 = calcPhi(originLat)
    q1 = calcPhi(standardParallel1)
    q2 = calcPhi(standardParallel2)
    nc = (pow(m1, 2.0) - pow(m2, 2.0)) / (q2 - q1)
    cc = pow(m1, 2.0) + nc * q1
    rho0 = RE_NAD83 * Math.sqrt(cc - nc * q0) / nc
    rho = Math.sqrt(pow(easting, 2.0) + pow((rho0 - northing), 2.0))
    q = (cc - pow((rho * nc / RE_NAD83), 2.0)) / nc
    beta = Math.asin(q / (1.0 - Math.log((1.0 - E_NAD83) / (1.0 + E_NAD83)) * (1.0 - pow(E_NAD83, 2.0))/(2.0 * E_NAD83)))
    a = 1.0 / 3.0 * pow(E_NAD83, 2.0) + 31.0 / 180.0 * pow(E_NAD83, 4.0) + 517.0 / 5040.0 * pow(E_NAD83, 6.0)
    b = 23.0/360.0 * pow(E_NAD83, 4.0) + 251.0 / 3780.0 * pow(E_NAD83, 6.0)
    c = 761.0/45360.0 * pow(E_NAD83, 6.0)
    theta = Math.atan2(easting, (rho0 - northing))

    lat = (beta + a * Math.sin(2.0 * beta) + b * Math.sin(4.0 * beta) + c * Math.sin(6.0 * beta))/D2R
    lon = centralMeridian + (theta / D2R) / nc
    [lat, lon]
  end