Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/test/relphot.flatcorr.dvo
===================================================================
--- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/test/relphot.flatcorr.dvo	(revision 38308)
+++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/test/relphot.flatcorr.dvo	(revision 38309)
@@ -21,6 +21,7 @@
 if (not($?PLOT)) set PLOT = 0
 
-$USE_CELL_OFFSETS = 1
-$SET_UC_RESID = 1
+$USE_CELL_OFFSETS = 1; # turn on non-zero cell offsets
+$SET_UC_RESID = 1; # write a non-zero resid value to the ubercal table (should not have an impact)
+
 
 # set various global variables
@@ -35,5 +36,6 @@
     echo ADDSTAR_RADIUS 1.0
     echo CATMODE                 SPLIT
-    echo CATFORMAT               PS1_V2
+    echo CATFORMAT               PS1_V5
+    # echo CATFORMAT               PS1_V2 -- must be at least PS1_V3 for ubercaldist to work
     echo SKY_DEPTH               4
     echo CAMERA                  gpc1
@@ -64,4 +66,6 @@
     echo "10110 GPC1.r.XY10          dep  $zpt_nominal:r $klam_nominal:r 0.000     -     - 0.0000     0     1   0.015 0.000 1.000  0.000   0x0000 0x0000  0xe0440130 0x1003bc88"
     echo "10111 GPC1.r.XY11          dep  $zpt_nominal:r $klam_nominal:r 0.000     -     - 0.0000     0     1   0.015 0.000 1.000  0.000   0x0000 0x0000  0xe0440130 0x1003bc88"
+    echo "2011  2MASS_J              ref   0.000  0.000 0.000     -     - 0.0000     0     6   0.080 0.000 2.000  0.000   0x0000 0x0000  0x0000 0x0000"
+    echo "3001  SYNTH.g              ref   0.000  0.000 0.000     -     - 0.0000     0     1     nan 0.000 2.000  0.600   0x0000 0x0000  0x0000 0x0000"
   output stdout      
 
@@ -77,4 +81,7 @@
   $CHIP_DX = 1000
   $CHIP_DY = 1000
+  $DEFECT_FRAC = 0.03
+  $OFFSET_FRAC_UC = 0.0
+  $OFFSET_FRAC_NC = 0.5
 
   # we have two sets of images: ubercaled and not-ubercaled
@@ -149,5 +156,18 @@
   init
   mkzptfile
-  mkstars 500
+  mkstars 1000
+
+  catdir $catdir
+
+  mk2mass $fileroot $catdir
+  ck2mass addstar
+
+  mksynth $fileroot $catdir
+  cksynth addstar
+
+  # to add a valid tycho test, I need to generate a zpt table for tycho as well
+  mktycho $fileroot $catdir
+  cktycho addstar
+  break
 
   tapPLAN {((mjd_uc[]*4 + mjd_nc[]*4) + (mjd_uc[] + mjd_nc[]) + (mjd_uc[] + mjd_nc[])) / 3}
@@ -179,4 +199,5 @@
     ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot
   end
+  ck2mass relphot
 
   tapDONE
@@ -241,5 +262,5 @@
 
   skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} 
-  mextract ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001)
+  mextract ra dec mag xccd yccd dbflags photflags where (abs(time - $MJD_IMAGE) < 0.0001)
  
   match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2
@@ -252,16 +273,22 @@
   set dm = mag - stars_mag_m
 
+  set good_photom = (not(photflags & 0x02))
+  foreach field dm xccd yccd
+    subset $field\_good = $field if good_photom
+  end
+
   if ($PLOT) 
     dev -n 1
-    lim mag dm; clear; box; plot mag dm
+    lim mag dm; clear; box; plot mag dm; plot -c red -pt 7 mag dm where (photflags & 0x00000002)
+    echo "red points are known bad photometry"
   end
 
   # reproduce eddie's esoteric padding
-  set cell_xbin = max(0 , min(1 , int((xccd +  8) / ($CHIP_DX / $NCELL_X))))
-  set cell_ybin = max(0 , min(1 , int((yccd + 10) / ($CHIP_DY / $NCELL_Y))))
-  subset dm00 = dm if (cell_xbin == 0) && (cell_ybin == 0)
-  subset dm01 = dm if (cell_xbin == 0) && (cell_ybin == 1)
-  subset dm10 = dm if (cell_xbin == 1) && (cell_ybin == 0)
-  subset dm11 = dm if (cell_xbin == 1) && (cell_ybin == 1)
+  set cell_xbin = max(0 , min(1 , int((xccd_good +  8) / ($CHIP_DX / $NCELL_X))))
+  set cell_ybin = max(0 , min(1 , int((yccd_good + 10) / ($CHIP_DY / $NCELL_Y))))
+  subset dm00 = dm_good if (cell_xbin == 0) && (cell_ybin == 0)
+  subset dm01 = dm_good if (cell_xbin == 0) && (cell_ybin == 1)
+  subset dm10 = dm_good if (cell_xbin == 1) && (cell_ybin == 0)
+  subset dm11 = dm_good if (cell_xbin == 1) && (cell_ybin == 1)
 
   # uncorrected values behave like this:
@@ -298,21 +325,48 @@
 
   if ("$MODE" == "setphot_uc") 
-    vstat -q dm
-    sprintf line "setphot_uc  ZP %8.2f $FILTER %6.3f %6.3f %6.3f: %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL $MEAN $SIGMA dm[]
+    vstat -q dm_good
+    sprintf line "setphot_uc  ZP %8.2f $FILTER %6.3f %6.3f %6.3f: %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL $MEAN $SIGMA dm_good[]
     tapOK {abs($MEAN) < 0.005} "$line"
+
+    subset tmp = dm if (dbflags & 0x00008000)
+    tapOK {tmp[] == dm[]} "ubercal detections have ubercal bit set"
+
+    imextract MCAL FLAGS where (abs(time - $MJD_IMAGE) < 0.0001)
+    subset tmp = FLAGS if (FLAGS & 0x0000200)
+    tapOK {tmp[] == FLAGS[]} "ubercal images have ubercal bit set"
+    
     return
   end
 
   if ("$MODE" == "setphot_nc") 
-    vstat -q dm
-    sprintf line "setphot_nc  ZP %8.2f $FILTER %6.3f %6.3f %6.3f: %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL $MEAN $SIGMA dm[]
+    vstat -q dm_good
+    sprintf line "setphot_nc  ZP %8.2f $FILTER %6.3f %6.3f %6.3f: %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL $MEAN $SIGMA dm_good[]
     tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "$line"
+
+    subset tmp = dm if (not(dbflags & 0x00008000))
+    tapOK {tmp[] == dm[]} "non-ubercal detections do NOT have ubercal bit set"
+
+    imextract FLAGS where (abs(time - $MJD_IMAGE) < 0.0001)
+    subset tmp = FLAGS if (not(FLAGS & 0x0000200))
+    tapOK {tmp[] == FLAGS[]} "non-ubercal images do NOT have ubercal bit set"
+    
     return
   end
 
   if ("$MODE" == "relphot") 
-    vstat -q dm
-    sprintf line "relphot     ZP %8.2f $FILTER %6.3f %6.3f %6.3f: %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL $MEAN $SIGMA dm[]
+    vstat -q dm_good
+    sprintf line "relphot     ZP %8.2f $FILTER %6.3f %6.3f %6.3f: %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL $MEAN $SIGMA dm_good[]
     tapOK {abs($MEAN) < 0.005} "$line"
+
+    mextract ra dec dbflags
+    subset rs = ra  if (dbflags & 0x8000)
+    subset ds = dec if (dbflags & 0x8000)
+    avmatch rs ds 0.1 $FILTER\:flags $FILTER\:uc_dist 
+
+    subset tmp = rs if ($FILTER\:flags & 0x8)
+    tapOK {tmp[] == rs[]} "objects with ubercal detections have ubercal bit set (tmp[] vs rs[])"
+ 
+    subset tmp = rs if ($FILTER\:uc_dist < 0.1)
+    tapOK {tmp[] == rs[]} "objects with ubercal detections have ubercal dist of 0.0"
     return
   end
@@ -377,8 +431,15 @@
     $options = -mjd $MJD 
     $options = $options -radec $ra $dec 
-    $options = $options -type PS1_V2 
+
+    #$options = $options -type PS1_V2 
+    #$options = $options -type PS1_V3 
+    #$options = $options -type PS1_V4 
+    $options = $options -type PS1_V5 
+
     $options = $options -coords 
     $options = $options -photcode GPC1.$FILTER.XY$ix\$iy 
     $options = $options -no-noise 
+    $options = $options -flags
+    $options = $options -bad-psfqf-frac $DEFECT_FRAC
     $options = $options -size $CHIP_DX $CHIP_DY 
     $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} 
@@ -429,5 +490,42 @@
   end
 
-  write -f "%10.6f %10.6f %7.3f  %6.1f %6.1f %7.3f" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi
+  # create bad detections / detections to ignore or avoid, save to a file?
+
+  # create 3 types of things to ignore / avoid
+  # * badPhot (deviate the photometry)
+  # * badFlag (set a bad flag bit)
+  # * poorFlag (set a bad flag bit)
+  set badPhot  = rnd(stars_r) <  $DEFECT_FRAC
+  set badFlag  = rnd(stars_r) <  $DEFECT_FRAC
+  set poorFlag = rnd(stars_r) <  $DEFECT_FRAC
+
+  # create deviant photometry (1%)
+  # I'm setting the flag value 0x02 (EXTMODEL) if the object is deviated.  This is a simple way to pass the info
+  set stars_mi = stars_mi + 0.2*badPhot
+  set stars_fl = zero(stars_mi) + 0x00000080*badFlag + 0x00000100*poorFlag + 0x00000002*badPhot
+
+  # bad flags that we sould respect:
+  # PM_SOURCE_MODE_FAIL             = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero)
+  # PM_SOURCE_MODE_SATSTAR          = 0x00000080, ///< Source model peak is above saturation
+  # PM_SOURCE_MODE_BADPSF           = 0x00000400, ///< Failed to get good estimate of object's PSF
+  # PM_SOURCE_MODE_DEFECT           = 0x00000800, ///< Source is thought to be a defect
+  # PM_SOURCE_MODE_SATURATED        = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail)
+  # PM_SOURCE_MODE_CR_LIMIT         = 0x00002000, ///< Source has crNsigma above limit
+  # PM_SOURCE_MODE_MOMENTS_FAILURE  = 0x00008000, ///< could not measure the moments
+  # PM_SOURCE_MODE_SKY_FAILURE      = 0x00010000, ///< could not measure the local sky
+  # PM_SOURCE_MODE_SKYVAR_FAILURE   = 0x00020000, ///< could not measure the local sky variance
+  # PM_SOURCE_MODE_SIZE_SKIPPED     = 0x10000000, ///< size could not be determined
+
+  # poor flags that we should respect:
+  # PM_SOURCE_MODE_POOR             = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?)
+  # PM_SOURCE_MODE_PAIR             = 0x00000020, ///< Source fitted with a double psf
+  # PM_SOURCE_MODE_BLEND            = 0x00000100, ///< Source is a blend with other sources
+  # PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N
+  # PM_SOURCE_MODE_BLEND_FIT        = 0x00400000, ///< source was fitted as a blend
+  # PM_SOURCE_MODE_ON_SPIKE         = 0x20000000, ///< peak lands on diffraction spike
+  # PM_SOURCE_MODE_ON_GHOST         = 0x40000000, ///< peak lands on ghost or glint
+  # PM_SOURCE_MODE_OFF_CHIP         = 0x80000000, ///< peak lands off edge of chip
+
+  write -f "%10.6f %10.6f %7.3f  %6.1f %6.1f %7.3f 0x%x" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi stars_fl
 end
 
@@ -450,5 +548,174 @@
   set stars_ra  = $RA_CENTER  + $RA_RANGE  * (rnd(tmp) - 0.5) / 3600 / dcos ($DEC_CENTER)
   set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) / 3600
-  set stars_mag = 15.0 + 4.0 * rnd(tmp) 
+  set stars_mag = 10.0 + 10.0 * rnd(tmp) 
+  set stars_dmag = 0.02 + 0.1*rnd(tmp)
+end
+
+macro mk2mass
+  if ($0 != 3)
+    echo "USAGE: mk2mass (fileroot) (catdir)"
+    break
+  end
+
+  local ROOT CATDIR
+
+  $ROOT    = $1
+  $CATDIR  = $2
+
+  set stars_flags = 0x00c00000*(rnd(stars_ra) < 0.1) + 0x7*(rnd(stars_ra) < 0.8)
+  write $ROOT.2mass.dat -f "%10.6f %10.6f %6.3f %6.3f 0x%x" stars_ra stars_dec stars_mag stars_dmag stars_flags
+
+  tapEXEC addstar -ref $ROOT.2mass.dat -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass -photcode 2MASS_J
+
+  set r_2mass = stars_ra
+  set d_2mass = stars_dec
+  set m_2mass = stars_mag
+  set dm_2mass = stars_dmag
+  set fl_2mass = stars_flags
+end
+
+macro ck2mass
+  if ($0 != 2)
+    echo "USAGE: ck2mass (stage)"
+    break
+  end
+
+  if ("$1" == "addstar") 
+    skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} 
+    mextract ra dec mag xccd yccd dbflags photflags where (photcode == 2011)
+    tapOK {ra[] == r_2mass[]} "read in as many 2mass stars as we generated"
+    
+    match2d -closest ra dec r_2mass d_2mass 0.001 -index1 index1 -index2 index2
+    subset tmp = index1 if (index1 >= 0)
+    tapOK {tmp[] == ra[]} "2mass sources match those generated"
+  end
+
+  if ("$1" == "relphot")
+    subset rs = r_2mass if (fl_2mass & 0x00c00000)
+    subset ds = d_2mass if (fl_2mass & 0x00c00000)
+    avmatch rs ds 0.1 flags
+    subset tmp = rs if (flags & 0x02000000); # extended in 2mass
+    tapOK {tmp[] == rs[]} "2mass extended objects properly marked"
+    
+    subset rs = r_2mass if (fl_2mass & 0x07)
+    subset ds = d_2mass if (fl_2mass & 0x07)
+    avmatch rs ds 0.1 flags
+    subset tmp = rs if (flags & 0x08000000); # extended in 2mass
+    tapOK {tmp[] == rs[]} "2mass good objects properly marked"
+  end
+end
+
+macro mksynth
+  if ($0 != 3)
+    echo "USAGE: mksynth (fileroot) (catdir)"
+    break
+  end
+
+  local ROOT CATDIR
+
+  $ROOT    = $1
+  $CATDIR  = $2
+ 
+  subset r_synth  = stars_ra    if (stars_mag < 17)
+  subset d_synth  = stars_dec   if (stars_mag < 17)
+  subset m_synth  = stars_mag   if (stars_mag < 17)
+  subset dm_synth = stars_dmag  if (stars_mag < 17)
+  subset fl_synth = stars_flags if (stars_mag < 17)
+
+  write $ROOT.synth.dat -f "%10.6f %10.6f %6.3f %6.3f 0x%x" r_synth d_synth m_synth dm_synth fl_synth
+
+  tapEXEC addstar -ref $ROOT.synth.dat -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass -photcode SYNTH.g
+end
+
+macro cksynth
+  if ($0 != 2)
+    echo "USAGE: cksynth (stage)"
+    break
+  end
+
+  if ("$1" == "addstar") 
+    skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} 
+    mextract ra dec mag xccd yccd dbflags photflags where (photcode == 3001)
+    tapOK {ra[] == r_synth[]} "read in as many synth stars as we generated"
+    
+    match2d -closest ra dec r_synth d_synth 0.001 -index1 index1 -index2 index2
+    subset tmp = index1 if (index1 >= 0)
+    tapOK {tmp[] == ra[]} "synth sources match those generated"
+  end
+
+  if ("$1" == "relphot")
+    subset rs = r_synth if (m_synth < 13.5)
+    subset ds = d_synth if (m_synth < 13.5)
+    avmatch rs ds 0.1 g:flags
+    subset tmp = rs if (g:flags & 0x04); # synth used for photometry
+    tapOK {tmp[] == rs[]} "synth objects used for photometry properly marked"
+  end
+end
+
+macro mktycho
+  if ($0 != 3)
+    echo "USAGE: mktycho (fileroot) (catdir)"
+    break
+  end
+
+  local ROOT CATDIR
+
+  $ROOT    = $1
+  $CATDIR  = $2
+
+  subset r_tycho  = stars_ra    if (stars_mag < 12)
+  subset d_tycho  = stars_dec   if (stars_mag < 12)
+  subset m_tycho  = stars_mag   if (stars_mag < 12)
+  subset dm_tycho = stars_dmag  if (stars_mag < 12)
+  subset fl_tycho = stars_flags if (stars_mag < 12)
+
+  write $ROOT.tycho.dat -f "%10.6f %10.6f %6.3f %6.3f 0x%x" r_tycho d_tycho m_tycho dm_tycho fl_tycho
+
+  tapEXEC addstar -ref $ROOT.tycho.dat -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass -photcode TYCHO.g
+
+  $KAPA = kapa -noX
+  dev -n tmpdev
+  section default -imtool none
+  resize 50 50
+  region $RA_CENTER $DEC_CENTER 0.2
+  create tmp 0 1
+  cdensify tmpmap tmp tmp
+  set map = zero(tmpmap) + 0.2
+  close
+  $KAPA = kapa
+
+  local filter  
+  foreach filter g r i z y
+    if ("$filter" == "g")
+      wd map $ROOT.synthzpts.fits -bitpix -32 -bzero 0.0 -bscale 1.0
+    end
+    keyword map EXTNAME -w "map_$filter"
+    wd map $ROOT.synthzpts.fits -bitpix -32 -bzero 0.0 -bscale 1.0 -extend
+  end
+end
+
+macro cktycho
+  if ($0 != 2)
+    echo "USAGE: cktycho (stage)"
+    break
+  end
+
+  if ("$1" == "addstar") 
+    skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} 
+    mextract ra dec mag xccd yccd dbflags photflags where (photcode == 3001)
+    tapOK {ra[] == r_tycho[]} "read in as many tycho stars as we generated"
+    
+    match2d -closest ra dec r_tycho d_tycho 0.001 -index1 index1 -index2 index2
+    subset tmp = index1 if (index1 >= 0)
+    tapOK {tmp[] == ra[]} "tycho sources match those generated"
+  end
+
+  if ("$1" == "relphot")
+    subset rs = r_tycho if (m_tycho < 13.5)
+    subset ds = d_tycho if (m_tycho < 13.5)
+    avmatch rs ds 0.1 g:flags
+    subset tmp = rs if (g:flags & 0x04); # tycho used for photometry
+    tapOK {tmp[] == rs[]} "tycho objects used for photometry properly marked"
+  end
 end