
macro init
  rd a /data/pikake.1/eugene/refcat.20120202/allsky.images/allsky.totals.20120705.v1.ncount_235.fits
  getcoords a ra dec
end

macro getcounts
  if ($0 != 3)
    echo "USAGE: getcounts (Dmin) (Dmax)"
    break
  end

  set mask = (dec >= $1) && (dec < $2)
  set am = a * mask
  stats -q am
end

macro t1.even
  delete -q Nd Ad
  for d -38 90 4.0
    getcounts $d {$d + 4.0}
    concat $TOTAL Nd
    $A = {2*$pi*(dsin($d+4.0) - dsin($d))*(180/$pi)^2}
    concat $A Ad
    fprintf %e $TOTAL
  end
end

# nd[2], nd[3] ~ 2.8M deg^-2
# nd @ pole ~ 500k deg^-2
# nd @ mid ~ 1.5M deg^-2

# will have ~ 1.5e9 in 1.75 deg slice at -30
# height is ~ 3.0 @ +20

# slope is ~ 0.03, from 1.75 @ -32, pole slice is ~66 - 90
# 1.75 @ -32, 4.75 @ 66, plus pole

macro t2.slope
  $D0o = -33
  $dDo = 1.8
  $mD = 0.035

  $D0 = $D0o
  $dD = $dDo

  delete -q Dx Nx Ax 

  concat $D0 Dx
  while ($D0 < 68)
    $D1 = $D0 + $dD
    getcounts $D0 $D1

    concat $TOTAL Nx
    $A = {2*$pi*(dsin($D1) - dsin($D0))*(180/$pi)^2}
    concat $A Ax

    $D0 = $D1
    $dD = $dDo + $mD*($D1 - $D0o)
    echo $dD
    concat $D0 Dx
  end

  getcounts $D0 90
  concat $TOTAL Nx
  $A = {2*$pi*(dsin(90) - dsin($D0))*(180/$pi)^2}
  concat $A Ax
end

vlist         Do -33.0 -31.2 -29.3 -27.4 -25.4 -23.3 -21.2 -19.0 -16.7 -14.3 -11.9
vlist -append Do  -9.3  -6.7  -4.0  -1.2   1.7   4.7   7.9  11.1  14.4  17.9  21.5
vlist -append Do  25.2  29.0  33.0  37.1  41.4  45.8  50.3  55.0  59.9  65.0  90.0

macro t3.fixed

  delete -q No
  for i 1 Do[]
    getcounts Do[$i-1] Do[$i]
    concat $TOTAL No
  end
  concat NAN No
end
