Viewing contents of file '../idllib/contrib/buie/basphote.pro'
;+
; NAME:
; basphote
; PURPOSE: (one line)
; Circular aperture photometry extraction from images.
; DESCRIPTION:
;
; CATEGORY:
; CCD data processing
; CALLING SEQUENCE:
; Basphote,gain,image,exptime,xloc,yloc,radius,sky1,sky2,logfile[,objnum]
;
; INPUTS:
; gain : Gain of the CCD. Photons per count (DN).
; image : CCD image array.
; exptime : Exposure time for the image in seconds.
; xloc, yloc : Current location of the cursor in image coordinates.
; radius : Current aperture radius in pixels.
; sky1 : Inner radius of the sky annulus (in pixels).
; sky2 : Outer radius of the sky annulus (in pixels).
; logfile : Name of the photometry log file.
;
; OPTIONAL INPUT PARAMETERS:
; objnum : Starting serial number for object in frame, default=0
; If supplied as a scalar it will be incremented by
; the number of positions supplied in xloc/yloc.
; If given as a vector, will specify the exact number to
; and there will be no incrementing when done.
;
; KEYWORD INPUT PARAMETERS: (default on flags is false):
; AIRMASS - Optional airmass value.
; ALTLOG - Flag, if set, output is generated in an alternate format.
; BAD - Flag (or array) marking data as bad, default=good
; BOXMRAD - Size of the box to look for local max in. Default=radius.
; If boxmrad is negative, then the call to BOXM that finds the
; local max is suppressed. In effect, basphote assumes that
; the input location is already the maximum. You still need
; to provide a non-zero number for boxmrad so that a local
; area is defined for the fwhm calculation.
; DT - Delta-time, in seconds, between any two frames.
; EXACT - Flag, if true: take position as exact; otherwise find it.
; FNAME - File name of image.
; if ALTLOG set use the actual filename for the image on disk.
; if not set this should be an 8 character code that relates
; to the image (Default is blank).
; FILTER - Filter code, required if ALTLOG is set.
; JD - Julian date of observation (mid-time), required if ALTLOG set.
; NAME - Object name(s), required if ALTLOG set.
; NOLOG - If set, no output logfile information is generated.
; The default is to generate a log file.
; NOMEXT - Optional nominal extinction.
; PRINTALL- Print all objects even if off chip (if printing enabled)
; PSCALE - Plate scale in arc-sec per pixel, required if ALTLOG is set.
; SILENT - Flag, if true --> Do not generate any screen output.
; ZPOINT - Optional zero point.
; KEYWORD OUTPUT PARAMETERS:
; ERR - Optional return of the magnitude error.
; FLERR - Uncertainty (1 sigma) of object flux.
; FLUX - Object flux (photons per second)
; FWHM - FWHM of object image(s), in arcsec if PSCALE provided, otherwise
; returned in pixels
; ONCHIP - Byte array of flags that indicate if object was on-chip.
; OUTJD - Optional output of Julian dates (for cubes). This will be a
; vector of length rank(image). Contents are computed from the
; input Julian date (JD=) and the input delta-time (DT=).
; MAG - Optional return of the instrumental magnitude.
; MAX - Optional return of peak signal in object.
; SKYMEAN - Optional return of sky signal, counts/pixel.
; SKYERR - Optioanl return of sky signal uncertainty, counts/pixel.
; XCEN - Optional output of centroid x-position(s).
; YCEN - Optional output of centroid y-position(s).
;
; OUTPUTS:
; All output is sent to the screen and the logfile. Selected variables are
; returned via optional keywords.
;
; COMMON BLOCKS:
; None.
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
; Ported by Doug Loucks, Lowell Observatory, 1992 Oct from the
; C-language version written by Marc Buie.
;
; 1/6/93, Marc W. Buie, Lowell Observatory, added the alternate logfile
; format to the program. Streamlined the required input parameters.
; Calling sequence is now different.
;
; 1/25/93, MWB, changed default on boxmrad to be the LAST value in the
; radius vector.
;
; 2/3/93, DWL. Moved the file format declarations outside of the main loop.
; 2/4/93, DWL. Modified to accept all valid combinations of scalar and
; vector values for xloc, yloc, radius, sky1, and sky2 inputs.
; 2/4/93, DWL. Added checks for 'seeing' sub-array extraction.
; 2/4/93, DWL, Added keywords XCEN and YCEN to allow return of the refined
; center(s) to the calling program. Also, xcen, ycen, xmax, and ymax
; will be vectors if xloc and yloc are vectors.
; 2/10/93,DWL, Removed XMAX and YMAX keywords.
; 2/10/93,DWL, Added checks for off-chip and too-near-edge conditions.
; 4/21/93,DWL, Added keywords MAG and ERR. The instrumental magnitude and
; uncertainty are returned to the caller if these keywords are present.
; If the inputs xcen and ycen are vectors, these variables will also be
; returned as vectors.
; 5/21/93, MWB, Fixed variables that are undefined if off-chip condition
; is found. Dummy values are set to prevent printouts from dying.
; 5/26/93, MWB, Added FLUX, FLUXERR, and NOLOG keywords. Fixed bug on
; computation of sky background error, errant gain factor removed.
; 9/30/93, MWB, Modified usage of NAME to allow vector input.
; 2/2/1994, DWL, Mods to allow 3-D image cube photometry. Also added
; OUTJD keyword.
; 4/28/94, MWB, Fixed bug that recomputed a new centroid position if
; the EXACT flag was set. This had NO effect on photometry.
; 8/30/94, MWB, Added ONCHIP and PRINTALL keywords.
; 5/18/95, MWB, Added FWHM optional output keyword.
; 6/16/95, MWB, If PSCALE negative, fwhm not computed.
; 1/8/96, MWB, Added override on serial number, objnum can now be a vector.
; 7/13/96, MWB, changed FWHM calculation to remove function fitting.
; 10/31/96, MWB, added BAD= keyword for flag pass through
; 3/17/97, MWB, added MAX keyword
; 97/09/10, MWB, removed extraneous sub-expression evaluation (5x faster).
; 98/03/04, MWB, Added SKYMEAN, SKYERR return keywords.
; 98/09/21, MWB, added suppression of BOXM call through a negative BOXMRAD.
;-
PRO Basphote, gain, image, exptime, xloc, yloc, radius, sky1, sky2, $
logfile, in_objnum, $
AIRMASS=in_airmass, ALTLOG=altlog, BAD=bad, $
BOXMRAD=boxmrad, DT=in_dt, $
ERR=out_magerr, EXACT=exact, FILTER=filter, FLERR=out_fluxerr, $
FLUX=out_flux, FNAME=fname, FWHM=out_fwhm, $
JD=in_jd, OUTJD=out_jd, MAX=out_max, $
MAG=out_imag, NAME=name, NOLOG=nolog, $
NOMEXT=in_nomext, ONCHIP=out_onchip, SKYMEAN=out_skymean, $
SKYERR=out_skyerr, SILENT=silent, $
PSCALE=in_pscale, XCEN=out_xcen, YCEN=out_ycen, $
ZPOINT=in_zpoint ;, TIMEIT=timeit
;common com_basphote,info
;timeit = keyword_set(timeit)
;if timeit then cputime,time0
;IF timeit and n_elements(info) eq 0 THEN BEGIN
; info = { t: fltarr(20), ncalls: 0L }
; ; t[0] = startup overhead (to start of main FOR loop)
; ; t[1] = boxm
; ; t[2] = 1st center
; ; t[3] = sky 42%
; ; t[4] = final centrod call
; ; t[5] = FWHM
; ; t[6] = I/O
;ENDIF
;IF timeit THEN info.ncalls = info.ncalls+1
;Command line verification, 9 parameters are required, the tenth is optional.
IF N_PARAMS() LT 8 OR N_PARAMS() GT 10 THEN BEGIN
MESSAGE, /INFO, $
'basphote,gain,image,exptime,xloc,yloc,radius,sky1,sky2,logfile[,objnum]'
RETURN
ENDIF
IF badpar(gain,[2,3,4,5],0,CALLER='BASPHOTE (gain): ') THEN return
IF badpar(image,[1,2,3,4,5],[2,3],CALLER='BASPHOTE (image): ', $
RANK=im_rank, DIMEN=im_dimen) THEN return
IF badpar(exptime,[2,3,4,5],0,CALLER='BASPHOTE (exptime): ') THEN return
IF badpar(xloc,[2,3,4,5],[0,1],CALLER='BASPHOTE (xloc): ') THEN return
IF badpar(yloc,[2,3,4,5],[0,1],CALLER='BASPHOTE (yloc): ') THEN return
IF badpar(radius,[2,3,4,5],[0,1],CALLER='BASPHOTE (radius): ') THEN return
IF badpar(sky1,[2,3,4,5],[0,1],CALLER='BASPHOTE (sky1): ') THEN return
IF badpar(sky2,[2,3,4,5],[0,1],CALLER='BASPHOTE (sky2): ') THEN return
IF NOT KEYWORD_SET(nolog) THEN $
IF badpar(logfile,7,0,CALLER='BASPHOTE (logfile): ') THEN return
IF badpar(in_objnum,[0,2,3],[0,1],CALLER='BASPHOTE (objnum): ',DEF=0, $
RANK=num_rank, NPTS=num_npts) THEN return
IF badpar(boxmrad,[0,2,3,4,5],0,CALLER='BASPHOTE (boxmrad): ', $
DEF=radius[n_elements(radius)-1] ) THEN return
IF badpar(exact,[0,1,2,3],0,CALLER='BASPHOTE (exact): ',DEF=0) THEN return
IF badpar(printall,[0,1,2,3],0,CALLER='BASPHOTE (printall): ',DEF=0) THEN return
IF KEYWORD_SET(in_pscale) THEN pscale=in_pscale ELSE pscale=1.0
IF badpar(bad,[0,2,3],[0,1],CALLER='BASPHOTE (BAD): ',default=0) THEN error=1
IF KEYWORD_SET(altlog) THEN BEGIN
error=0
IF badpar(fname, 7, 0,CALLER='BASPHOTE (fname): ') THEN error=1
IF badpar(filter, 7, 0,CALLER='BASPHOTE (filter): ') THEN error=1
IF badpar(in_jd, 5, 0,CALLER='BASPHOTE (jd): ') THEN error=1
IF badpar(name, 7,[0,1],CALLER='BASPHOTE (name): ') THEN error=1
IF badpar(in_pscale,[4,5],0,CALLER='BASPHOTE (in_pscale): ') THEN error=1
IF error THEN RETURN
ENDIF
IF badpar( in_dt, [0,2,3,4,5], 0, CALLER='BASPHOTE (dt): ', DEFAULT=0.0D0) $
THEN RETURN
IF NOT KEYWORD_SET(fname) THEN fname=''
IF NOT KEYWORD_SET(name) THEN name=''
check=[N_ELEMENTS(xloc),N_ELEMENTS(radius), $
N_ELEMENTS(sky1),N_ELEMENTS(sky2),N_ELEMENTS(name),N_ELEMENTS(bad)]
z=WHERE(check ne 1, count)
IF N_ELEMENTS(xloc) NE N_ELEMENTS(yloc) THEN BEGIN
MESSAGE,'xloc and yloc must be the same length.', /INFO
RETURN
ENDIF
IF num_rank EQ 1 AND num_npts NE N_ELEMENTS(xloc) THEN BEGIN
MESSAGE,'serial number vector must match length of xloc and yloc', /INFO
RETURN
ENDIF
IF COUNT NE 0 THEN BEGIN
IF MIN(check[z]) NE MAX(check[z]) THEN BEGIN
IF N_ELEMENTS(xloc) EQ 1 THEN BEGIN
MESSAGE, /INFO, $
'radius, sky1, sky2, and name must be the same length or scalar.'
ENDIF ELSE BEGIN
MESSAGE, /INFO, $
'radius, sky1, sky2, name, and bad must match the length of xloc or be scalar.'
ENDELSE
RETURN
ENDIF
ENDIF
; Check for a few more keywords.
IF KEYWORD_SET( in_airmass ) THEN airmass=in_airmass ELSE airmass=0.0
IF KEYWORD_SET( in_nomext ) THEN nomext=in_nomext ELSE nomext=0.0
IF KEYWORD_SET( in_zpoint ) THEN zpoint=in_zpoint ELSE zpoint=0.0
IF KEYWORD_SET( in_jd ) THEN jd=in_jd ELSE jd=0.0D0
; Input image array size needs to be known.
xsize = im_dimen[ 0 ]
ysize = im_dimen[ 1 ]
IF im_rank EQ 3 THEN zsize=im_dimen[2] ELSE zsize=1
; Set-up the local copy of the Julian date.
IF zsize GT 1 THEN BEGIN
l_jd = jd + DINDGEN( zsize ) * ( DOUBLE( in_dt ) / 86400.0D0 )
ENDIF ELSE BEGIN
l_jd = jd
ENDELSE
; Make local copies of the inputs xloc, yloc, radius, sky1, and sky2.
; If any are vectors, the others of type scalar will be expanded into
; vectors having the same length. At this point, length requirements have
; been verified.
l_xloc = xloc
l_yloc = yloc
l_radius = radius
l_sky1 = sky1
l_sky2 = sky2
l_name = name
l_bad = bad
; Information about these inputs is needed to determine which are scalars
; and which are vectors.
xloc_stat = SIZE( xloc )
radius_stat = SIZE( radius )
sky1_stat = SIZE( sky1 )
sky2_stat = SIZE( sky2 )
name_stat = SIZE( name )
bad_stat = SIZE( bad )
;
largest = MAX( check )
IF largest GT 1 THEN BEGIN
; At least one is a vector with two or more elements. Make those of type
; scalar into vectors of the same length.
IF xloc_stat[0] EQ 0 THEN BEGIN
l_xloc = REPLICATE( xloc, largest )
l_yloc = REPLICATE( yloc, largest )
ENDIF
IF radius_stat[0] EQ 0 THEN l_radius = REPLICATE( radius, largest )
IF sky1_stat[0] EQ 0 THEN l_sky1 = REPLICATE( sky1, largest )
IF sky2_stat[0] EQ 0 THEN l_sky2 = REPLICATE( sky2, largest )
IF name_stat[0] EQ 0 THEN l_name = REPLICATE( name, largest )
IF bad_stat[0] EQ 0 THEN l_bad = REPLICATE( bad, largest )
ENDIF
l_xcen = FLTARR( largest, zsize )
l_ycen = FLTARR( largest, zsize )
localmax = FLTARR( largest, zsize )
imag = FLTARR( largest, zsize )
magerr = FLTARR( largest, zsize )
flux = FLTARR( largest, zsize )
fluxerr = FLTARR( largest, zsize )
fwhm = FLTARR( largest, zsize )
skymean = FLTARR( largest, zsize )
skyerr = FLTARR( largest, zsize )
l_onchip = BYTARR( largest, zsize )
; Open the log file for appending the new results.
IF NOT KEYWORD_SET(nolog) THEN $
OPENW, logfile_lu, logfile, /APPEND, /GET_LUN
; Define the log file formats.
IF KEYWORD_SET( altlog ) THEN BEGIN
fmt1 = '(a,1x,"''",a,"''",1x,a,1x,f13.5,1x,f8.3,1x,f6.2,1x,f7.3,1x,f7.3,' + $
'1x,f7.3,1x,i4.4,1x,f8.3,1x,f8.3,1x,f5.2,1x,f7.1,1x,f8.4,1x,f7.4,1x,i1)'
ENDIF ELSE BEGIN
fmt1 = '(A8, 1X, I4.4, 1X, F8.3, 1X, F8.3, 1X, F8.3, 1X, F5.2, 1X, F6.3, ' + $
'1X, F4.1, 1X, F5.1, 1X, F8.4, 1X, F6.4)'
ENDELSE
fmt2 = '(14X, F8.3, 1X, F8.3, 1X, F8.2, 1X, F6.2, 1X, F11.2, 1X, F8.2, ' + $
'1X, F5.2)'
boxsize = FIX( abs(boxmrad) + 0.5 )
;IF timeit THEN BEGIN
; cputime,time1
; info.t[0]=info.t[0]+(time1-time0)
;ENDIF
; Do the requested photometry for all inputs.
FOR frame=0, zsize-1 DO BEGIN
; Cull out the working image plane for processing. This is needed to
; avoid expensive array copies later.
IF zsize eq 1 THEN BEGIN
wimage = temporary(image)
ENDIF ELSE BEGIN
wimage = image[*,*,frame]
ENDELSE
;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
FOR i=0, largest-1 DO BEGIN
; if timeit then begin
; cputime,time1
; time2=time1
; time3=time1
; time4=time1
; time5=time1
; time6=time1
; time7=time1
; time8=time1
; time9=time1
; time10=time1
; endif
IF num_rank EQ 0 THEN BEGIN
objnum = in_objnum + i
ENDIF ELSE BEGIN
objnum = in_objnum[i]
ENDELSE
area = !PI * l_radius[i]*l_radius[i]
xcen = l_xloc[i]
ycen = l_yloc[i]
xcen1 = l_xloc[i]
ycen1 = l_yloc[i]
autosky = l_sky2[i] GT 0
; Set some initial variables to values to be printed if off-chip
; conditions
; are encountered.
photons = 0.0
tmp_skymean = 0.0
tmp_skyerr = 0.0
skyphot = 0.0
counts = 0.0
photerr = 0.0
sknfrac = 0.0
imag[i, frame] = 99.9999
magerr[i, frame] = 0.0
flux[i, frame] = 0.0
fluxerr[i, frame] = 0.0
fwhm[i, frame] = -1.0
skymean[i, frame] = 0.0
skyerr[i,frame] = 0.0
localmax[i, frame] = 0.0
; First, find the maximum within a box centered on the last known
;object position.
IF frame EQ 0 THEN BEGIN
;Use the provided input position.
objxloc = FIX( l_xloc[i] + 0.5 )
objyloc = FIX( l_yloc[i] + 0.5 )
ENDIF ELSE BEGIN
;Use the refined center from the last frame.
objxloc = FIX( l_xcen[i, frame-1] + 0.5 )
objyloc = FIX( l_ycen[i, frame-1] + 0.5 )
ENDELSE
; if timeit then cputime,time2
if boxmrad gt 0.0 then begin
boxm, wimage, objxloc, objyloc, boxsize, boxsize, xmax, ymax
endif else begin
xmax = objxloc
ymax = objyloc
endelse
; if timeit then cputime,time3
offchip = MIN([xmax, ymax, xsize-1-xmax, ysize-1-ymax]) LT 0
IF offchip THEN BEGIN
onchip = 0B
GOTO, skipcalc
ENDIF
localmax[i,frame] = wimage[ xmax, ymax ]
IF exact THEN BEGIN
; Use the provided position as the location of the aperture.
xcen = l_xloc[i]
ycen = l_yloc[i]
ENDIF ELSE BEGIN
; Refine the center by computing a center of mass of object aperture,
; no sky correction. Use the location of the brightest pixel near the
; provided location for the center.
tmp_skymean = 0.0
too_near_edge = ( xmax + l_radius[i] GT xsize - 0.5 ) OR $
( ymax + l_radius[i] GT ysize - 0.5 ) OR $
( xmax - l_radius[i] LT -0.5 ) OR ( ymax - l_radius[i] LT -0.5 )
IF too_near_edge THEN BEGIN
xcen1 = xmax
ycen1 = ymax
onchip = 0B
GOTO, skipcalc
ENDIF
centrod, wimage, xmax, ymax, l_radius[i], 0.0, 0.0, $
tmp_skymean, xcen, ycen, counts
ENDELSE
; if timeit then cputime,time4
; Compute the sky signal.
IF autosky THEN BEGIN
;Sky is to be determined.
Getannul, wimage, xcen, ycen, l_sky1[i], l_sky2[i], skybuf
Robomean, skybuf, 3.0, 0.5, tmp_skymean, avgdev, tmp_skyerr, $
var, skew, kurt, nsky
IF nsky NE 0 THEN BEGIN
tmp_skyerr = tmp_skyerr / SQRT( FLOAT( nsky ) )
ENDIF ELSE BEGIN ; This should never happen, warn if it does.
tmp_skyerr = 0.0
MESSAGE,'WARNING: Sky annulus is empty or entirely off chip',/INFO
ENDELSE
ENDIF ELSE BEGIN
; Sky is being provided externally.
tmp_skymean = l_sky1[i]
tmp_skyerr = ABS( l_sky2[i] )
ENDELSE
; if timeit then cputime,time5
; Compute the final result for the object.
too_near_edge = ( xcen + l_radius[i] GT xsize - 0.5 ) OR $
( ycen + l_radius[i] GT ysize - 0.5 ) OR $
( xcen - l_radius[i] LT -0.5 ) OR ( ycen - l_radius[i] LT -0.5 )
IF too_near_edge THEN BEGIN
xcen1 = xcen
ycen1 = ycen
onchip = 0B
GOTO, skipcalc
ENDIF
; if timeit then cputime,time6
centrod, wimage, xcen, ycen, l_radius[i], 0.0, 0.0, tmp_skymean, $
xcen1, ycen1, counts
; if timeit then cputime,time7
IF exact THEN BEGIN
xcen1 = xcen
ycen1 = ycen
ENDIF
photons = counts * gain
skyphot = tmp_skymean * gain
sigphotsky = area * tmp_skyerr * gain
varphotobj = photons + sigphotsky * sigphotsky
photerr = SQRT( varphotobj + sigphotsky * sigphotsky )
flux[i, frame] = photons / exptime
fluxerr[i, frame] = photerr / exptime
skymean[i,frame] = tmp_skymean
skyerr[i,frame] = tmp_skyerr
onchip = 1B
IF varphotobj NE 0.0 THEN BEGIN
sknfrac = sigphotsky * sigphotsky / varphotobj
ENDIF ELSE BEGIN
sknfrac = 0.0
ENDELSE
IF photons GT 0.0 THEN BEGIN
imag[i, frame] = -2.5 * ALOG10( flux[i, frame] ) + 24.0
magerr[i, frame] = 1.085736205 * photerr / photons
ENDIF ELSE BEGIN
imag[i, frame] = 99.999
magerr[i, frame] = 0.0
ENDELSE
; if timeit then cputime,time8
; Find out what the seeing was on the object.
; First, extract a sub-box around the max.
lx = xmax - boxsize
ly = ymax - boxsize
hx = xmax + boxsize
hy = ymax + boxsize
IF lx LT 0 THEN lx = 0
IF ly LT 0 THEN ly = 0
IF hx GE xsize THEN hx = xsize - 1
IF hy GE ysize THEN hy = ysize - 1
sub = wimage[ lx : hx, ly : hy ]
IF pscale GT 0.0 THEN BEGIN
xsum=fltarr(hx-lx+1)
ysum=fltarr(hy-ly+1)
FOR ii=0,hy-ly DO xsum = xsum+sub[*,ii]-tmp_skymean
FOR ii=0,hx-lx DO ysum = ysum+sub[ii,*]-tmp_skymean
xwd = total(xsum)/max(xsum)
ywd = total(ysum)/max(ysum)
tmp_fwhm = mean([xwd,ywd])
; radp, sub, xcen1-lx, ycen1-ly, r, im, tmp_fwhm
ENDIF ELSE BEGIN
tmp_fwhm = 1.0
ENDELSE
fwhm[i,frame] = tmp_fwhm*abs(pscale)
if fwhm[i,frame] gt 99.994 or fwhm[i,frame] lt -9.98 then fwhm[i,frame] = -1.0
; if timeit then cputime,time9
skipcalc:
IF ( NOT KEYWORD_SET(nolog) ) AND (printall OR onchip) THEN BEGIN
IF KEYWORD_SET( altlog ) THEN BEGIN
; New log file format
PRINTF, logfile_lu, FORMAT=fmt1, $
fname, l_name[i], filter, l_jd[frame], exptime, gain, l_radius[i], $
l_sky1[i], l_sky2[i], objnum, xcen1, ycen1, fwhm[i,frame], $
localmax[i,frame], imag[i, frame], magerr[i, frame], l_bad[i]
ENDIF ELSE BEGIN
; Old log file format
PRINTF, logfile_lu, FORMAT=fmt1, $
fname, objnum, xcen1, ycen1, exptime, gain, l_radius[i], $
l_sky1[i], l_sky2[i], imag[i, frame], magerr[i, frame]
PRINTF, logfile_lu, FORMAT=fmt2, $
xcen, ycen, skyphot, tmp_skyerr, photons, photerr, sknfrac
ENDELSE
ENDIF
; Send the results to the screen if wanted.
IF ( NOT KEYWORD_SET(silent) ) AND (printall OR onchip) THEN BEGIN
PRINT, ' '
PRINT, ' Image ', fname, ' Object ', l_name[i]
PRINT, 'Serial # ', objnum, $
', ', gain, ' e-/DN', $
', ', exptime, ' sec', $
', obj=', l_radius[i], $
FORMAT='(A,I4.4, A,G0.2,A, A,G0.1,A, A,G0.1,$)'
IF autosky THEN BEGIN
PRINT, ', sky (',l_sky1[i],',',l_sky2[i],')', $
FORMAT='(A,G0.1,A,G0.1,A)'
ENDIF ELSE BEGIN
PRINT, ', fixed sky.'
ENDELSE
PRINT, ' Position x y Peak counts= ', $
localmax[i,frame], ' (', localmax[i,frame]-tmp_skymean, $
' above sky)', FORMAT='(A,G0.2,A,G0.2,A)'
PRINT, ' Input', l_xloc[i], l_yloc[i], $
' Object ', counts, ' +/- ', photerr/gain, ' counts', $
FORMAT='(A,I10,I9,A,F9.2,A,F6.2,A)'
IF offchip THEN BEGIN
PRINT, ' \\\\\\\\ off-chip position ////////'
GOTO, continue
ENDIF
PRINT, ' Maximum', xmax, ymax, photons/exptime, ' +/- ', $
photerr/exptime, ' phot/sec', FORMAT='(A,I10,I9,12X,F9.2,A,F6.2,A)'
PRINT, ' Light #1', xcen, ycen, $
' Sky ', tmp_skymean, ' +/- ', tmp_skyerr, ' counts/pix', $
FORMAT='(A,F10.3,F9.3,A,F9.2,A,F6.2,A)'
PRINT, ' Light #2', xcen1, ycen1, $
tmp_skymean/exptime, ' +/- ', tmp_skyerr/exptime, $
' counts/pix/sec', FORMAT='(A,F10.3,F9.3,12X,F9.2,A,F6.2,A)'
PRINT, ' (sky noise ratio ', sknfrac, ')', $
tmp_skymean*gain/exptime*area, ' +/- ', tmp_skyerr*gain/exptime*area, $
' phot/sec', $
FORMAT='(A,F0.3,A,11X,F9.2,A,F6.2,A)'
IF KEYWORD_SET(in_pscale) THEN BEGIN
PRINT,' FWHM: ',fwhm[i,frame]/pscale,' pixels, ', $
fwhm[i,frame],' arc-sec', $
FORMAT='(a,1x,f5.2,1x,a,2x,f5.2,1x,a)'
ENDIF ELSE BEGIN
PRINT,' FWHM: ',fwhm[i,frame],' pixels', $
FORMAT='(a,1x,f6.2,1x,a)'
ENDELSE
PRINT, ' '
PRINT, '===> ', imag[i, frame], ' +/- ', magerr[i, frame], $
' <=== Instrumental magnitude', FORMAT='(A,F8.4,A,F8.4,A)'
IF (nomext NE 0 OR zpoint NE 0) AND airmass LT 9.0 THEN BEGIN
smag = imag[i, frame] - nomext * airmass + zpoint
PRINT, '===> ', smag, ' +/- ', magerr[i, frame], $
' <=== Provisional standard magnitude', FORMAT='(A,F8.4,A,F8.4,A)'
ENDIF
ENDIF
l_xcen[i, frame] = xcen1
l_ycen[i, frame] = ycen1
l_onchip[i, frame] = onchip
; if timeit then cputime,time10
continue:
; IF timeit THEN BEGIN
; IF time2-time1 GT 0. THEN info.t[0] = info.t[0] + (time2-time1)
; IF time3-time2 GT 0. THEN info.t[1] = info.t[1] + (time3-time2)
; IF time4-time3 GT 0. THEN info.t[2] = info.t[2] + (time4-time3)
; IF time5-time4 GT 0. THEN info.t[3] = info.t[3] + (time5-time4)
; IF time6-time5 GT 0. THEN info.t[0] = info.t[0] + (time6-time5)
; IF time7-time6 GT 0. THEN info.t[4] = info.t[4] + (time7-time6)
; IF time8-time7 GT 0. THEN info.t[0] = info.t[0] + (time8-time7)
; IF time9-time8 GT 0. THEN info.t[5] = info.t[5] + (time9-time8)
; IF time10-time9 GT 0. THEN info.t[6] = info.t[6] + (time10-time9)
; ENDIF
ENDFOR
;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
; For simple 2-d image, resurrect original copy.
IF zsize eq 1 THEN BEGIN
image = temporary(wimage)
ENDIF
ENDFOR
;if timeit then cputime,time1
IF NOT KEYWORD_SET(nolog) THEN FREE_LUN, logfile_lu
; Return some computations if the keywords are defined.
IF MAX( [ largest, zsize ] ) EQ 1 THEN BEGIN
out_xcen = l_xcen[0,0]
out_ycen = l_ycen[0,0]
out_imag = imag[0,0]
out_magerr = magerr[0,0]
out_flux = flux[0,0]
out_fluxerr = fluxerr[0,0]
out_fwhm = fwhm[0,0]
out_skymean = skymean[0,0]
out_skyerr = skyerr[0,0]
out_max = localmax[0,0]
out_onchip = l_onchip[0,0]
out_sky = skymean[0,0]
ENDIF ELSE BEGIN
IF (largest GT 1) AND (zsize GT 1) THEN BEGIN
out_xcen = l_xcen
out_ycen = l_ycen
out_imag = imag
out_magerr = magerr
out_flux = flux
out_fluxerr = fluxerr
out_fwhm = fwhm
out_skymean = skymean
out_skyerr = skyerr
out_max = localmax
out_onchip = l_onchip
ENDIF ELSE BEGIN
out_xcen = l_xcen[*]
out_ycen = l_ycen[*]
out_imag = imag[*]
out_magerr = magerr[*]
out_flux = flux[*]
out_fluxerr = fluxerr[*]
out_fwhm = fwhm[*]
out_skymean = skymean[*]
out_skyerr = skyerr[*]
out_max = localmax[*]
out_onchip = l_onchip[*]
ENDELSE
ENDELSE
out_jd = l_jd
IF num_rank EQ 0 THEN BEGIN
; Update serial number to point at the next number.
in_objnum = in_objnum + largest
ENDIF
;if timeit then begin
; cputime,time2
; info.t[0] = info.t[0] + (time2-time1)
; IF info.ncalls eq 100 THEN BEGIN
; alltime = total(info.t)
; print,'BASPHOTE: ',alltime,info.t[0:6]/alltime*100, $
; format='(a,f6.1,7(1x,f4.1))'
; info.ncalls = 0
; info.t[*] = 0.0
; ENDIF
;endif
END