IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 12, 2015, 6:18:23 PM (11 years ago)
Author:
eugene
Message:

merge changes from EAM dev branch ohana.20150429

Location:
trunk/Ohana
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/addstar/test

  • trunk/Ohana/src/addstar/test/relphot.flatcorr.dvo

    r33653 r38441  
    1414# * no other mechanism to tell DVO about a flat-field correction
    1515
     16$TAP_BREAK = 1
     17$TAP_VERBOSE = 1
    1618input tap.dvo
    1719
    1820# set globals
    19 if (not($?PLOT)) set PLOT = 0
     21if (not($?PLOT)) set PLOT = 0
     22
     23$USE_CELL_OFFSETS = 1; # turn on non-zero cell offsets
     24$SET_UC_RESID = 1; # write a non-zero resid value to the ubercal table (should not have an impact)
     25$DO_2MASS = 1
     26$DO_SYNTH = 1
     27$DO_TYCHO = 1
     28$NC_CLOUDS = 0.08
     29$UC_CLOUDS = 0.02
    2030
    2131# set various global variables
    2232macro init
     33  # generate a fake ptolemy.rc (to make this a contained test)
     34  exec rm -f ptolemy.rc
     35  output ptolemy.rc
     36    echo "\#"
     37    echo PHOTCODE_FILE dvo.photcodes
     38    echo GSCFILE GSCregions.tbl
     39    echo ZERO_PT 25.0
     40    echo ADDSTAR_RADIUS 1.0
     41    echo CATMODE                 SPLIT
     42    echo CATFORMAT               PS1_V5
     43    # echo CATFORMAT               PS1_V2 -- must be at least PS1_V3 for ubercaldist to work
     44    echo SKY_DEPTH               4
     45    echo CAMERA                  gpc1
     46    echo CAMERA_CONFIG           dvo.layout
     47    echo "input dvo.config"
     48  output stdout
     49
     50  # we need to define a photcode table consistent with the fake construction
     51  $zpt_nominal:g = 24.0
     52  $zpt_nominal:r = 24.5
     53
     54  # klam is negative, so klam*(airmass - 1) increases the zero point
     55  $klam_nominal:g = -0.15
     56  $klam_nominal:r = -0.10
     57
     58  # generate a fake dvo.photcode
     59  exec rm -f dvo.photcodes
     60  output dvo.photcodes
     61    echo "\# code  name                type    zero  slope offset c1    c2   slope   zero  equiv  sys scale   scale  sys     poor   bad     poor   bad"
     62    echo "1     g                    sec   0.000  0.000 0.000     1     3 0.0000     0  1051   0.000 0.000 0.000  0.000   0x0000 0x0000  0x0000 0x0000"
     63    echo "2     r                    sec   0.000  0.000 0.000     2     3 0.0000     0  1052   0.000 0.000 0.000  0.000   0x0000 0x0000  0x0000 0x0000"
     64    echo "10000 GPC1.g.XY00          dep  $zpt_nominal:g $klam_nominal:g 0.000     -     - 0.0000     0     1   0.015 0.000 1.000  0.000   0x0000 0x0000  0xe0440130 0x1003bc88"
     65    echo "10001 GPC1.g.XY01          dep  $zpt_nominal:g $klam_nominal:g 0.000     -     - 0.0000     0     1   0.015 0.000 1.000  0.000   0x0000 0x0000  0xe0440130 0x1003bc88"
     66    echo "10010 GPC1.g.XY10          dep  $zpt_nominal:g $klam_nominal:g 0.000     -     - 0.0000     0     1   0.015 0.000 1.000  0.000   0x0000 0x0000  0xe0440130 0x1003bc88"
     67    echo "10011 GPC1.g.XY11          dep  $zpt_nominal:g $klam_nominal:g 0.000     -     - 0.0000     0     1   0.015 0.000 1.000  0.000   0x0000 0x0000  0xe0440130 0x1003bc88"
     68    echo "10100 GPC1.r.XY00          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"
     69    echo "10101 GPC1.r.XY01          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"
     70    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"
     71    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"
     72    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"
     73    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"
     74    echo "2020  TYCHO_B              ref   0.000  0.000 0.000     -     - 0.0000     0     -   0.000 0.000 0.000  0.000   0x0000 0x0000  0x0000 0x0000"
     75  output stdout     
     76
    2377  $RA_CENTER = 10.0
    2478  $DEC_CENTER = 20.0
    … …  
    3286  $CHIP_DX = 1000
    3387  $CHIP_DY = 1000
    34 
    35   # images are loaded into dvo with GPC1 photcodes, so we need to get the nominal zps for those filters
    36   $zpt_nominal:g = 24.58
    37   $zpt_nominal:r = 24.80
    38   $zpt_nominal:i = 24.74
    39   $zpt_nominal:z = 24.26
    40   $zpt_nominal:y = 23.41
    41 
    42   # klam is negative, so klam*(airmass - 1) increase the zero point
    43   $klam_nominal:g = -0.15
    44   $klam_nominal:r = -0.10
    45   $klam_nominal:i = -0.04
    46   $klam_nominal:z = -0.03
    47   $klam_nominal:y = -0.03
     88  $DEFECT_FRAC = 0.05
     89  $OFFSET_FRAC_UC = 0.0
     90  $OFFSET_FRAC_NC = 0.5
    4891
    4992  # we have two sets of images: ubercaled and not-ubercaled
    5093
    51   # sequence to count images (only used in this function)
    52   create tmpseq 0 9
    53 
    54   $filt_uc:n = 9
    55   $filt_uc:0 = g
    56   $filt_uc:1 = g
    57   $filt_uc:2 = g
    58   $filt_uc:3 = r
    59   $filt_uc:4 = r
    60   $filt_uc:5 = r
    61   $filt_uc:6 = i
    62   $filt_uc:7 = i
    63   $filt_uc:8 = i
    64 
    65   # airmass slopes for these sequnece
     94  # the number of exposures is defined by filt_uc:n
     95  list filt_uc -split g g g g r r r r
     96  $filt_uc:n = 4
     97
     98  # tmpseq is used to generate vectors in this function
     99  create tmpseq 0 $filt_uc:n
     100
     101  # exptime and airmasses for uc images
     102  set exptime_uc = 10.0 + zero(tmpseq)
     103  set airmass_uc =  1.3 + zero(tmpseq)
    66104  set klam_uc = zero(tmpseq)
    67   for i 0 tmpseq[]
    68     klam_uc[$i] = $klam_nominal:$filt_uc:$i
    69   end   
    70 
    71   # mjd and zpt values for ubercal'ed imaged
    72   set exptime_uc = 10.0 + zero(tmpseq)
    73   set airmass_uc = 1.3 + zero(tmpseq)
     105  set zpt_uc = zero(tmpseq)
    74106
    75107  # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0)
    76108  # note that K is defined as a negative value (is this sensible?)
    77   set zpt_uc = 25.0 + 2.5*log(exptime_uc) + klam_uc*(airmass_uc - 1.0) + tmpseq*0.0050 - 0.0025
    78 
    79   set mjd_uc = zero(zpt_uc)
    80   mjd_uc[0] = 55000.01
    81   mjd_uc[1] = 55000.02
    82   mjd_uc[2] = 55000.03
    83   mjd_uc[3] = 55015.01
    84   mjd_uc[4] = 55015.02
    85   mjd_uc[5] = 55015.03
    86   mjd_uc[6] = 55025.01
    87   mjd_uc[7] = 55025.02
    88   mjd_uc[8] = 55025.03
     109
     110  # zero points and airmass slopes for these sequnece
     111  for i 0 $filt_uc:n
     112    klam_uc[$i] = $klam_nominal:$filt_uc:$i
     113    zpt_uc[$i]  = $zpt_nominal:$filt_uc:$i + 2.5*log(exptime_uc[$i]) + klam_uc[$i]*(airmass_uc[$i] - 1.0) + $UC_CLOUDS*(rnd(klam_uc[$i]) - 0.5)
     114  end   
     115  set mjd_uc = 55000.00 + 0.01*(tmpseq % 3) + 10.0*int(tmpseq / 3)
    89116
    90117  # mjd and zpt values for not-ubercal'ed images
    91   $filt_nc:n = 9
    92   $filt_nc:0 = g
    93   $filt_nc:1 = g
    94   $filt_nc:2 = g
    95   $filt_nc:3 = r
    96   $filt_nc:4 = r
    97   $filt_nc:5 = r
    98   $filt_nc:6 = i
    99   $filt_nc:7 = i
    100   $filt_nc:8 = i
     118  list filt_nc -split g g g g r r r r
     119  $filt_nc:n = 4
     120
     121  # tmpseq is used to generate vectors in this function
     122  create tmpseq 0 $filt_nc:n
     123
     124  # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0)
     125  set exptime_nc = 10.0 + zero(tmpseq)
     126  set airmass_nc =  1.6 + zero(tmpseq)
     127  set klam_nc = zero(tmpseq)
     128  set zpt_nc = zero(tmpseq)
    101129
    102130  # airmass slopes for these sequnece
    103   set klam_nc = zero(tmpseq)
    104   for i 0 tmpseq[]
    105     klam_nc[$i] = $klam_nominal:$filt_nc:$i
     131  for i 0 $filt_nc:n
     132    klam_nc[$i] = $klam_nominal:$filt_nc:$i
     133    zpt_nc[$i]  = $zpt_nominal:$filt_nc:$i + 2.5*log(exptime_nc[$i]) + klam_nc[$i]*(airmass_nc[$i] - 1.0) + $NC_CLOUDS*(rnd(klam_nc[$i]) - 0.5)
    106134  end   
    107 
    108   # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0)
    109   set exptime_nc = 15.0 + zero(tmpseq)
    110   set airmass_nc = 1.6 + zero(tmpseq)
    111 
    112   # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0)
    113   set zpt_nc = 25.0 + 2.5*log(exptime_nc) + klam_nc*(airmass_nc - 1.0) - tmpseq*0.0050 + 0.0025
    114 
    115   set mjd_nc = zero(zpt_nc)
    116   mjd_nc[0] = 55000.11
    117   mjd_nc[1] = 55000.12
    118   mjd_nc[2] = 55000.13
    119   mjd_nc[3] = 55015.11
    120   mjd_nc[4] = 55015.12
    121   mjd_nc[5] = 55015.13
    122   mjd_nc[6] = 55025.11
    123   mjd_nc[7] = 55025.12
    124   mjd_nc[8] = 55025.13
     135  set mjd_nc = 55000.10 + 0.01*(tmpseq % 3) + 10.0*int(tmpseq / 3)
     136
     137  delete tmpseq
    125138end
    126139
    … …  
    139152  mkdir $rootdir
    140153
     154  break -auto off
    141155  exec rm -rf $catdir
     156  exec rm -rf $fileroot*
     157  break -auto on
    142158
    143159  init
    144160  mkzptfile
    145   mkstars 500
     161  mkstars 1000
     162
     163  catdir $catdir
     164
     165  mk2mass $fileroot $catdir
     166  ck2mass addstar
     167 
     168  mksynth $fileroot $catdir
     169  cksynth addstar
     170
     171  # to add a valid tycho test, I need to generate a zpt table for tycho as well
     172  mktycho $fileroot $catdir
     173  cktycho addstar
    146174
    147175  tapPLAN {((mjd_uc[]*4 + mjd_nc[]*4) + (mjd_uc[] + mjd_nc[]) + (mjd_uc[] + mjd_nc[])) / 3}
    … …  
    149177  # generate the basic images and check they were correctly ingested by dvo
    150178  mksequence $fileroot $catdir
    151   for i 0 mjd_uc[] 3
    152     ckexposure catdir.test mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw
    153   end
    154   for i 0 mjd_nc[] 3
    155     ckexposure catdir.test mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw
     179  for i 0 mjd_uc[]
     180    ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw
     181  end
     182  for i 0 mjd_nc[]
     183    ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw
    156184  end
    157185
    158186  # run setphot on the db and check that the images now match the expected values
    159   exec setphot -update -ubercal testzpt.fits -D CATDIR $catdir
    160   for i 0 mjd_uc[] 3
    161     ckexposure catdir.test mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc
    162   end
    163   for i 0 mjd_nc[] 3
    164     ckexposure catdir.test mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc
     187  tapEXEC setphot -update -ubercal testzpt.fits -D CATDIR $catdir
     188  for i 0 mjd_uc[]
     189    ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc
     190  end
     191  for i 0 mjd_nc[]
     192    ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc
    165193  end
    166194
    167195  # run relphot on the db and check that the images now match the expected values
    168   exec relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR catdir.test -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& tmp.log
    169   for i 0 mjd_uc[] 3
    170     ckexposure catdir.test mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot
    171   end
    172   for i 0 mjd_nc[] 3
    173     ckexposure catdir.test mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot
    174   end
     196  tapEXEC relphot -images g,r -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -D IMAGE_OFFSET 0.5 -update -synthphot_zpts $fileroot.synthzpts.fits
     197  for i 0 mjd_uc[]
     198    ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot_uc
     199  end
     200  for i 0 mjd_nc[]
     201    ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc
     202  end
     203  ck2mass relphot
     204  cksynth relphot
     205  cktycho relphot
    175206
    176207  tapDONE
     208end
     209
     210macro try.relphot
     211  if ($0 != 3)
     212    echo "go (fileroot) (catdir)"
     213    break
     214  end
     215
     216  local i fileroot rootdir catdir
     217
     218  $fileroot = $1
     219  $catdir = $2
     220
     221  # run relphot on the db and check that the images now match the expected values
     222  tapEXEC relphot -images g,r -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -D IMAGE_OFFSET 0.5 -update -synthphot_zpts $fileroot.synthzpts.fits -nloop 24
     223
     224  for i 0 mjd_uc[]
     225    # ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot_uc
     226      ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot_uc
     227  end
     228  for i 0 mjd_nc[]
     229    # ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc
     230      ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc
     231  end
    177232end
    178233
    … …  
    215270  $AIRMASS     = $5
    216271  $FILTER      = $6
    217   $MODE        = $7
     272  $FUNC_MODE   = $7
    218273
    219274  # XXX need a function to extract the nominal zpt for a given filter / photcode from the db
    … …  
    235290
    236291  skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2}
    237   mextract ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001)
     292  mextract ra dec mag xccd yccd dbflags photflags where (abs(time - $MJD_IMAGE) < 0.0001)
    238293 
    239294  match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2
    … …  
    246301  set dm = mag - stars_mag_m
    247302
     303  set good_photom = (not(photflags & 0x02))
     304  foreach field dm xccd yccd
     305    subset $field\_good = $field if good_photom
     306  end
     307
    248308  if ($PLOT)
    249309    dev -n 1
    250     lim mag dm; clear; box; plot mag dm
    251   end
    252 
    253   set cell_xbin = int(xccd / ($CHIP_DX / $NCELL_X))
    254   set cell_ybin = int(yccd / ($CHIP_DY / $NCELL_Y))
    255   subset dm00 = dm if (cell_xbin == 0) && (cell_ybin == 0)
    256   subset dm01 = dm if (cell_xbin == 0) && (cell_ybin == 1)
    257   subset dm10 = dm if (cell_xbin == 1) && (cell_ybin == 0)
    258   subset dm11 = dm if (cell_xbin == 1) && (cell_ybin == 1)
     310    lim mag dm; clear; box; plot mag dm; plot -c red -pt 7 mag dm where (photflags & 0x00000002)
     311    dev -n 2
     312    lim xccd dm; clear; box; plot xccd dm; plot -c red -pt 7 xccd dm where (photflags & 0x00000002)
     313    dev -n 3
     314    lim ra dm; clear; box; plot ra dm; plot -c red -pt 7 ra dm where (photflags & 0x00000002)
     315    echo "red points are known bad photometry"
     316  end
     317
     318  # reproduce eddie's esoteric padding
     319  set cell_xbin = max(0 , min(1 , int((xccd_good +  8) / ($CHIP_DX / $NCELL_X))))
     320  set cell_ybin = max(0 , min(1 , int((yccd_good + 10) / ($CHIP_DY / $NCELL_Y))))
     321  subset dm00 = dm_good if (cell_xbin == 0) && (cell_ybin == 0)
     322  subset dm01 = dm_good if (cell_xbin == 0) && (cell_ybin == 1)
     323  subset dm10 = dm_good if (cell_xbin == 1) && (cell_ybin == 0)
     324  subset dm11 = dm_good if (cell_xbin == 1) && (cell_ybin == 1)
    259325
    260326  # uncorrected values behave like this:
    … …  
    273339  $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0)
    274340
    275   if ("$MODE" == "raw")
     341  if ("$FUNC_MODE" == "raw")
    276342    for ix 0 $NCELL_X
    277343      for iy 0 $NCELL_Y
    … …  
    280346 
    281347        set dm_adjust = dm$ix\$iy - $ZPT_NOMINAL + $ZPT_REAL_NORM + cell_off[$ix][$iy]
     348        # echo set dm_adjust = dm$ix\$iy - $ZPT_NOMINAL + $ZPT_REAL_NORM + cell_off[$ix][$iy]
     349
    282350        vstat -q dm_adjust
    283         tapOK {abs($MEAN) < 0.005} "addstar raw ZP $MJD_IMAGE $FILTER"
     351        sprintf line "addstar raw ZP %8.2f $FILTER %6.3f %6.3f %6.3f: %7.4f %7.4f : cell%d%d %6.3f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL $MEAN $SIGMA $ix $iy cell_off[$ix][$iy] dm_adjust[]
     352        tapOK {abs($MEAN) < 0.005} "$line"
    284353      end
    285354    end
    … …  
    287356  end
    288357
    289   if ("$MODE" == "setphot_uc")
    290     vstat -q dm
    291     tapOK {abs($MEAN) < 0.005} "setphot_uc  ZP $MJD_IMAGE $FILTER"
     358  if ("$FUNC_MODE" == "setphot_uc")
     359    vstat -q dm_good
     360    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[]
     361    tapOK {abs($MEAN) < 0.005} "$line"
     362
     363    subset tmp = dm if (dbflags & 0x00008000)
     364    tapOK {tmp[] == dm[]} "ubercal detections have ubercal bit set"
     365
     366    imextract MCAL FLAGS where (abs(time - $MJD_IMAGE) < 0.0001)
     367    subset tmp = FLAGS if (FLAGS & 0x0000200)
     368    tapOK {tmp[] == FLAGS[]} "ubercal images have ubercal bit set"
     369   
    292370    return
    293371  end
    294372
    295   if ("$MODE" == "setphot_nc")
    296     vstat -q dm
    297     tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "setphot_nc  ZP $MJD_IMAGE $FILTER"
     373  if ("$FUNC_MODE" == "setphot_nc")
     374    vstat -q dm_good
     375    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[]
     376    tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "$line"
     377
     378    subset tmp = dm if (not(dbflags & 0x00008000))
     379    tapOK {tmp[] == dm[]} "non-ubercal detections do NOT have ubercal bit set"
     380
     381    imextract FLAGS where (abs(time - $MJD_IMAGE) < 0.0001)
     382    subset tmp = FLAGS if (not(FLAGS & 0x0000200))
     383    tapOK {tmp[] == FLAGS[]} "non-ubercal images do NOT have ubercal bit set"
     384   
    298385    return
    299386  end
    300387
    301   if ("$MODE" == "relphot")
    302     vstat -q dm
    303     tapOK {abs($MEAN) < 0.005} "relphot    ZP $MJD_IMAGE $FILTER"
     388  if (("$FUNC_MODE" == "relphot_uc") || ("$FUNC_MODE" == "relphot_nc"))
     389    vstat -q dm_good
     390    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[]
     391
     392    if ("$FUNC_MODE" == "relphot_uc")
     393      tapOK {abs($MEAN) < 0.001} "$line"
     394      tapOK {abs($SIGMA) < 0.001} "$line"
     395    end
     396    if ("$FUNC_MODE" == "relphot_nc")
     397      # XXX NOTE relaxed condition : bad photometry points are partly clipped but still cause a bias
     398      tapOK {abs($MEAN) < 0.020} "$line"
     399      tapOK {abs($SIGMA) < 0.020} "$line"
     400    end
     401
     402    mextract ra dec mag dbflags
     403    subset rs = ra  if ((dbflags & 0x8000) && (mag > 13.5))
     404    subset ds = dec if ((dbflags & 0x8000) && (mag > 13.5))
     405    avmatch rs ds 0.1 $FILTER\:flags $FILTER\:uc_dist
     406
     407    subset tmp = rs if (($FILTER\:flags & 0x8) || ($FILTER\:flags & 0x4))
     408    tapOK {tmp[] == rs[]} "objects with ubercal detections have ubercal bit set or used SYNTH photometry (tmp[] vs rs[])"
     409 
     410    subset tmp = rs if (($FILTER\:uc_dist < 0.1) || ($FILTER\:flags & 0x4))
     411    tapOK {tmp[] == rs[]} "objects with ubercal detections have ubercal dist of 0.0 or used SYNTH photometry"
    304412    return
    305413  end
     414
     415  echo "ERROR: unknown mode for ckexposure : $FUNC_MODE"
    306416
    307417end
    … …  
    316426
    317427  # we have defined a set of mjd and zero point values for ubercal'ed images
    318   # for i 0 mjd_uc[]
    319   for i 0 mjd_uc[] 3
     428  for i 0 mjd_uc[]
    320429    mkexposure $1.uc.$i $RA_CENTER $DEC_CENTER zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] mjd_uc[$i] $filt_uc:$i $2
    321430  end
    322431
    323432  # we have a second set of mjd and zero point values for non-ubercal'ed images
    324   # for i 0 mjd_nc[]
    325   for i 0 mjd_nc[] 3
     433  for i 0 mjd_nc[]
    326434    mkexposure $1.nc.$i $RA_CENTER $DEC_CENTER zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2
    327435  end
    … …  
    334442 end
    335443
    336  local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR
     444 local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR cmd
    337445
    338446 $ROOT    = $1
    … …  
    366474    $options = -mjd $MJD
    367475    $options = $options -radec $ra $dec
    368     $options = $options -type PS1_V2
     476
     477    #$options = $options -type PS1_V2
     478    #$options = $options -type PS1_V3
     479    #$options = $options -type PS1_V4
     480    $options = $options -type PS1_V5
     481
    369482    $options = $options -coords
    370483    $options = $options -photcode GPC1.$FILTER.XY$ix\$iy
    371484    $options = $options -no-noise
     485    $options = $options -flags
     486    $options = $options -bad-psfqf-frac $DEFECT_FRAC
    372487    $options = $options -size $CHIP_DX $CHIP_DY
    373488    $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY}
    374489    $options = $options -airmass $AIRMASS
    375490    $options = $options -exptime $EXPTIME
    376     exec mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options
    377     echo mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options
    378    
     491
     492    tapEXEC mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options
     493
    379494    # the fake images have inconsistent ra,dec and airmass,sidtime values
    380     echo addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass
    381     exec addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass >& tmp.log
     495    tapEXEC addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass
    382496  end
    383497 end
    … …  
    405519  subset stars_mr = stars_mag if valid
    406520
    407   set cell_x = int(stars_x / ($CHIP_DX / $NCELL_X))
    408   set cell_y = int(stars_y / ($CHIP_DX / $NCELL_X))
     521  # note eddie's esoteric padding rule:
     522  set cell_x = max(0 , min(1 , int((stars_x +  8) / ($CHIP_DX / $NCELL_X))))
     523  set cell_y = max(0 , min(1 , int((stars_y + 10) / ($CHIP_DX / $NCELL_X))))
    409524
    410525  # eddie says: M_real = M_inst + zpt + ff_offset
    … …  
    418533  end
    419534
    420   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
     535  # create bad detections / detections to ignore or avoid, save to a file?
     536
     537  # create 3 types of things to ignore / avoid
     538  # * badPhot (deviate the photometry)
     539  # * badFlag (set a bad flag bit)
     540  # * poorFlag (set a bad flag bit)
     541  set badPhot  = rnd(stars_r) <  $DEFECT_FRAC
     542  set badFlag  = rnd(stars_r) <  $DEFECT_FRAC
     543  set poorFlag = rnd(stars_r) <  $DEFECT_FRAC
     544
     545  # create deviant photometry (1%)
     546  # I'm setting the flag value 0x02 (EXTMODEL) if the object is deviated.  This is a simple way to pass the info
     547  # note that outlier rejection only occurs for Nmeas > 5
     548  set stars_mi = stars_mi + 0.2*badPhot
     549  set stars_fl = zero(stars_mi) + 0x00000080*badFlag + 0x00000100*poorFlag + 0x00000002*badPhot
     550
     551  # bad flags that we sould respect:
     552  # PM_SOURCE_MODE_FAIL             = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero)
     553  # PM_SOURCE_MODE_SATSTAR          = 0x00000080, ///< Source model peak is above saturation
     554  # PM_SOURCE_MODE_BADPSF           = 0x00000400, ///< Failed to get good estimate of object's PSF
     555  # PM_SOURCE_MODE_DEFECT           = 0x00000800, ///< Source is thought to be a defect
     556  # PM_SOURCE_MODE_SATURATED        = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail)
     557  # PM_SOURCE_MODE_CR_LIMIT         = 0x00002000, ///< Source has crNsigma above limit
     558  # PM_SOURCE_MODE_MOMENTS_FAILURE  = 0x00008000, ///< could not measure the moments
     559  # PM_SOURCE_MODE_SKY_FAILURE      = 0x00010000, ///< could not measure the local sky
     560  # PM_SOURCE_MODE_SKYVAR_FAILURE   = 0x00020000, ///< could not measure the local sky variance
     561  # PM_SOURCE_MODE_SIZE_SKIPPED     = 0x10000000, ///< size could not be determined
     562
     563  # poor flags that we should respect:
     564  # PM_SOURCE_MODE_POOR             = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?)
     565  # PM_SOURCE_MODE_PAIR             = 0x00000020, ///< Source fitted with a double psf
     566  # PM_SOURCE_MODE_BLEND            = 0x00000100, ///< Source is a blend with other sources
     567  # PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N
     568  # PM_SOURCE_MODE_BLEND_FIT        = 0x00400000, ///< source was fitted as a blend
     569  # PM_SOURCE_MODE_ON_SPIKE         = 0x20000000, ///< peak lands on diffraction spike
     570  # PM_SOURCE_MODE_ON_GHOST         = 0x40000000, ///< peak lands on ghost or glint
     571  # PM_SOURCE_MODE_OFF_CHIP         = 0x80000000, ///< peak lands off edge of chip
     572
     573  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
    421574end
    422575
    … …  
    439592  set stars_ra  = $RA_CENTER  + $RA_RANGE  * (rnd(tmp) - 0.5) / 3600 / dcos ($DEC_CENTER)
    440593  set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) / 3600
    441   set stars_mag = 15.0 + 4.0 * rnd(tmp)
     594  set stars_mag = 10.0 + 10.0 * rnd(tmp)
     595  set stars_dmag = 0.02 + 0.1*rnd(tmp)
     596end
     597
     598macro mk2mass
     599  if ($0 != 3)
     600    echo "USAGE: mk2mass (fileroot) (catdir)"
     601    break
     602  end
     603
     604  if (not($DO_2MASS)) return
     605
     606  local ROOT CATDIR
     607
     608  $ROOT    = $1
     609  $CATDIR  = $2
     610
     611  set stars_flags = 0x00c00000*(rnd(stars_ra) < 0.1) + 0x7*(rnd(stars_ra) < 0.8)
     612  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
     613
     614  tapEXEC addstar -ref $ROOT.2mass.dat -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass -photcode 2MASS_J
     615
     616  set r_2mass = stars_ra
     617  set d_2mass = stars_dec
     618  set m_2mass = stars_mag
     619  set dm_2mass = stars_dmag
     620  set fl_2mass = stars_flags
     621end
     622
     623macro ck2mass
     624  if ($0 != 2)
     625    echo "USAGE: ck2mass (stage)"
     626    break
     627  end
     628
     629  if (not($DO_2MASS)) return
     630
     631  if ("$1" == "addstar")
     632    skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2}
     633    mextract ra dec mag xccd yccd dbflags photflags where (photcode == 2011)
     634    tapOK {ra[] == r_2mass[]} "read in as many 2mass stars as we generated"
     635   
     636    match2d -closest ra dec r_2mass d_2mass 0.001 -index1 index1 -index2 index2
     637    subset tmp = index1 if (index1 >= 0)
     638    tapOK {tmp[] == ra[]} "2mass sources match those generated"
     639  end
     640
     641  if ("$1" == "relphot")
     642    subset rs = r_2mass if (fl_2mass & 0x00c00000)
     643    subset ds = d_2mass if (fl_2mass & 0x00c00000)
     644    avmatch rs ds 0.1 flags
     645    subset tmp = rs if (flags & 0x02000000); # extended in 2mass
     646    tapOK {tmp[] == rs[]} "2mass extended objects properly marked"
     647   
     648    subset rs = r_2mass if (fl_2mass & 0x07)
     649    subset ds = d_2mass if (fl_2mass & 0x07)
     650    avmatch rs ds 0.1 flags
     651    subset tmp = rs if (flags & 0x08000000); # extended in 2mass
     652    tapOK {tmp[] == rs[]} "2mass good objects properly marked"
     653  end
     654end
     655
     656macro mksynth
     657  if ($0 != 3)
     658    echo "USAGE: mksynth (fileroot) (catdir)"
     659    break
     660  end
     661
     662  if (not($DO_SYNTH)) return
     663
     664  local ROOT CATDIR
     665
     666  $ROOT    = $1
     667  $CATDIR  = $2
     668 
     669  subset r_synth  = stars_ra    if (stars_mag < 17)
     670  subset d_synth  = stars_dec   if (stars_mag < 17)
     671  subset m_synth  = stars_mag   if (stars_mag < 17)
     672  subset dm_synth = stars_dmag  if (stars_mag < 17)
     673  subset fl_synth = stars_flags if (stars_mag < 17)
     674
     675  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
     676
     677  tapEXEC addstar -ref $ROOT.synth.dat -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass -photcode SYNTH.g
     678end
     679
     680macro cksynth
     681  if ($0 != 2)
     682    echo "USAGE: cksynth (stage)"
     683    break
     684  end
     685
     686  if (not($DO_SYNTH)) return
     687
     688  if ("$1" == "addstar")
     689    skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2}
     690    mextract ra dec mag xccd yccd dbflags photflags where (photcode == 3001)
     691    tapOK {ra[] == r_synth[]} "read in as many synth stars as we generated"
     692   
     693    match2d -closest ra dec r_synth d_synth 0.001 -index1 index1 -index2 index2
     694    subset tmp = index1 if (index1 >= 0)
     695    tapOK {tmp[] == ra[]} "synth sources match those generated"
     696  end
     697
     698  if ("$1" == "relphot")
     699    subset rs = r_synth if (m_synth < 13.5)
     700    subset ds = d_synth if (m_synth < 13.5)
     701    avmatch rs ds 0.1 g:flags
     702    subset tmp = rs if (g:flags & 0x04); # synth used for photometry
     703    tapOK {tmp[] == rs[]} "synth objects used for photometry properly marked"
     704  end
     705end
     706
     707macro mktycho
     708  if ($0 != 3)
     709    echo "USAGE: mktycho (fileroot) (catdir)"
     710    break
     711  end
     712
     713  if (not($DO_TYCHO)) return
     714
     715  local ROOT CATDIR
     716
     717  $ROOT    = $1
     718  $CATDIR  = $2
     719
     720  subset r_tycho  = stars_ra    if (stars_mag < 12)
     721  subset d_tycho  = stars_dec   if (stars_mag < 12)
     722  subset m_tycho  = stars_mag   if (stars_mag < 12)
     723  subset dm_tycho = stars_dmag  if (stars_mag < 12)
     724  subset fl_tycho = stars_flags if (stars_mag < 12)
     725
     726  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
     727
     728  tapEXEC addstar -ref $ROOT.tycho.dat -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass -photcode TYCHO_B
     729
     730  $KAPA = kapa -noX
     731  dev -n tmpdev
     732  # resize 100 100
     733  resize 50 50
     734  region $RA_CENTER $DEC_CENTER 0.2
     735  section default -imtool none
     736  # causes trouble for this to come before region...
     737  create tmp 0 1
     738  cdensify tmpmap tmp tmp
     739  set map = zero(tmpmap) + 0.2
     740  close
     741  $KAPA = kapa
     742
     743  local filter 
     744  foreach filter g r i z y
     745    if ("$filter" == "g")
     746      wd map $ROOT.synthzpts.fits -bitpix -32 -bzero 0.0 -bscale 1.0
     747    end
     748    keyword map EXTNAME -w "map_$filter"
     749    wd map $ROOT.synthzpts.fits -bitpix -32 -bzero 0.0 -bscale 1.0 -extend
     750  end
     751end
     752
     753macro cktycho
     754  if ($0 != 2)
     755    echo "USAGE: cktycho (stage)"
     756    break
     757  end
     758
     759  if (not($DO_TYCHO)) return
     760
     761  if ("$1" == "addstar")
     762    skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2}
     763    mextract ra dec mag xccd yccd dbflags photflags where (photcode == 2020)
     764    tapOK {ra[] == r_tycho[]} "read in as many tycho stars as we generated"
     765   
     766    match2d -closest ra dec r_tycho d_tycho 0.001 -index1 index1 -index2 index2
     767    subset tmp = index1 if (index1 >= 0)
     768    tapOK {tmp[] == ra[]} "tycho sources match those generated"
     769  end
     770
     771  if ("$1" == "relphot")
     772    subset rst = r_tycho if (m_tycho < 12)
     773    subset dst = d_tycho if (m_tycho < 12)
     774    # XXX : should not restrict to 'g'
     775    avmatch rst dst 0.1 g:flags
     776    subset tmp = rst if (g:flags & 0x40); # tycho used for photometry
     777    tapOK {tmp[] == rst[]} "tycho objects used for photometry properly marked"
     778
     779    # non-tycho stars with synth phot should be repaired
     780    subset rsn = r_synth if (m_synth > 12) && (m_synth < 13.5)
     781    subset dsn = d_synth if (m_synth > 12) && (m_synth < 13.5)
     782    # XXX : should not restrict to 'g'
     783    avmatch rsn dsn 0.1 g:flags
     784    subset tmp = rsn if (g:flags & 0x80); # synth photometry repaired
     785    tapOK {tmp[] == rsn[]} "tycho objects used for photometry properly marked"
     786  end
    442787end
    443788
    … …  
    456801
    457802  # define the season boundaries
    458   keyword phu TS0_0000 -wf 55000.0
    459   keyword phu TS1_0000 -wf 55010.0
    460   keyword phu TS0_0001 -wf 55010.0
    461   keyword phu TS1_0001 -wf 55020.0
    462   keyword phu TS0_0002 -wf 55020.0
    463   keyword phu TS1_0002 -wf 55030.0
     803  keyword phu S0_MJD0 -wf 55000.0
     804  keyword phu S1_MJD0 -wf 55010.0
     805  keyword phu S0_MJD1 -wf 55010.0
     806  keyword phu S1_MJD1 -wf 55020.0
     807  keyword phu S0_MJD2 -wf 55020.0
     808  keyword phu S1_MJD2 -wf 55030.0
    464809
    465810  # create a single flat-field correction map for a single chip
    466811  mcreate cell_off 2 2
    467   cell_off[0][0] =  0.01
    468   cell_off[1][0] = -0.01
    469   cell_off[0][1] =  0.02
    470   cell_off[1][1] = -0.02
     812  if ($USE_CELL_OFFSETS)
     813    cell_off[0][0] =  0.05
     814    cell_off[1][0] = -0.03
     815    cell_off[0][1] =  0.08
     816    cell_off[1][1] = -0.10
     817  end
    471818
    472819  # in the correction table, we have an image of the full exposure
    … …  
    481828    for iyc 0 $NCELL_Y
    482829      for ix 0 $NCHIP_X
     830        set tcell_off = cell_off
     831        if ($ix < 4)
     832          rotate tcell_off 180
     833        end
    483834        for ixc 0 $NCELL_X
    484           concat cell_off[$ixc][$iyc] image_offset
     835          concat tcell_off[$ixc][$iyc] image_offset
    485836        end
    486837      end
    … …  
    498849  set mjd_obs = mjd_uc
    499850  set zp = zpt_uc
    500   set resid = 0.02*rnd(zpt_uc) - 0.01
     851
     852  if ($SET_UC_RESID)
     853    set resid = 0.02*rnd(zpt_uc) - 0.01
     854  else
     855    set resid = zero(zpt_uc)
     856  end
    501857
    502858  delete -q flatcorr
Note: See TracChangeset for help on using the changeset viewer.