proj = nil
sigma = 2.0
pts = []
for line in ARGF
  case line
  when /^#proj\s+/ then
    proj = $'
    next
  when /^#sigma\s+/ then
    sigma = $'.to_f
    next
  when /^#/ then next
  end
  lon, lat, x, y = line.strip.split(/[\/\s,]+/).map{|s| s.to_f}
  pxy = (`echo #{lon} #{lat} | proj #{proj}`)
  px, py = pxy.strip.split(/\s+/).map{|s| s.to_f}
  pt = {:lon => lon, :lat => lat, :x => x, :y => y, :px => px, :py => py}
  pts.push pt
end

def stat(pts, zx, zy, sig)
  ss = sx = sy = 0.0
  for pt in pts
    wt = (sig ** -2)
    ss += wt
    sx += pt[zx] * wt
    sy += pt[zy] * wt
  end
  sxoss = sx / ss
  st2 = b = 0.0
  for pt in pts
    t = (pt[zx] - sxoss) / sig
    st2 += t * t
    b += t * pt[zy] / sig
  end
  b /= st2
  a = (sy - sx * b) / ss
  siga = Math.sqrt((1.0 + sx * sx / (ss * st2)) / ss)
  sigb = Math.sqrt(1.0 / st2)
  chi2 = 0.0
  for pt in pts
    e2 = ((pt[zy] - a - b * pt[zx]) / sig) ** 2
    STDERR.puts [e2, pt].inspect if e2 > 4.0
    chi2 += e2
  end
  {:a => a, :b => b, :siga => siga, :sigb => sigb, :chi2 => chi2}
end

def repstat h
  printf("a:%-10.1f (%-10.2f) b:%13.4e (%13.3e) a/b:%-8.0f 1/b:%-6.0f chi2:%g\n",
    h[:a], h[:siga], h[:b], h[:sigb], h[:a]/h[:b], 1.0/h[:b], h[:chi2])
end

repstat stat(pts, :px, :y, sigma)
repstat stat(pts, :py, :x, sigma)