## ----------------------------------------------------- ##
## Reducing the (scanty) 2002 April/May Lick/Gemini data ##
## ----------------------------------------------------- ##
                        steve, 30 April 2002
                        only worrying about K' data for now


########################################
## First swap FTS extension for .fits ##
########################################

  ls > filenames
  awk '{ printf("mv %s %s.fits\n",$1,$1)}' filenames > runme
  :.,$s/\.FTS\.fits/\.fits/g
  source runme


####################
## Make Sky Flats ##
####################

Skip darks and twi flats.  Make sky flat from science frames.
Handle standards and science frames separately.  First do 
science frames.

1.  Make list of science frames: list.science.
2.  Check numbers: iterstat @list.science > list.mean
3.  Plot up means in IDL  -->  all look useable.
4.  Combine:

imcombine ("@list.science",
"skyflat_k_10x12", rejmask="", plfile="", sigma="", logfile="STDOUT",
combine="average", reject="minmax", project=no, outtype="real",
offsets="none", masktype="none", maskvalue=0., blank=0., scale="none",
zero="median", weight="none", statsec=" ", expname="", lthreshold=INDEF,
hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes, lsigma=3., hsigma=3.,
rdnoise="6", gain="3.8", snoise="0.", sigscale=0.1, pclip=-0.5, grow=0.)

5.  Check it:

                display skyflat_k_10x12 1
                display 02MYB165 2
                display 02MYB050 3
                beep

6.  Normalize.
7.  Flatten science frames:  imarith @list.science / skyflat_k_10x12_nrm @list.science//ff


##############
## Make BPM ##
##############

imcopy skyflat_k_10x12_nrm bpm_K
imreplace bpm_K 0 lower=INDEF upper=0.5
imreplace bpm_K 1 lower=0.5 upper=1.5
imreplace bpm_K 0 lower=1.5 upper=INDEF
imcopy bpm_K.fits bpm_K.pl

** Handle standards identically but separately **

imhead @list.std
iterstat @list.std > list.stdmean

imcombine ("@list.std",
"skyflat_k_10x1", rejmask="", plfile="", sigma="", logfile="STDOUT",
combine="average", reject="minmax", project=no, outtype="real",
offsets="none", masktype="none", maskvalue=0., blank=0., scale="none",
zero="median", weight="none", statsec=" ", expname="", lthreshold=INDEF,
hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes, lsigma=3., hsigma=3.,
rdnoise="6", gain="3.8", snoise="0.", sigscale=0.1, pclip=-0.5, grow=0.)

display skyflat_k_10x1 1
display 02MYB042 2
display 02MYB196 3
beep

iterstat skyflat_k_10x1
imarith skyflat_k_10x1 / 2856.697 skyflat_k_10x1_nrm
display skyflat_k_10x1_nrm 4

imcopy skyflat_k_10x1_nrm bpm_Kstd
imreplace bpm_Kstd 0 lower=INDEF upper=0.5
imreplace bpm_Kstd 1 lower=0.5 upper=1.5
imreplace bpm_Kstd 0 lower=1.5 upper=INDEF
imcopy bpm_Kstd.fits bpm_Kstd.pl

imarith @list.std / skyflat_k_10x12_nrm @list.std//ff


############
## Dimsum ##
############

1.  Shoot crap off into subdirectories.  Since we targeted multiple radio
    sources in single pointings, the data will be organized by pointing.
    (I'm pretty sure the pointings are unique, so no spectral index in the
    directory name.)

source list.mkdir
ls *ff.fits > list.reduce
vi list.reduce
source list.reduce
:.,$s/02MY/mv 02MY/g

2.  Run reduce in each.

cd J1434+3427/
ls *fits > list.objects
vi list.objects
:.,$s/.fits//g

reduce ("@list.objects",
"J1434+3427", no, ref_img="02MYB121ff", prefix="ss_", chk=yes, fp_xslm=yes,
fp_fixpix=yes, fp_xzap=yes, fp_badpixupd=yes, mk_shifts=yes, fp_registar=yes,
masking=yes, maskdereg=yes, mp_xslm=yes, mp_fixpix=yes, mp_xzap=yes,
mp_badpixupd=yes, mp_registar=yes, nmed=6, bpimage="bpm_K.pl", mfac=2,
secpexp=1., ncoavg=1, rjct=1, nskymin=3, forcescale=yes, premask=yes,
pixin=yes, scale_x=1., scale_y=1., ang_a2x=0., senseab=yes, ang_y2N=0.,
rotation=yes, nsig_obj_msk=1.1, nsig_cr_msk=1.5, ilist="")


################
## Astrometry ##
################


cd /sleater/sdawson/lick/gemini/02apr/n4_reduced/J1435+3459

display J1435+3459_mp 1
display dss_plate_finder 2

imexam > ref.xy
vi ref.xy

tvmark 2 ref.xy mark=circle radii=5 color=206 label=no number=yes

imexam > image.xy
vi image.xy

xyeq yes dss_plate_finder "" ref.xy iraf 1 2 0 > ref.radec

cd,'/sleater/sdawson/lick/gemini/02apr/n4_reduced/J1435+3459'
matchpos

ccmap ("image.xyradec",
"image.db", solutions="", images="J1435+3459_mp", results="image.summary",
xcolumn=1, ycolumn=2, lngcolumn=3, latcolumn=4, xmin=INDEF, xmax=INDEF,
ymin=INDEF, ymax=INDEF, lngunits="hours", latunits="degrees",
insystem="j2000", refpoint="coords", lngref="INDEF", latref="INDEF",
refsystem="INDEF", lngrefunits="", latrefunits="", projection="tan",
fitgeometry="general", function="polynomial", xxorder=2, xyorder=2,
xxterms="half", yxorder=2, yyorder=2, yxterms="half", reject=INDEF,
update=yes, pixsystem="logical", verbose=yes, interactive=yes,
graphics="stdgraph", cursor="")

imexam


################
## Photometry ##
################

1.  First figure out what DIMSUM does by way of combining images.  Find some
    counts in sky-subtracted single frame, then compare to the combined frame.
    We're going to use tasks in noao.digiphot.daophot to do this, so you're
    going to need to set up photometry parameters as follows.

    a)	Measure your FWHM for a couple of stars (could use daoedit "a", cut and paste
	into a file, then "tstat FILE c5").  Measure sky and RMS in image (could use
	iterstat).  I find:

		FWHM = 4.0  (from stars in a couple of images)
   		SKY  = 0.08 (iterstat mean on std image)
		RMS  = 1.40 (iterstat RMS on std image)
	
	Things like readnoise and gain were taken from the Gemini manual.

    b)	Now set each of the following with ":e" in the noao.digiphot.daophot.phot:

	DATAPARS:

	datamin = SKY - 7*RMS = -9.72
	datamax = saturation 150000 electrons / 24 electrons/DN = 6250 DN (for Gemini LW chip)
      	readnoise = 40 electrons / sqrt(2) (based on sqrt(2) improvement over listed 8-read LW value)
      	gain = 24 electrons/DN (for Gemini LW chip)
		
       (scale = 1.)             Image scale in units per pixel
     (fwhmpsf = 4.)             FWHM of the PSF in scale units
    (emission = yes)            Features are positive ?
       (sigma = 1.4)            Standard deviation of background in counts
     (datamin = -9.72)          Minimum good data value
     (datamax = 6250.)          Maximum good data value
       (noise = "poisson")      Noise model
     (ccdread = "")             CCD readout noise image header keyword
        (gain = " ")            CCD gain image header keyword
   (readnoise = 28.)            CCD readout noise in electrons
       (epadu = 24.)            Gain in electrons per count
    (exposure = "ITIME")        Exposure time image header keyword
     (airmass = "AIRMASS")      Airmass image header keyword
      (filter = "FILTER")       Filter image header keyword
     (obstime = "TIME_OBS")     Time of observation image header keyword
       (itime = INDEF)          Exposure time
    (xairmass = INDEF)          Airmass
     (ifilter = " ")            Filter
       (otime = " ")            Time of observation
        (mode = "ql")           

	CENTERPARS:

	cbox = 2*FWHM = 8.0

  (calgorithm = "centroid")     Centering algorithm
        (cbox = 8.)             Centering box width in scale units
  (cthreshold = 0.)             Centering threshold in sigma above background
  (minsnratio = 1.)             Minimum signal-to-noise ratio for centering alg
    (cmaxiter = 10)             Maximum iterations for centering algorithm
    (maxshift = 1.)             Maximum center shift in scale units
       (clean = no)             Symmetry clean before centering
      (rclean = 1.)             Cleaning radius in scale units
       (rclip = 2.)             Clipping radius in scale units
      (kclean = 3.)             K-sigma rejection criterion in skysigma
    (mkcenter = no)             Mark the computed center
        (mode = "ql")           

	FITSKYPARS:

	annulus = 4*FWHM = 16.0
	dannulu = 3*FWHM = 12.0

  (salgorithm = "median")       Sky fitting algorithm
     (annulus = 16.)            Inner radius of sky annulus in scale units
    (dannulus = 12.)            Width of sky annulus in scale units
    (skyvalue = 0.)             User sky value
    (smaxiter = 10)             Maximum number of sky fitting iterations
     (sloclip = 0.)             Lower clipping factor in percent
     (shiclip = 0.)             Upper clipping factor in percent
    (snreject = 50)             Maximum number of sky fitting rejection iterati
   (sloreject = 3.)             Lower K-sigma rejection limit in sky sigma
   (shireject = 3.)             Upper K-sigma rejection limit in sky sigma
       (khist = 3.)             Half width of histogram in sky sigma
     (binsize = 0.10000000149012) Binsize of histogram in sky sigma
      (smooth = no)             Boxcar smooth the histogram
       (rgrow = 0.)             Region growing radius in scale units
       (mksky = no)             Mark sky annuli on the display
        (mode = "ql")           

	PHOTPARS:

	apertures = the apertures you want

   (weighting = "constant")     Photometric weighting scheme
   (apertures = "4.,6.,8.,10.,12.") List of aperture radii in scale units
        (zmag = 25.)            Zero point of magnitude scale
     (mkapert = yes)            Draw apertures on the display
        (mode = "ql")           

2.  We're going to look at what's going on in the standard image.  Use daofind
    to pick off some objects:

	image   =              FS27_mp  Input image(s)
	output  =             FS27.coo  Output coordinate file(s) (default: image.coo.?)
	(starmap=                     ) Output density enhancement image(s)
	(skymap =                     ) Output sky image(s)
	(datapar=                     ) Data dependent parameters
	(findpar=                     ) Object detection parameters
	(boundar=              nearest) Boundary extension (constant|nearest|reflect|wrap)
	(constan=                   0.) Constant for boundary extension
	(interac=                   no) Interactive mode ?
	(verify =            )_.verify) Verify critical daofind parameters ?
	(update =            )_.update) Update critical daofind parameters ?
	(verbose=           )_.verbose) Print daofind messages ?
	(graphic=          )_.graphics) Graphics device
	(display=           )_.display) Display device
	(icomman=                     ) Image cursor: [x y wcs] key [cmd]
	(gcomman=                     ) Graphics cursor: [x y wcs] key [cmd]
	(mode   =                   ql)

    ...then tvmark them: tvmark 1 FS27.coo.  Daofind kind of sucks, so you'll want
    to do some editing of the coordinate file.

3.  Do the same in one of the sky--subtracted images.  Fuck around until you've got
    coordinate lists with the same objects in the same order.  (Don't forget to 
    change the SKY and FWHM, etc.) Check it, for one sky-subtracted frame:

       (scale = 1.)             Image scale in units per pixel
     (fwhmpsf = 4.)             FWHM of the PSF in scale units
    (emission = yes)            Features are positive ?
       (sigma = 12.)            Standard deviation of background in counts
     (datamin = -83.)           Minimum good data value
     (datamax = 6250.)          Maximum good data value
       (noise = "poisson")      Noise model
     (ccdread = "")             CCD readout noise image header keyword
        (gain = " ")            CCD gain image header keyword
   (readnoise = 28.)            CCD readout noise in electrons
       (epadu = 24.)            Gain in electrons per count
    (exposure = "ITIME")        Exposure time image header keyword
     (airmass = "AIRMASS")      Airmass image header keyword
      (filter = "FILTER")       Filter image header keyword
     (obstime = "TIME_OBS")     Time of observation image header keyword
       (itime = INDEF)          Exposure time
    (xairmass = INDEF)          Airmass
     (ifilter = " ")            Filter
       (otime = " ")            Time of observation
        (mode = "ql")           

	display FS27_mp 1
	tvmar 1 FS27_mp.coo 
	display ss_02MYB192ff.fits 2
	tvmar 2 ss_02MYB192ff.coo

4.  Now run phot with those sources in the combine frame and the individual frame
    and compare fluxes.  Don't forget to change parameters.  What I get for one
    of the ten objects in each case:

FS27_mp                557.331   400.361   10   FS27.coo               10      \
   557.199    400.205    -0.132  -0.156  0.020   0.021          0    NoError   \
   0.5367987      1.618507       0.8384053      1650   8        0    NoError   \
   10.            INDEF          Kp                     05:34:12               \
   4.00     617.8139      50.57147   590.6672      20.572 0.023 0    NoError  *\
   6.00     771.1645      113.3344   710.3267      20.371 0.028 0    NoError  *\
   8.00     858.4654      201.4037   750.3521      20.312 0.036 0    NoError  *\
   10.00    939.3286      314.4582   770.5278      20.283 0.045 0    NoError  *\
   12.00    1076.496      452.6089   833.5359      20.198 0.051 0    NoError  *

ss_02MYB192ff          238.134   200.315   10   ss_02MYB192ff.coo      10      \
   238.255    200.157    0.121   -0.158  0.009   0.007          0    NoError   \
   1.140791       12.1612        6.777579       411    3        0    NoError   \
   10.            INDEF          Kp                     05:36:27               \
   2.00     559.1352      12.80232   544.5305      20.660 0.089 0    NoError  *\
   3.00     721.0471      28.60178   688.4184      20.405 0.106 0    NoError  *\
   4.00     805.3165      50.55906   747.6392      20.316 0.133 0    NoError  *\
   5.00     866.1939      78.83659   776.2578      20.275 0.165 0    NoError  *\
   6.00     1006.775      113.3066   877.5158      20.142 0.181 0    NoError  *

    The fluxes are the same!  That means that DIMSUM is averaging.  Also, I independently
    grabbed some fluxes with "p" in ATV on the combined image and got the same result.
    That means that the flux values are straight counts in the image; phot isn't doing
    some screwy thing with the exposure time.  Sweet!

5.  Now to concentrate on just the standard star.  Center on it, and grab the flux
    for a ton of aperatures:

FS27_mp                456.387   315.910   1    FS27.coo               1       \
   456.344    315.897    -0.043  -0.013  0.006   0.006          0    NoError   \
   0.1790568      1.04927        0.4157563      1654   5        0    NoError   \
   10.            INDEF          Kp                     05:34:12               \
   4.00     4940.655      50.53754   4931.605      18.268 0.004 0    NoError  *\
   6.00     5643.696      113.4907   5623.374      18.125 0.004 0    NoError  *\
   8.00     5930.317      201.57     5894.224      18.074 0.004 0    NoError  *\
   10.00    6099.444      314.4081   6043.147      18.047 0.005 0    NoError  *\
   12.00    6229.538      452.5243   6148.51       18.028 0.005 0    NoError  *\
   14.00    6301.636      615.6567   6191.398      18.021 0.006 0    NoError  *\
   16.00    6347.73       804.581    6203.664      18.018 0.007 0    NoError  *\
   18.00    6421.383      1018.452   6239.022      18.012 0.008 0    NoError  *\
   20.00    6516.447      1257.208   6291.335      18.003 0.009 0    NoError  *\
   22.00    6570.181      1520.629   6297.902      18.002 0.010 0    NoError  *\
   24.00    6611.678      1809.585   6287.659      18.004 0.011 0    NoError  *\
   26.00    6623.013      2123.667   6242.756      18.012 0.013 0    NoError  *\
   28.00    6679.888      2463.368   6238.805      18.012 0.015 0    NoError  *

   The apertures fluxes flatten out at about 6100 (where the magnitude errors become
   comparable to the differences in successive magnitudes).  From there on out, the
   flux is 6240 +/- 50 (the error here is just the RMS).

6.  Calculate zero point.  We know the following for the standard:

	K' = 13.102 mag
	flux = 6240 counts
	time = 10 seconds

    The K' magnitude is the interpolated value from http://www.ast.cam.ac.uk/~mhoenig/datadisk/faintstd_kprime.html.
    The flux was calculated above.  The exposure time comes from the fact that each individual 
    pointing was 10 seconds, and we learned in steps 1 - 4 that DIMSUM does an average, so
    that the combined image should be thought of as having the same exposure time as the
    individual images.  The zero point obeys:

	mag = -2.5 log (flux / time) + zpt 	--->	zpt = mag + 2.5 log (flux / time)

    Therefore, the K' zero point is 20.09.  (Should do error propagation...)

7.  Measure some objects.

    Now you can set it up so that phot returns the right magnitudes!  With the exposure time
    set to 10 seconds (under datapars) and the zero point set to 20.09 (under photpars), phot
    now returns:

FS27_mp                456.387   315.910   1    FS27.coo               1       \
   456.344    315.897    -0.043  -0.013  0.006   0.006          0    NoError   \
   0.1790568      1.04927        0.4157563      1654   5        0    NoError   \
   10.            INDEF          Kp                     05:34:12               \
   4.00     4940.655      50.53754   4931.605      13.358 0.004 0    NoError  *\
   6.00     5643.696      113.4907   5623.374      13.215 0.004 0    NoError  *\
   8.00     5930.317      201.57     5894.224      13.164 0.004 0    NoError  *\
   10.00    6099.444      314.4081   6043.147      13.137 0.005 0    NoError  *\
   12.00    6229.538      452.5243   6148.51       13.118 0.005 0    NoError  *\
   14.00    6301.636      615.6567   6191.398      13.111 0.006 0    NoError  *\
   16.00    6347.73       804.581    6203.664      13.108 0.007 0    NoError  *\
   18.00    6421.383      1018.452   6239.022      13.102 0.008 0    NoError  *\
   20.00    6516.447      1257.208   6291.335      13.093 0.009 0    NoError  *\
   22.00    6570.181      1520.629   6297.902      13.092 0.010 0    NoError  *\
   24.00    6611.678      1809.585   6287.659      13.094 0.011 0    NoError  *\
   26.00    6623.013      2123.667   6242.756      13.102 0.013 0    NoError  *\
   28.00    6679.888      2463.368   6238.805      13.102 0.015 0    NoError  * 

    So when you go to play with some objects, just change the exposure time...