
if (not($?VERBOSE)) set VERBOSE = 0
if (not($?SHOW_PLOTS)) set SHOW_PLOTS = 1

macro checkfit
  if ($0 != 3)
    echo "USAGE: checkfit (file) (number)"
    break
  end

  data $1
  read -fits SkyChip.xgal FR_MINOR_MIN FR_MINOR_MAX FR_MINOR_DEL FR_MAJOR_MIN FR_MAJOR_MAX FR_MAJOR_DEL GAL_CHISQ GAL_FLUX NPIX

  delete -q Xi Yi Xv Yv Zv

  local ix iy

  $N = 0
  for ix 0 5
    for iy 0 5
      $Vx = FR_MAJOR_MIN[$2] + FR_MAJOR_DEL[$2]*$ix
      $Vy = FR_MINOR_MIN[$2] + FR_MINOR_DEL[$2]*$iy

      concat $ix Xi
      concat $iy Yi
      concat $Vx Xv
      concat $Vy Yv
 
      # echo concat GAL_CHISQ:$N[0] Zv
      concat GAL_CHISQ:$N[$2] Zv
      if ($VERBOSE)
        fprintf "%1d %1d %5.3f %5.3f %7.3f" Xi[$N] Yi[$N] Xv[$N] Yv[$N] Zv[$N]
      end
      $N ++
    end
  end

  vstat -q Zv
  set Zv = Zv * NPIX[$2] / $MIN

  if ($VERBOSE)
    fit2d Xv Yv Zv 2 -clip 5 3
  else
    fit2d Xv Yv Zv 2 -clip 5 3  -q
  end

  $DET = 1 / (4*$CX0Y2*$CX2Y0 - $CX1Y1^2)
  $Xm = ($CX1Y1*$CX0Y1 - 2*$CX0Y2*$CX1Y0) * $DET
  $Ym = ($CX1Y1*$CX1Y0 - 2*$CX2Y0*$CX0Y1) * $DET
  # echo $Xm $Ym

  applyfit2d Xv Yv Zvf
  set dZ = Zv - Zvf

  if ($SHOW_PLOTS)
    lim -n 0 Xv Yv; clear; box -ticks 1111
    vstat -q Zv
    $Range = $MAX - $MIN
    $Midpt = 0.5*($MAX + $MIN)
    zplot -pt  2 Xv Yv Zv  {$Midpt - 0.55*$Range} {$Midpt - 0.3*$Range} -c blue
    zplot -pt  7 Xv Yv Zvf {$Midpt - 0.55*$Range} {$Midpt - 0.3*$Range} -c red
    
    dot $Xm $Ym -pt 7 -sz 3 -c red

    lim -n 1 Xv Yv; clear; box -ticks 1111
    vstat -q dZ
    $Range = $MAX - $MIN
    $Midpt = 0.5*($MAX + $MIN)
    zplot -pt 10 Xv Yv dZ {$Midpt - 0.8*$Range} {$Midpt + 0.6*$Range}
    
    dot $Xm $Ym -pt 7 -sz 3 -c red
  end

  $A = $CX2Y0
  $B = $CX1Y0 + 2*$CX2Y0*$Xm + $CX1Y1*$Ym
  $C = -1

  $dXm = (sqrt($B^2 - 4*$A*$C) - $B) / (2*$A)

  $A = $CX0Y2
  $B = $CX0Y1 + 2*$CX0Y2*$Ym + $CX1Y1*$Xm
  $C = -1

  $dYm = (sqrt($B^2 - 4*$A*$C) - $B) / (2*$A)

  echo ($Xm, $Ym) +/- ($dXm, $dYm)
end
