include Math

DEG2RAD = PI / 180.0
RAD2DEG = 180.0 / PI
PID4 = PI * 0.25
PID2 = PI * 0.5

def makebox xy, d
  x, y = xy
  result = []
  result.push [1, 1]
  result.push [x, 1]
  result.push [x, y]
  result.push [1, y]
  result.push [1, 1]
  result
end

def lambert fnam, conf
  slat1 = conf[:standard][1]
  slat2 = conf[:standard][3]
  slon = conf[:standard][0]
  rx = conf[:basepoint][0]
  ry = conf[:basepoint][1]
  rlon = conf[:basepoint][2]
  rlat = conf[:basepoint][3]
  d = conf[:distance][0]
  #
  e_re = 6371.0e3
  sign = (slat2 >= 0.0) ? 1.0 : -1.0
  m1 = cos(slat1 * DEG2RAD)
  t1 = tan(PID4 - sign * slat1 * DEG2RAD * 0.5)
  n = af = nil
  if slat1 == slat2
    n = cos(PID2 - sign * slat1)
    af = e_re * m1 / n / (t1 ** n)
  else
    m2 = cos(slat2 * DEG2RAD)
    t2 = tan(PID4 - sign * slat2 * DEG2RAD * 0.5)
    n = log(m1 / m2) / log(t1 / t2)
    af = e_re * m1 / n / (t1 ** n)
  end
  r0 = af * (tan(PID4 - sign * rlat * DEG2RAD * 0.5) ** n)
  x0 = rx -        r0 * sin((rlon - slon) * DEG2RAD * n) / d
  y0 = ry - sign * r0 * cos((rlon - slon) * DEG2RAD * n) / d
  xy_list = makebox(conf[:size], d)
  result = []
  for x, y in xy_list
    dx = x - x0
    dy = sign * (y - y0)
    lat = sign * (PID2 - atan(
      (sqrt(dx * dx + dy * dy) * d / af) ** (1.0 / n)
      ) * 2.0) * RAD2DEG
    lon = (atan2(dx, dy) / n) * RAD2DEG + slon
    lon -= 360.0 if lon > 180.0
    result.push [lat, lon]
  end
  open(fnam, 'w') { |fp|
    for lat, lon in result
      fp.printf "%-8.1f %-8.1f\n", lon, lat
    end
  }
end

lambert 'mf.txt', {
  :size => [721, 577],
  :basepoint => [489.0, 409.0, 140.0, 30.0],
  :distance => [5000.0, 5000.0],
  :standard => [140.0, 30.0, 140.0, 60.0]
}
lambert 'rf.txt', {
  :size => [551, 801],
  :basepoint => [320.0, 850.0, 140.0, 30.0],
  :distance => [2000.0, 2000.0],
  :standard => [140.0, 30.0, 140.0, 60.0]
}
lambert 'lf.txt', {
  :size => [1581, 1301],
  :basepoint => [1121.0, 901.0, 140.0, 30.0],
  :distance => [2000.0, 2000.0],
  :standard => [140.0, 30.0, 140.0, 60.0]
}