From 11f3517438a371e3a1b38b43694251d9e3407a8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E6=9D=8E=E7=90=BC?= Date: Thu, 2 Nov 2017 11:35:13 -0700 Subject: [PATCH 01/11] idl stage4 geom --- pro/kcwi_alt_cals.pro | 360 ++++++++++++++ pro/kcwi_apply_dgeom.pro | 226 +++++++++ pro/kcwi_apply_geom.pro | 295 ++++++++++++ pro/kcwi_central_wave.pro | 107 +++++ pro/kcwi_cfg__define.pro | 194 ++++++++ pro/kcwi_cfg_string.pro | 39 ++ pro/kcwi_dgeom__define.pro | 130 +++++ pro/kcwi_extract_arcs.pro | 217 +++++++++ pro/kcwi_fit_center.pro | 504 ++++++++++++++++++++ pro/kcwi_flatten_cube.pro | 134 ++++++ pro/kcwi_geom__define.pro | 198 ++++++++ pro/kcwi_geom_resid.pro | 74 +++ pro/kcwi_get_atlas.pro | 36 ++ pro/kcwi_get_imname.pro | 75 +++ pro/kcwi_group_dgeom.pro | 123 +++++ pro/kcwi_group_geom.pro | 143 ++++++ pro/kcwi_plot_arcfits.pro | 356 ++++++++++++++ pro/kcwi_print_cfgs.pro | 125 +++++ pro/kcwi_print_info.pro | 106 +++++ pro/kcwi_read_atlas.pro | 79 ++++ pro/kcwi_read_cfg.pro | 193 ++++++++ pro/kcwi_read_cfgs.pro | 112 +++++ pro/kcwi_read_ppar.pro | 178 +++++++ pro/kcwi_read_proc.pro | 187 ++++++++ pro/kcwi_reverse_geom.pro | 146 ++++++ pro/kcwi_set_dgeom.pro | 154 ++++++ pro/kcwi_set_geom.pro | 302 ++++++++++++ pro/kcwi_solve_arcs.pro | 692 +++++++++++++++++++++++++++ pro/kcwi_solve_arcs_iter.pro | 887 +++++++++++++++++++++++++++++++++++ pro/kcwi_solve_dgeom.pro | 285 +++++++++++ pro/kcwi_solve_geom.pro | 97 ++++ pro/kcwi_stage4geom.pro | 442 +++++++++++++++++ pro/kcwi_trace_cbars.pro | 567 ++++++++++++++++++++++ pro/kcwi_verify_cfg.pro | 57 +++ pro/kcwi_verify_dirs.pro | 120 +++++ pro/kcwi_verify_geom.pro | 56 +++ pro/kcwi_verify_ppar.pro | 56 +++ pro/kcwi_write_dgeom.pro | 82 ++++ pro/kcwi_write_geom.pro | 100 ++++ 39 files changed, 8234 insertions(+) create mode 100644 pro/kcwi_alt_cals.pro create mode 100644 pro/kcwi_apply_dgeom.pro create mode 100644 pro/kcwi_apply_geom.pro create mode 100644 pro/kcwi_central_wave.pro create mode 100644 pro/kcwi_cfg__define.pro create mode 100644 pro/kcwi_cfg_string.pro create mode 100644 pro/kcwi_dgeom__define.pro create mode 100644 pro/kcwi_extract_arcs.pro create mode 100644 pro/kcwi_fit_center.pro create mode 100644 pro/kcwi_flatten_cube.pro create mode 100644 pro/kcwi_geom__define.pro create mode 100644 pro/kcwi_geom_resid.pro create mode 100644 pro/kcwi_get_atlas.pro create mode 100644 pro/kcwi_get_imname.pro create mode 100644 pro/kcwi_group_dgeom.pro create mode 100644 pro/kcwi_group_geom.pro create mode 100644 pro/kcwi_plot_arcfits.pro create mode 100644 pro/kcwi_print_cfgs.pro create mode 100644 pro/kcwi_print_info.pro create mode 100644 pro/kcwi_read_atlas.pro create mode 100644 pro/kcwi_read_cfg.pro create mode 100644 pro/kcwi_read_cfgs.pro create mode 100644 pro/kcwi_read_ppar.pro create mode 100644 pro/kcwi_read_proc.pro create mode 100644 pro/kcwi_reverse_geom.pro create mode 100644 pro/kcwi_set_dgeom.pro create mode 100644 pro/kcwi_set_geom.pro create mode 100644 pro/kcwi_solve_arcs.pro create mode 100644 pro/kcwi_solve_arcs_iter.pro create mode 100644 pro/kcwi_solve_dgeom.pro create mode 100644 pro/kcwi_solve_geom.pro create mode 100644 pro/kcwi_stage4geom.pro create mode 100644 pro/kcwi_trace_cbars.pro create mode 100644 pro/kcwi_verify_cfg.pro create mode 100644 pro/kcwi_verify_dirs.pro create mode 100644 pro/kcwi_verify_geom.pro create mode 100644 pro/kcwi_verify_ppar.pro create mode 100644 pro/kcwi_write_dgeom.pro create mode 100644 pro/kcwi_write_geom.pro diff --git a/pro/kcwi_alt_cals.pro b/pro/kcwi_alt_cals.pro new file mode 100644 index 0000000..0d6e8e8 --- /dev/null +++ b/pro/kcwi_alt_cals.pro @@ -0,0 +1,360 @@ +; +; Copyright (c) 2017, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_ALT_CALS +; +; PURPOSE: +; Find calibration files in an alternate directory. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; Result = KCWI_ALT_CALS(Kcfg, Adir, Ppar) +; +; INPUTS: +; Kcfg - KCWI_CFG struct giving configuration of target +; Adir - Alternate calibration source directory +; Ppar - KCWI_PPAR struct for pipeline params +; +; KEYWORDS: +; +; OUTPUTS: +; Full file specification of alternate calibration file. +; +; SIDE EFFECTS: +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2018-JUN-27 Initial version +;- +function kcwi_alt_cals, kcfg, adir, ppar, $ + bias=bias, dark=dark, flat=flat, geom=geom, dgeom=dgeom, $ + prof=prof, rr=rr, drr=drr, std=std + ; + ; setup + pre = 'KCWI_ALT_CALS' + cfile = '' + ; + ; verify inputs + if kcwi_verify_cfg(kcfg,/init) ne 0 then return,cfile + if kcwi_verify_ppar(ppar,/init) ne 0 then return,cfile + ; + ; test adir + if not file_test(adir,/directory,/read) then begin + kcwi_print_info,ppar,pre,'Alt cal dir not accessable: '+adir, $ + /error + return,cfile + endif + ; + ; log + kcwi_print_info,ppar,pre,systime(0) + ; + ; bias + if keyword_set(bias) then begin + kcwi_print_info,ppar,pre,'Searching for mbiases in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_mbias.fit*', count=nf) + ; + ; found some biases + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['CCDMODE','AMPMODE','XBINSIZE','YBINSIZE', $ + 'GAINMUL','NAXIS1','NAXIS2'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist,count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest ccd temperature to match + endif else begin + tdel = abs(mcfg.tmpa1 - kcfg.tmpa1) + t = (where(tdel eq min(tdel)))[0] + cfile = mcfg[t].obsdir + mcfg[t].obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching biases found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate biases found' + ; + ; darks + endif else if keyword_set(dark) then begin + kcwi_print_info,ppar,pre,'Searching for mdarks in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_mdark.fit*', count=nf) + ; + ; found some darks + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['CCDMODE','AMPMODE','XBINSIZE','YBINSIZE','GAINMUL'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist,count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest exposure time to match + endif else begin + tdel = abs(mcfg.exptime - kcfg.exptime) + tind = where(tdel eq min(tdel),ntind) + ; + ; same exposure time, use ccd temperature to match + if ntind gt 1 then begin + zcfg = mcfg[tind] + zdel = abs(zcfg.tmpa1 - kcfg.tmpa1) + zind = (where(zdel eq min(zdel)))[0] + mcfg = zcfg[zind] + endif else $ + mcfg = mcfg[tind] + cfile = mcfg.obsdir + mcfg.obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching darks found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate darks found' + ; + ; flats + endif else if keyword_set(flat) then begin + kcwi_print_info,ppar,pre,'Searching for mflats in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_mflat.fit*', count=nf) + ; + ; found some flats + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['CCDMODE','AMPMODE','XBINSIZE','YBINSIZE', $ + 'GAINMUL','GRATID','GRANGLE', 'FILTNUM', $ + 'CAMANG','IFUNUM'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist, $ + imgtype='cflat',count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest bench temperature to match + endif else begin + tdel = abs(mcfg.tmpa8 - kcfg.tmpa8) + t = (where(tdel eq min(tdel)))[0] + cfile = mcfg[t].obsdir + mcfg[t].obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching flats found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate flats found' + ; + ; geom files + endif else if keyword_set(geom) then begin + kcwi_print_info,ppar,pre,'Searching for geom files in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_geom.fit*', count=nf) + ; + ; found some geom files + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['XBINSIZE','YBINSIZE','GRATID','GRANGLE', $ + 'FILTNUM','CAMANG','IFUNUM'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist,count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest bench temperature to match + endif else begin + tdel = abs(mcfg.tmpa8 - kcfg.tmpa8) + t = (where(tdel eq min(tdel)))[0] + cfile = mcfg[t].obsdir + mcfg[t].obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching geom files found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate geom files found' + ; + ; direct geom files + endif else if keyword_set(dgeom) then begin + kcwi_print_info,ppar,pre,'Searching for dgeom files in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_dgeom.fit*', count=nf) + ; + ; found some dgeom files + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['XBINSIZE','YBINSIZE','CAMANG','IFUNUM'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist,count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest bench temperature to match + endif else begin + tdel = abs(mcfg.tmpa8 - kcfg.tmpa8) + t = (where(tdel eq min(tdel)))[0] + cfile = mcfg[t].obsdir + mcfg[t].obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching dgeom files found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate dgeom files found' + ; + ; slice profiles + endif else if keyword_set(prof) then begin + kcwi_print_info,ppar,pre,'Searching for profs in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_prof.fit*', count=nf) + ; + ; found some profs + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['CCDMODE','AMPMODE','XBINSIZE','YBINSIZE', $ + 'GAINMUL','GRATID','GRANGLE', 'FILTNUM', $ + 'CAMANG','IFUNUM'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist,count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest bench temperature to match + endif else begin + tdel = abs(mcfg.tmpa8 - kcfg.tmpa8) + t = (where(tdel eq min(tdel)))[0] + cfile = mcfg[t].obsdir + mcfg[t].obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching profs found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate profs found' + ; + ; relative response files + endif else if keyword_set(rr) then begin + kcwi_print_info,ppar,pre,'Searching for rrs in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_rr.fit*', count=nf) + ; + ; found some relative response files + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['CCDMODE','AMPMODE','XBINSIZE','YBINSIZE', $ + 'GAINMUL','GRATID','GRANGLE', 'FILTNUM', $ + 'CAMANG','IFUNUM'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist,count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest bench temperature to match + endif else begin + tdel = abs(mcfg.tmpa8 - kcfg.tmpa8) + t = (where(tdel eq min(tdel)))[0] + cfile = mcfg[t].obsdir + mcfg[t].obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching rrs found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate rrs found' + ; + ; direct relative response files + endif else if keyword_set(drr) then begin + kcwi_print_info,ppar,pre,'Searching for drrs in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_drr.fit*', count=nf) + ; + ; found some direct relative response files + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['CCDMODE','AMPMODE','XBINSIZE','YBINSIZE', $ + 'GAINMUL','CAMANG','IFUNUM'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist,count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest bench temperature to match + endif else begin + tdel = abs(mcfg.tmpa8 - kcfg.tmpa8) + t = (where(tdel eq min(tdel)))[0] + cfile = mcfg[t].obsdir + mcfg[t].obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching drrs found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate drrs found' + ; + ; inverse sensitivity files + endif else if keyword_set(std) then begin + kcwi_print_info,ppar,pre,'Searching for invsens files in '+adir + cfgs = kcwi_read_cfgs(adir, filespec='*_invsens.fit*', count=nf) + ; + ; found some ivnsens files + if nf gt 0 then begin + ; + ; matching criteria + tlist = ['CCDMODE','AMPMODE','XBINSIZE','YBINSIZE', $ + 'GAINMUL','GRATID','GRANGLE', 'FILTNUM', $ + 'CAMANG','IFUNUM'] + mcfg = kcwi_match_cfg(cfgs,kcfg,ppar,tlist,count=nm) + ; + ; found some matches + if nm gt 0 then begin + ; + ; only one? + if nm eq 1 then begin + cfile = mcfg.obsdir + mcfg.obsfname + ; + ; more than one, use closest bench temperature to match + endif else begin + tdel = abs(mcfg.tmpa8 - kcfg.tmpa8) + t = (where(tdel eq min(tdel)))[0] + cfile = mcfg[t].obsdir + mcfg[t].obsfname + endelse + endif else $ + kcwi_print_info,ppar,pre, $ + 'No matching invsens files found' + endif else $ + kcwi_print_info,ppar,pre,'No alternate invsens files found' + endif + return,cfile +end ; kcwi_alt_cals diff --git a/pro/kcwi_apply_dgeom.pro b/pro/kcwi_apply_dgeom.pro new file mode 100644 index 0000000..efd8d67 --- /dev/null +++ b/pro/kcwi_apply_dgeom.pro @@ -0,0 +1,226 @@ +; +; Copyright (c) 2016, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_APPLY_DGEOM +; +; PURPOSE: +; Apply the final geometric transformation to an input image. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_APPLY_DGEOM, Img, Hdr, Kdgeom, Ppar, Dimg, Dhdr +; +; INPUTS: +; Img - Stage1 processed object image to apply geometry to +; Hdr - the corresponding header of the input image +; Kdgeom - KCWI_DGEOM struct after KCWI_SOLVE_DGEOM has been run +; +; INPUT KEYWORDS: +; VERBOSE - extra output +; +; OUTPUTS: +; Dimg - 2-D direct mode image +; Dhdr - A fits header with WCS info from Kdgeom +; +; SIDE EFFECTS: +; None. +; +; PROCEDURE: +; Apply the geomtric transformation in Kdgeom to Img and copy Hdr to +; Dhdr and update WCS keywords. +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2016-NOV-09 Initial Revision +;- +pro kcwi_apply_dgeom,img,hdr,kdgeom,ppar,dimg,dhdr +; +; startup +pre = 'KCWI_APPLY_DGEOM' +q = '' +; +; Check structs +if kcwi_verify_geom(kdgeom,/init) ne 0 then return +if kcwi_verify_ppar(ppar) ne 0 then begin + ppar = {kcwi_ppar} +endif +; +; get image original size +sz = size(img,/dim) +; +; output image +onx = 180 / kdgeom.xbinsize +ony = 300 / kdgeom.ybinsize +dimg = fltarr(onx,ony) +npix = fltarr(onx,ony) +maximx = 0 +; +; image number +imgnum = sxpar(hdr,'FRAMENO') +object = sxpar(hdr,'TARGNAME',count=ntarg) +if ntarg le 0 then $ + object = sxpar(hdr,'OBJECT') +imgtyp = sxpar(hdr,'CALTYPE') +; +; log +kcwi_print_info,ppar,pre,'Slicing and dicing image '+strn(imgnum)+': '+object+'...' +; +; loop over slices +for i=0,23 do begin + x0 = kdgeom.x0[i] + x1 = kdgeom.x1[i] + y0 = kdgeom.y0[i] + y1 = kdgeom.y1[i] + w = kdgeom.wid[i] + nxx = (x1-x0) + 1 + if nxx gt maximx then maximx = nxx + sub = img[x0:x1,y0:y1] + subr = rot(sub,kdgeom.angles[i],cubic=-0.5) + si = size(subr) + x = lindgen(si[1],si[2]) + y = x/si[1] + x = x MOD si[1] + res = interpolate(subr,x-kdgeom.xoff[i],y) + yc = si[2]/2 + if i gt 0 then begin + yp0 = yp1 + yp1 = yp0 + w*2 + endif else begin + yp0 = 0 + yp1 = w*2 + endelse + dimg[0:nxx-1,yp0:yp1] += res[0:nxx-1,yc-w:yc+w] + npix[0:nxx-1,yp0:yp1] += 1. + if ppar.verbose eq 1 then $ + print,strn(i)+' ',format='($,a)' +endfor +; +; normalize +t = where(npix eq 0., nt) +if nt gt 0 then $ + npix[t] = 1. +dimg = dimg/npix +maximy = yp1 +; +dimg = transpose(dimg[0:maximx, 0:maximy]) +if ppar.verbose eq 1 then begin + print,'Done.',format='($,a)' + print,'' +endif +; +; update header +dhdr = hdr +; +; image dimensions +sxaddpar,dhdr,'NAXIS',2 +sxaddpar,dhdr,'NAXIS1',maximy +sxaddpar,dhdr,'NAXIS2',maximx +; +; pixel scales +sxaddpar,dhdr,'PXSCL', kdgeom.pxscl*kdgeom.xbinsize,' Pixel scale along slice' +sxaddpar,dhdr,'SLSCL', kdgeom.slscl,' Pixel scale purpendicular to slices' +; +; geometry origins +sxaddpar,dhdr, 'CBARSFL', kdgeom.cbarsfname,' Continuum bars image' +sxaddpar,dhdr, 'ARCFL', kdgeom.arcfname, ' Arc image' +sxaddpar,dhdr, 'CBARSNO', kdgeom.cbarsimgnum,' Continuum bars image number' +sxaddpar,dhdr, 'ARCNO', kdgeom.arcimgnum, ' Arc image number' +sxaddpar,dhdr, 'GEOMFL', kdgeom.geomfile,' Geometry file' +; +; get sky coords +rastr = sxpar(hdr,'RA',count=nra) +decstr = sxpar(hdr,'DEC',count=ndec) +if nra ne 1 or ndec ne 1 then begin + rastr = sxpar(hdr,'TARGRA',count=nra) + decstr = sxpar(hdr,'TARGDEC',count=ndec) +endif +if nra eq 1 and ndec eq 1 then begin + radec_parse,rastr,decstr,':',ra,dec +endif else begin + ra = -99.d0 + dec = -99.d0 +endelse +; +; Position Angle ( = SKYPA) in degrees +; Plus an offset between rotator and IFU (may be zero) +skypa = sxpar(hdr,'ROTPOSN',count=npa) + sxpar(hdr,'ROTREFAN') +crota = -(skypa + kdgeom.rotoff) / !RADEG +sxaddpar,dhdr,'IFUPA',skypa,' IFU position angle (degrees)' +sxaddpar,dhdr,'IFUROFF',kdgeom.rotoff,' IFU-SKYPA offset (degrees)' +; +; pixel scales +cdelt1 = -kdgeom.pxscl ; RA degrees per px (column) +cdelt2 = kdgeom.pxscl*kdgeom.xbinsize ; Dec degrees per slice (row) +; +; did we get good coords? +if nra ne 1 or ndec ne 1 or npa ne 1 then begin + ; + ; no good coords + ; a warning for objects + if strcmp(imgtyp,'object') eq 1 then begin + kcwi_print_info,ppar,pre,'no coords for image',imgnum,imgtyp, $ + format='(a,2x,a,2x,a)',/warning + ; otherwise just info + endif else begin + kcwi_print_info,ppar,pre,'no coords for image',imgnum,imgtyp, $ + format='(a,2x,a,2x,a)' + endelse + ; + ; zero coords + ra = 0. + dec = 0. + ; + ; nominal CD matrix (no rotation) + CD11 = cdelt1*cos(0.) + CD12 = abs(cdelt2)*sign(cdelt1)*sin(0.) + CD21 = -abs(cdelt1)*sign(cdelt2)*sin(0.) + CD22 = cdelt2*cos(0.) +endif else begin + ; + ; calculate CD matrix + CD11 = cdelt1*cos(crota) ; RA degrees per column + CD12 = abs(cdelt2)*sign(cdelt1)*sin(crota) ; RA degress per row + CD21 = -abs(cdelt1)*sign(cdelt2)*sin(crota) ; DEC degress per column + CD22 = cdelt2*cos(crota) ; DEC degrees per row +endelse +; +; get reference pixels +if ppar.crpix1 le 0. then $ + crpix1 = maximy/2. $ ; spatial slice direction +else crpix1 = ppar.crpix1 +if ppar.crpix2 le 0. then $ + crpix2 = maximx/2. $ ; spatial slit direction +else crpix2 = ppar.crpix2 +; +; WCS keywords +sxaddpar,dhdr,'WCSDIM',2,' number of dimensions in WCS' +sxaddpar,dhdr,'WCSNAME','KCWI' +sxaddpar,dhdr,'EQUINOX',2000. +sxaddpar,dhdr,'RADESYS','FK5' +sxaddpar,dhdr,'CTYPE1','RA---TAN' +sxaddpar,dhdr,'CTYPE2','DEC--TAN' +sxaddpar,dhdr,'CUNIT1','deg',' RA units' +sxaddpar,dhdr,'CUNIT2','deg',' DEC units' +sxaddpar,dhdr,'CNAME1','KCWI RA',' RA name' +sxaddpar,dhdr,'CNAME2','KCWI DEC',' DEC name' +sxaddpar,dhdr,'CRVAL1',ra,' RA zeropoint' +sxaddpar,dhdr,'CRVAL2',dec,' DEC zeropoint' +sxaddpar,dhdr,'CRPIX1',crpix1,' RA reference pixel' +sxaddpar,dhdr,'CRPIX2',crpix2,' DEC reference pixel' +sxaddpar,dhdr,'CD1_1',cd11,' RA degrees per column pixel' +sxaddpar,dhdr,'CD2_1',cd21,' DEC degrees per column pixel' +sxaddpar,dhdr,'CD1_2',cd12,' RA degrees per row pixel' +sxaddpar,dhdr,'CD2_2',cd22,' DEC degrees per row pixel' +sxaddpar,dhdr,'LONPOLE',180.0,' Native longitude of Celestial pole' +sxaddpar,dhdr,'LATPOLE',0.0,' Celestial latitude of native pole' +sxaddpar,dhdr,'HISTORY',' '+kdgeom.progid+' '+systime(0,kdgeom.timestamp) +sxaddpar,dhdr,'HISTORY',' '+pre+' '+systime(0) +; +return +end diff --git a/pro/kcwi_apply_geom.pro b/pro/kcwi_apply_geom.pro new file mode 100644 index 0000000..adb4b91 --- /dev/null +++ b/pro/kcwi_apply_geom.pro @@ -0,0 +1,295 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_APPLY_GEOM +; +; PURPOSE: +; Apply the final geometric transformation to an input image. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_APPLY_GEOM,Img, Hdr, Kgeom, Ppar, Cube, CubeHdr +; +; INPUTS: +; Img - Stage1 processed object image to apply geometry to +; Hdr - the corresponding header of the input image +; Kgeom - KCWI_GEOM struct after KCWI_SOLVE_ARCS has been run +; +; INPUT KEYWORDS: +; WARPOUT - output warped image set +; +; OUTPUTS: +; Cube - 3-D data cube image [slice,x,wavelength] +; CubeHdr - A fits header with WCS info from Kgeom +; +; SIDE EFFECTS: +; None. +; +; PROCEDURE: +; Apply the geomtric transformation in Kgeom to Img and copy Hdr to +; CubeHdr and update WCS keywords. +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-JUL-31 Initial Revision +; 2013-AUG-07 Added padding of input image to extend output slices +; 2013-OCT-02 Re-ordered output cube axes and added WCS +; 2015-APR-25 Added CWI flexure hooks (MM) +; 2016-OCT-05 Removed CWI flexure routines +;- +pro kcwi_apply_geom,img,hdr,kgeom,ppar,cube,chdr,warpout=warpout +; +; startup +pre = 'KCWI_APPLY_GEOM' +q = '' +; +; Check structs +if kcwi_verify_geom(kgeom,/init) ne 0 then return +if kcwi_verify_ppar(ppar) ne 0 then begin + ppar = {kcwi_ppar} +endif +; +; get image original size +sz = size(img,/dim) +; +; get the padding for image +pad = kgeom.ypad +; +; get low trim wavelength +wave0 = kgeom.wave0out +; +; get high trim pixel +lastpix = long ( ( kgeom.wave1out - kgeom.wave0out ) / kgeom.dwout ) +; +; pad image to include full slices +pimg = fltarr(sz[0],sz[1]+pad+pad) +; +; insert original +y0 = pad +y1 = (sz[1]-1) + pad +pimg[*,y0:y1] = img +; +; get slice size +slsize = fix(5.*kgeom.refdelx) + 1 +; +; log values +kcwi_print_info,ppar,pre,'Slice dimensions (x,y)',slsize+1l,lastpix+1l, $ + format='(a,2i9)' +; +; image number +imgnum = sxpar(hdr,'FRAMENO') +; +; object name +object = sxpar(hdr,'TARGNAME',count=ntarg) +if ntarg le 0 then $ + object = sxpar(hdr,'OBJECT') +; +; image type +imgtyp = sxpar(hdr,'CALTYPE') +; +; mask +nasmask = (strpos(sxpar(hdr,'BNASNAM'),'Mask') ge 0) +; +; log +kcwi_print_info,ppar,pre,'Slicing and dicing image '+strn(imgnum)+': '+object+'...' +; +; loop over slices +for i=0,23 do begin + ; + ; coefficients + kwx = kgeom.kwx[*,*,i] + kwy = kgeom.kwy[*,*,i] + ; + ; perform the resampling here + warp = poly_2d(pimg,kwx,kwy,2,cubic=-0.5) + ; + ; output warped image if requested + if keyword_set(warpout) and $ + strcmp(strtrim(sxpar(hdr,'CALTYPE'),2),'cbars') eq 1 and $ + strcmp(strtrim(sxpar(hdr,'BUNIT'),2),'electrons') eq 1 then $ + mwrfits,warp,'warp'+string(i,form='(i02)')+'.fits',hdr + ; + ; check dimensions + wsz = size(warp,/dim) + ; + ; are we too big? + if i eq 0 and slsize ge wsz[0] then $ + kcwi_print_info,ppar,pre,'slice size gt xsize',slsize,wsz[0], $ + format='(a,2i13)',/warning + if i eq 0 and lastpix ge wsz[1] then $ + kcwi_print_info,ppar,pre,'lastpix gt ysize',lastpix,wsz[1], $ + format='(a,2i13)',/warning + ; + ; trim slice + slice = warp[0:slsize,0:lastpix<(wsz[1]-1)] + if i eq 0 then begin + sz = size(slice,/dim) + cube = fltarr(24,sz[0],sz[1]) + endif + cube[i,*,*] = slice + if ppar.verbose eq 1 then $ + print,strn(i)+' ',format='($,a)' +endfor +if ppar.verbose eq 1 then begin + print,'Done.',format='($,a)' + print,'' +endif +; +; update header +chdr = hdr +; +; image dimensions +sxaddpar,chdr,'NAXIS',3 +sxaddpar,chdr,'NAXIS1',24 +sxaddpar,chdr,'NAXIS2',sz[0] +sxaddpar,chdr,'NAXIS3',sz[1],' length of data axis 3',after='NAXIS2' +; +; spatial scale and zero point +sxaddpar,chdr,'BARSEP',kgeom.refdelx,' separation of bars (binned pix)' +sxaddpar,chdr,'BAR0',kgeom.x0out,' first bar pixel position' +; +; wavelength ranges +if nasmask then begin + sxaddpar,chdr, 'WAVALL0', kgeom.waveallm0, ' Low inclusive wavelength' + sxaddpar,chdr, 'WAVALL1', kgeom.waveallm1, ' High inclusive wavelength' + sxaddpar,chdr, 'WAVGOOD0',kgeom.wavegoodm0, ' Low good wavelength' + sxaddpar,chdr, 'WAVGOOD1',kgeom.wavegoodm1, ' High good wavelength' +endif else begin + sxaddpar,chdr, 'WAVALL0', kgeom.waveall0, ' Low inclusive wavelength' + sxaddpar,chdr, 'WAVALL1', kgeom.waveall1, ' High inclusive wavelength' + sxaddpar,chdr, 'WAVGOOD0',kgeom.wavegood0, ' Low good wavelength' + sxaddpar,chdr, 'WAVGOOD1',kgeom.wavegood1, ' High good wavelength' +endelse +sxaddpar,chdr, 'WAVMID',kgeom.wavemid, ' middle wavelength' +; +; wavelength solution RMS +sxaddpar,chdr,'AVWVSIG',kgeom.avewavesig,' Avg. bar wave sigma (Ang)' +sxaddpar,chdr,'SDWVSIG',kgeom.stdevwavesig,' Stdev. bar wave sigma (Ang)' +; +; geometry solution RMS +xmo = moment(kgeom.xrsd,/nan) +ymo = moment(kgeom.yrsd,/nan) +sxaddpar,chdr, 'GEOXGSG', xmo[0], ' Global geometry X sigma (pix)' +sxaddpar,chdr, 'GEOYGSG', ymo[0], ' Global geometry Y sigma (pix)' +for i=0,23 do begin + sxaddpar,chdr, 'GEOXSG'+strn(i), kgeom.xrsd[i],' Avg. geometry X sigma (pix) slice '+strn(i),format='F7.3' + sxaddpar,chdr, 'GEOYSG'+strn(i), kgeom.yrsd[i],' Avg. geometry Y sigma (pix) slice '+strn(i),format='F7.3' +endfor +; +; pixel scales +sxaddpar,chdr,'PXSCL', kgeom.pxscl*kgeom.xbinsize,' Pixel scale along slice' +sxaddpar,chdr,'SLSCL', kgeom.slscl,' Pixel scale purpendicular to slices' +; +; geometry origins +sxaddpar,chdr, 'CBARSFL', kgeom.cbarsfname,' Continuum bars image' +sxaddpar,chdr, 'ARCFL', kgeom.arcfname, ' Arc image' +sxaddpar,chdr, 'CBARSNO', kgeom.cbarsimgnum,' Continuum bars image number' +sxaddpar,chdr, 'ARCNO', kgeom.arcimgnum, ' Arc image number' +sxaddpar,chdr, 'GEOMFL', kgeom.geomfile,' Geometry file' +; +; get sky coords +rastr = sxpar(hdr,'RA',count=nra) +decstr = sxpar(hdr,'DEC',count=ndec) +if nra ne 1 or ndec ne 1 then begin + rastr = sxpar(hdr,'TARGRA',count=nra) + decstr = sxpar(hdr,'TARGDEC',count=ndec) +endif +if nra eq 1 and ndec eq 1 then begin + radec_parse,rastr,decstr,':',ra,dec +endif else begin + ra = -99.d0 + dec = -99.d0 +endelse +; +; Position Angle ( = SKYPA) in degrees +; Plus an offset between rotator and IFU (may be zero) +skypa = sxpar(hdr,'ROTPOSN',count=npa) + sxpar(hdr,'ROTREFAN') +crota = -(skypa + kgeom.rotoff) / !RADEG +sxaddpar,chdr,'IFUPA',skypa,' IFU position angle (degrees)' +sxaddpar,chdr,'IFUROFF',kgeom.rotoff,' IFU-SKYPA offset (degrees)' +; +; pixel scales +cdelt1 = -kgeom.slscl ; RA degrees per px (column) +cdelt2 = kgeom.pxscl*kgeom.xbinsize ; Dec degrees per slice (row) +; +; did we get good coords? +if nra ne 1 or ndec ne 1 or npa ne 1 then begin + ; + ; no good coords + ; a warning for objects + if strcmp(imgtyp,'object') eq 1 then begin + kcwi_print_info,ppar,pre,'no coords for image',imgnum,imgtyp, $ + format='(a,2x,a,2x,a)',/warning + ; otherwise just info + endif else begin + kcwi_print_info,ppar,pre,'no coords for image',imgnum,imgtyp, $ + format='(a,2x,a,2x,a)' + endelse + ; + ; zero coords + ra = 0. + dec = 0. + ; + ; nominal CD matrix (no rotation) + CD11 = cdelt1*cos(0.) + CD12 = abs(cdelt2)*sign(cdelt1)*sin(0.) + CD21 = -abs(cdelt1)*sign(cdelt2)*sin(0.) + CD22 = cdelt2*cos(0.) +endif else begin + ; + ; calculate CD matrix + CD11 = cdelt1*cos(crota) ; RA degrees per column + CD12 = abs(cdelt2)*sign(cdelt1)*sin(crota) ; RA degress per row + CD21 = -abs(cdelt1)*sign(cdelt2)*sin(crota) ; DEC degress per column + CD22 = cdelt2*cos(crota) ; DEC degrees per row +endelse +; +; get reference pixels +if ppar.crpix1 le 0. then $ + crpix1 = 12. $ ; spatial slice direction: 24/2 +else crpix1 = ppar.crpix1 +if ppar.crpix2 le 0. then $ + crpix2 = sz[0]/2. $ ; spatial slit direction +else crpix2 = ppar.crpix2 +if ppar.crpix3 le 0. then $ + crpix3 = 1. $ ; wavelength direction +else crpix3 = ppar.crpix3 +; +; WCS keywords +sxaddpar,chdr,'WCSDIM',3,' number of dimensions in WCS' +sxaddpar,chdr,'WCSNAME','KCWI' +sxaddpar,chdr,'EQUINOX',2000. +sxaddpar,chdr,'RADESYS','FK5' +sxaddpar,chdr,'CTYPE1','RA---TAN' +sxaddpar,chdr,'CTYPE2','DEC--TAN' +sxaddpar,chdr,'CTYPE3','AWAV',' Air Wavelengths' +sxaddpar,chdr,'CUNIT1','deg',' RA units' +sxaddpar,chdr,'CUNIT2','deg',' DEC units' +sxaddpar,chdr,'CUNIT3','Angstrom',' Wavelength units' +sxaddpar,chdr,'CNAME1','KCWI RA',' RA name' +sxaddpar,chdr,'CNAME2','KCWI DEC',' DEC name' +sxaddpar,chdr,'CNAME3','KCWI Wavelength',' Wavelength name' +sxaddpar,chdr,'CRVAL1',ra,' RA zeropoint' +sxaddpar,chdr,'CRVAL2',dec,' DEC zeropoint' +sxaddpar,chdr,'CRVAL3',wave0,' Wavelength zeropoint' +sxaddpar,chdr,'CRPIX1',crpix1,' RA reference pixel' +sxaddpar,chdr,'CRPIX2',crpix2,' DEC reference pixel' +sxaddpar,chdr,'CRPIX3',crpix3,' Wavelength reference pixel' +sxaddpar,chdr,'CD1_1',cd11,' RA degrees per column pixel' +sxaddpar,chdr,'CD2_1',cd21,' DEC degrees per column pixel' +sxaddpar,chdr,'CD1_2',cd12,' RA degrees per row pixel' +sxaddpar,chdr,'CD2_2',cd22,' DEC degrees per row pixel' +sxaddpar,chdr,'CD3_3',kgeom.dwout,' Wavelength Angstroms per pixel' +sxaddpar,chdr,'LONPOLE',180.0,' Native longitude of Celestial pole' +sxaddpar,chdr,'LATPOLE',0.0,' Celestial latitude of native pole' +sxaddpar,chdr,'HISTORY',' '+kgeom.progid+' '+systime(0,kgeom.timestamp) +sxaddpar,chdr,'HISTORY',' '+pre+' '+systime(0) +; +return +end diff --git a/pro/kcwi_central_wave.pro b/pro/kcwi_central_wave.pro new file mode 100644 index 0000000..7f599fe --- /dev/null +++ b/pro/kcwi_central_wave.pro @@ -0,0 +1,107 @@ +; +; Copyright (c) 2016, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_CENTRAL_WAVE +; +; PURPOSE: +; Given the grating installed in KCWI and camera articulation +; angle and grating angle, return the approximate central +; wavelength falling onto the PCWI CCD. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; Result = KCWI_CENTRAL_WAVE(GratingNum, GratingAngle, CameraAngle) +; +; INPUTS: +; GratingNum - KCWI grating number (1 - 5) +; GratingAngle - Grating angle in degrees (0 - 360) +; CameraAngle - Camera articulation angle (-5 - 105) +; +; KEYWORDS: +; pwave - the peak wavelength in angstroms +; +; RETURNS: +; Approximate central wavelength on the PCWI CCD. Returns a negative +; value if the program encountered an error. +; +; SIDE EFFECTS: +; None. +; +; PROCEDURE: +; Uses formulae provided by Mat Matuszewski to calculate central and +; peak wavelengths. +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Matt Matuszewski (matmat@caltech.edu) +; 2013-AUG-12 Initial version +; 2014-AUG-14 Added Yellow grating +; 2014-SEP-10 Added MEDREZ (Richardson) grating +; 2016-JUN-30 Kludged from pcwi_central_wave +;- +function kcwi_central_wave, gid, grangle, camang, pwave=pwave, $ + help=help + +pre = 'KCWI_CENTRAL_WAVE' +wave = -1. + +if n_params(0) lt 3 or keyword_set(help) then begin + print,pre+': Info - Usage: Wave = '+pre+'( GratNum, GratAngle, CamAngle)' + print,pre+': Info - GratNum is 1-BH3, 2-BL, 3-BH2, 4-BM, 5-BH1' + print,pre+': Info - GratAngle is the grating angle in degrees' + print,pre+': Info - CamAngle is the camera articulation angle in degrees' + print,pre+': Info - wavelength returned in Angstroms' + return,wave +endif + +case gid of + 1: begin ; BH3 + corang = 180.0 + rho = 2.8017 + d0 = 7.714e-2 + d1 = 1.344e-4 + end + 2: begin ; BL + corang = 0.0 + rho = 0.870 + d0 = 9.392e-2 + d1 = 3.896e-5 + end + 3: begin ; BH2 + corang = 180.0 + rho = 3.2805 + d0 = 8.721e-2 + d1 = 1.560e-4 + end + 4: begin ; BM + corang = 0.0 + rho = 1.901 + d0 = 8.315e-2 + d1 = 9.123e-5 + end + 5: begin ; BH1 + corang = 180.0 + rho = 3.3 + d0 = 1.e-2 + d1 = 1.e-4 + end + ELSE: begin + print,pre+': Error - Unrecognized grating requested: '+ $ + strn(gid)+'. Returing -1.' + return,wave + end; else +endcase + +alpha = grangle - (corang + 13.0) +beta = camang - alpha + +pwave = ( sin( alpha*!DTOR ) - d0 ) / d1 +cwave = ( ( sin( alpha*!DTOR ) + sin( beta*!DTOR ) ) / rho ) * 10000. + +return,cwave +end diff --git a/pro/kcwi_cfg__define.pro b/pro/kcwi_cfg__define.pro new file mode 100644 index 0000000..aa701f3 --- /dev/null +++ b/pro/kcwi_cfg__define.pro @@ -0,0 +1,194 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_CFG__DEFINE +; +; PURPOSE: +; This procedure defines the structure for all the observation +; information for KCWI raw data read from the raw headers. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; kcfg = {KCWI_CFG} +; +; INPUTS: +; None +; +; OUTPUTS: +; None +; +; SIDE EFFECTS: +; defines the KCWI_CFG data structure +; +; PROCEDURE: +; Provide the automatic structure definition for the KCWI raw header +; observation structure KCWI_CFG. +; +; EXAMPLE: +; Instantiate and initialize a single cfg struct for KCWI: +; +; kcfg = {kcwi_cfg} +; kcfg = struct_init(kcfg) +; +; NOTES: +; Keep structure name on a line by itself (see struct_init.pro). +; Keep tag names 15 chars or less. +; All tags above the DRIVED PROPERTIES should exactly match header +; keywords (except case). +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-MAY-03 Initial version +;- +pro kcwi_cfg__define +; +tmp = { kcwi_cfg, $ +; +; Observation properties + observer:'', $ ; observer + telescop:'', $ ; telescope + instrume:'', $ ; instrument + object:'', $ ; object name + targname:'', $ ; target name + datepclr:'', $ ; UT start of observation (YYYY-MM-DDTHH:MM:SS) + daterend:'', $ ; UT end of observation (YYYY-MM-DDTHH:MM:SS) + el:-99.d0, $ ; Telescope elevation (deg) + az:-99.d0, $ ; Telescope azimuth (deg) + parang:-99.d0, $ ; Parallactic angle (deg) + epoch:-99.d0, $ ; Coordinate epoch + rotposn:-999.d0, $ ; Rotator user angle (deg) + rotrefan:-999.d0, $ ; Rotator reference angle (deg) + rotmode:'', $ ; Rotator mode + xposure:-9.0, $ ; shutter open duration (s) + telapse:-9.0, $ ; dark current duration (s) + airmass:-1.0, $ ; Airmass + ofname:'', $ ; Original file name + +; +; Configuration properties + hatpos:'', $ ; hatch position: Open or Closed + flimagin:'', $ ; dome lamp imaging mode: on or off + flspectr:'', $ ; dome lamp spectra mode: on or off + caltype:'', $ ; calibration type: bias, dark, arc, etc. + frameno:0l, $ ; image number + skyobs:0, $ ; sky observation? 0 - object, 1 - sky + shuffmod:0, $ ; is this a Nod & Shuffle observation? + bgratnam:'', $ ; Blue grating id + bgratnum:0, $ ; Blue graing number + bgrangle:0., $ ; Blue grating angle (degrees) + bgrenc:0l, $ ; Blue grating rotator encoder steps + bfiltnam:'', $ ; Blue filter id + bfiltnum:0, $ ; Blue filter number + bartang:0., $ ; Blue articulation angle + bartenc:0l, $ ; Blue cmaera articulation encoder steps + bcwave:0., $ ; Blue central wavelength (Ang) + bpwave:0., $ ; Blue peak wavelength (Ang) + bfocpos:0l, $ ; Blue focus stage encoder steps + bfocus:0., $ ; Blue focus (mm) + bnasnam:'', $ ; Blue mask position name + bnaspos:0, $ ; Blue mask position + shufrows:0, $ ; Number of CCD rows shuffled + calpnam:'', $ ; Cal position name ('Sky', 'Polar', 'Lens') + callang:0., $ ; Cal polarizer angle + ifunum:0, $ ; Slicer number (0-5, -1=unknown) + ifunam:'', $ ; Slicer name ("Small", etc.) + cwave:0., $ ; central wavelength (Ang) + gratanom:0., $ ; grating angle anomoly (degrees) + wave0:0., $ ; blue end of wavelength range (Ang) + wave1:0., $ ; red end of wavelength range (Ang) + dwav:0., $ ; average dispersion (Ang/pix) +; +; Temperatures + tmpa1:0., $ ; Blue CCD temperature (K) + tmpa8:0., $ ; Optical bench temperature (K) +; +; Cal Lamp properties + lmp0nam:'', $ ; Lamp 0 name + lmp0stat:0, $ ; Lamp 0 status (0 - off, 1 - on) + lmp0shst:0, $ ; Lamp 0 shutter (0 - closed, 1 - open) + lmp1nam:'', $ ; Lamp 0 name + lmp1stat:0, $ ; Lamp 0 status (0 - off, 1 - on) + lmp1shst:0, $ ; Lamp 0 shutter (0 - closed, 1 - open) + lmp3nam:'', $ ; Lamp 0 name + lmp3stat:0, $ ; Lamp 0 status (0 - off, 1 - on) +; +; CCD properties + biasrn1:3., $ ; bias read noise for amp 1 in electrons/pixel + biasrn2:3., $ ; bias read noise for amp 2 in electrons/pixel + biasrn3:3., $ ; bias read noise for amp 3 in electrons/pixel + biasrn4:3., $ ; bias read noise for amp 4 in electrons/pixel + gain1:0.145, $ ; gain for amp 1 in electrons/DN + gain2:0.145, $ ; gain for amp 2 in electrons/DN + gain3:0.145, $ ; gain for amp 3 in electrons/DN + gain4:0.145, $ ; gain for amp 4 in electrons/DN + gainmul:10, $ ; gain multiplier (1, 2, 5, 10) + nampsxy:'1 1', $ ; number of amplifiers in x and y + ampmode:'LL', $ ; amplifier mode: LL,LR,UL,UR,DUP,DLO,QUAD + nvidinp:0, $ ; number of amps (video inputs) + ccdmode:0, $ ; ccd mode: 0 - slow, 1 - fast +; +; Image geometry properties + naxis:0, $ ; number of data axes + naxis1:0, $ ; length of data axis 1 + naxis2:0, $ ; length of data axis 2 +; +; Wavelength values for reduced data + crval3:0., $ ; wavelength zeropoint + crpix3:0., $ ; wavelength reference pixel + cdelt3:0., $ ; wavelength dispersion Ang/px +; +; DERIVED PROPERTIES (above should match real FITS header keywords) + ra:-99.d0, $ ; RA + dec:-99.d0, $ ; Dec + juliandate:0.d0, $ ; Julian date of observation + date:'', $ ; UT start of observation (YYYY-MM-DDTHH:MM:SS) + binning:0, $ ; is image binned? + xbinsize:1, $ ; binning in x + ybinsize:1, $ ; binning in y + exptime:-9., $ ; exposure time (either xposure or telapse) +; +; Nod-and-shuffle rows (in trimmed image) + nsskyr0:0, $ ; Sky region row 0 (bottom, pix) + nsskyr1:0, $ ; Sky region row 1 (top, pix) + nsobjr0:0, $ ; Object region row 0 (bottom, pix) + nsobjr1:0, $ ; Object region row 1 (top, pix) +; +; Input file properties + obsfname:'', $ ; input observation FITS file name (sans dir) + obsdir:'', $ ; input directory for observation + obstype:'', $ ; observation type: 'zero', 'cal', 'obj', 'std' + imgnum:0L, $ ; image number +; +; Instrument state + imgtype:'', $ ; observation type: bias, dark, arc, etc. + gratid:'', $ ; Grating id name + gratnum:-1, $ ; Gratind id number + grangle:0., $ ; Grating angle (degrees) + grenc:0l, $ ; Grating rotator encoder steps + filter:'', $ ; Filter id name + filtnum:-1, $ ; filter id number + camang:0., $ ; Camera articulation angle (degrees) + campos:0l, $ ; Camera articulation encoder steps + focpos:0l, $ ; Camera focus encoder steps + focus:0., $ ; Camera focus (mm) + nasmask:0, $ ; is the Nod & Shuffle mask deployed? + lampname:'', $ ; illumination source +; +; Master group properties + groupnum:-1l, $ ; group image number + nimages:0, $ ; number of images in group + grouplist:'', $ ; range list of images in group + groupfile:'', $ ; filename of grouped image + grouppar:'', $ ; parameter file for grouped image +; +; Status + initialized:0, $ ; 0 - no, 1 - yes +; +; timestamp + timestamp:0.d0 $ ; timestamp for struct (unix seconds) + } +end diff --git a/pro/kcwi_cfg_string.pro b/pro/kcwi_cfg_string.pro new file mode 100644 index 0000000..f54097e --- /dev/null +++ b/pro/kcwi_cfg_string.pro @@ -0,0 +1,39 @@ +function kcwi_cfg_string,kcfg,long=long,delim=delim,bias=bias + nc = n_elements(kcfg) + cstr = strarr(nc) + if keyword_set(delim) then $ + del = ':' $ + else del = '' + for i=0,nc-1 do begin + cstr[i] = strn(kcfg[i].xbinsize)+strn(kcfg[i].ybinsize) + del + if keyword_set(bias) then begin + if keyword_set(delim) then $ + cstr[i] += strtrim(kcfg[i].ampmode,2) + del $ + else cstr[i] += string(strtrim(kcfg[i].ampmode,2), $ + format='(a-3)') + if keyword_set(delim) then $ + cstr[i] += string(kcfg[i].gainmul,form='(i2)') $ + + del $ + else cstr[i] += string(kcfg[i].gainmul,form='(i2)') + cstr[i] += strn(kcfg[i].ccdmode) + endif else begin + if keyword_set(delim) then $ + cstr[i] += string(strtrim(kcfg[i].bgratnam,2), $ + format='(a-3)')+ del $ + else cstr[i] += string(strtrim(kcfg[i].bgratnam,2), $ + format='(a-4)') + cstr[i] += strmid(kcfg[i].ifunam,0,1) + del + if keyword_set(long) then begin + cstr[i] += strtrim(kcfg[i].calpnam,2) + del + if keyword_set(delim) then $ + udel = del $ + else udel = ' ' + cstr[i] += strtrim(string(kcfg[i].callang, $ + format='(f9.1)'),2)+udel + endif + cstr[i] += strtrim(string(kcfg[i].bcwave, $ + format='(f10.1)'),2) + endelse + endfor + return, cstr +end diff --git a/pro/kcwi_dgeom__define.pro b/pro/kcwi_dgeom__define.pro new file mode 100644 index 0000000..17a8f8d --- /dev/null +++ b/pro/kcwi_dgeom__define.pro @@ -0,0 +1,130 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_DGEOM__DEFINE +; +; PURPOSE: +; This procedure defines the structure for all the information for +; KCWI direct mode geometric transformations. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; kdgeom = {KCWI_DGEOM} +; +; INPUTS: +; None +; +; OUTPUTS: +; None +; +; SIDE EFFECTS: +; defines the KCWI_DGEOM data structure +; +; PROCEDURE: +; Provide the automatic structure definition for the KCWI direct mode +; geometric transformation structure KCWI_DGEOM. +; +; EXAMPLE: +; Instantiate and initialize a single dgeom struct for KCWI: +; +; kdgeom = {kcwi_dgeom} +; kdgeom = struct_init(kdgeom) +; +; NOTES: +; Keep structure name on a line by itself (see struct_init.pro). +; Keep tag names 15 chars or less. +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2016-NOV-09 Initial version +;- +pro kcwi_dgeom__define +; +tmp = { kcwi_dgeom, $ +; +; output filename + geomfile:'', $ ; output save file for dgeom struct +; +; calibration images + cbarsimgnum:0l, $ ; cbars image number (0 - 9999) + cbarsjd:0.d0, $ ; cbars image julian date + cbarsfname:'', $ ; cbars image file name + arcimgnum:0l, $ ; arc image number + arcjd:0.d0, $ ; arc image julian date + arcfname:'', $ ; arc image file name +; +; IFU properties + ifunum: 0, $ ; Slicer number (1-3, Large, Medium, Small) + ifunam: '', $ ; Slicer name ("Large", etc.) +; +; filter properties + filter:'', $ ; filter id + filtnum:-1, $ ; filter number +; +; encoder positions + campos:0L, $ ; Camera encoder position + camang:0.0d, $ ; Camera articulation angle in degrees +; +; pixel scale + pxscl:0.00008096d0, $ ; degrees per unbinned spatial pixel +; +; slice scale + slscl:0.00075437d0, $ ; degrees per slice +; +; IFU rotator offset + rotoff:0.0, $ ; degrees from rotator PA +; +; CCD binning + xbinsize:1, $ ; binning in x (pixels) + ybinsize:1, $ ; binning in y (pixels) +; +; CCD size + nx:0, $ ; number of x pixels in image + ny:0, $ ; number of y pixels in image +; +; CCD readnoise + rdnoise:3., $ ; e-/pixel +; +; minimum pixels for arc finding + minpix:6, $ ; minimum number of illuminated pix in arc col +; +; slice angles + angles:fltarr(24), $ ; angle in degrees of each slice image +; +; ref slice for x offsets + refslice:11, $ ; reference slice (0 - 23, default is 11) +; +; x offsets + xoff:fltarr(24), $ ; x offset as determined from bars image +; +; starting x position of arc slice image + x0:fltarr(24), $ ; starting x pos in image of each slice (pixels) +; +; ending x position of arc slice image + x1:fltarr(24), $ ; ending x pos in image of each slice (pixels) +; +; starting y position of arc slice image + y0:fltarr(24), $ ; starting y pos in image of each slice (pixels) +; +; ending y position of arc slice image + y1:fltarr(24), $ ; ending y pos in image of each slice (pixels) +; +; width of arc slice image + wid:fltarr(24), $ ; width of arc slice image in pixels +; +; do we fit with a gaussian (or use double erf?) + do_gauss:0, $ ; 0 - no, 1 - yes +; +; Status + initialized:0, $ ; 0 - no, 1 - yes + status:-1, $ ; 0 - good fit, else not good + progid:'', $ ; program that last modified the dgeom file +; +; timestamp + timestamp:0.d0 $ ; timestamp for struct (unix seconds) + } +end diff --git a/pro/kcwi_extract_arcs.pro b/pro/kcwi_extract_arcs.pro new file mode 100644 index 0000000..66f1e13 --- /dev/null +++ b/pro/kcwi_extract_arcs.pro @@ -0,0 +1,217 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_EXTRACT_ARCS +; +; PURPOSE: +; Extract individual arc spectra along continuum bar image bars. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_EXTRACT_ARCS,ArcImg, Kgeom, Spec +; +; INPUTS: +; ArcImg - 2-D image of type 'arc' +; Kgeom - KCWI_GEOM struct constructed from KCWI_TRACE_CBARS run +; +; INPUT KEYWORDS: +; NAVG - number of pixels in x to average for each y pixel +; CCWINDIW- cross-correlation window in pixels for offset plots +; VERBOSE - extra output +; DISPLAY - set to display a plot of each spectrum after shifting +; +; OUTPUTS: +; Spec - a array of vectors, one for each bar defined by Kgeom +; +; SIDE EFFECTS: +; None. +; +; PROCEDURE: +; Use the Kgeom.[kx,ky] coefficients to de-warp the arc image, after +; which a spectrum is extracted by simple summing for each of the +; bars in the 'cbars' image using the positions from Kgeom.barx. The +; spectra are assumed aligned with dispersion along the y-axis of +; the image. The extraction window in x at each y pixel is determined +; by the keyword NAVG. Since the spectra are extracted using whole +; pixel windows, the reference x position for each spectrum is then +; centered on the window center. These are recorded in Kgeom.refx. +; +; EXAMPLE: +; Define the geometry from a 'cbars' image and use it to extract and +; display the spectra from an 'arc' image from the same calibration +; sequence. +; +; cbars = mrdfits('image7142_int.fits',0,chdr) +; kcwi_trace_cbars,cbars,Kgeom,/centroid +; arc = mrdfits('image7140_int.fits',0,ahdr) +; kcwi_extract_arcs,arc,kgeom,arcspec,ppar +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-JUL-23 Initial Revision +; 2013-JUL-24 Checks if NASMASK in place, computes bar offsets +; 2013-DEC-09 Default nav now 2, takes median of bar sample +;- +; +pro kcwi_extract_arcs,img,kgeom,spec,ppar, $ + navg=navg, ccwindow=ccwindow, $ + help=help +; +; startup +pre = 'KCWI_EXTRACT_ARCS' +q = '' +; +; check inputs +if n_params(0) lt 3 or keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', ArcImg, Kgeom, Spec' + return +endif +; +; check Ppar +if kcwi_verify_ppar(ppar,/init) ne 0 then return +; +; Check Kgeom +if kcwi_verify_geom(kgeom) ne 0 then return +; +; keywords +if keyword_set(navg) then $ + nav = (fix(navg/2)-1) > 1 $ +else nav = 8/kgeom.ybinsize +; +if keyword_set(ccwindow) then $ + ccwn = ccwindow $ +else ccwn = kgeom.ccwn +do_plots = ppar.display +; +; check size of input image +sz = size(img,/dim) +if sz[0] eq 0 then begin + kcwi_print_info,ppar,pre, $ + 'input image must be 2-D, preferrably an arc image.',/error + return +endif +nx = sz[0] +ny = sz[1] +; +; number of spectra +ns = n_elements(kgeom.barx) +; +; log +kcwi_print_info,ppar,pre,'Extracting arcs from image',kgeom.arcimgnum +; +; transform image +warp = poly_2d(img,kgeom.kx,kgeom.ky,2,cubic=-0.5) +; +; output spectra +spec = fltarr(ny,ns) +; +; bar offsets +boff = fltarr(ns) +; +; loop over x values, get sample, compute spectrum +for i=0,ns-1 do begin + xpix = fix(kgeom.barx[i]) + kgeom.refx[i] = float(xpix)+0.5 + ; + ; get sample + sub = warp[(xpix-nav):(xpix+nav),*] + ; + ; compute spectrum + vec = median(sub,dim=1) + ims_asym,vec,mn,sig,wgt,siglim=[1.5,3.0] + good = where(wgt eq 1,ngood) + if ngood gt 9 then $ + cont = min(vec[good]) $ + else cont = median(vec[good]) + cont = mn - 3.*sig + vec -= cont + spec(*,i) = vec +endfor +; +; get arc header +hdr = headfits(kgeom.arcfname) +sxaddpar,hdr,'HISTORY',' '+pre+' '+systime(0) +; +; geometry origins +sxaddpar,hdr, 'CBARSFL', kgeom.cbarsfname,' Continuum bars image' +sxaddpar,hdr, 'ARCFL', kgeom.arcfname, ' Arc image' +sxaddpar,hdr, 'CBARSNO', kgeom.cbarsimgnum,' Continuum bars image number' +sxaddpar,hdr, 'ARCNO', kgeom.arcimgnum, ' Arc image number' +sxaddpar,hdr, 'GEOMFL', kgeom.geomfile,' Geometry file' +; +; write out arcs +outfile = kcwi_get_imname(ppar,kgeom.arcimgnum,"_arcs",/reduced) +kcwi_print_info,ppar,pre,'Writing',outfile,/info,format='(a,1x,a)' +mwrfits,spec,outfile,hdr,/create,/iscale +; +; setup for display (if requested) +if do_plots ge 2 then begin + if !d.window lt 0 then $ + window,0,title=kcwi_drp_version() $ + else wset,0 + deepcolor + !p.background=colordex('white') + !p.color=colordex('black') + th=2 + si=1.75 + w=findgen(ny) + xrng=[-ccwn,ny+ccwn] + if kgeom.nasmask eq 1 then $ + xrng = [((1200./kgeom.ybinsize)-ccwn), $ + ((3000./kgeom.ybinsize)+ccwn)] +endif +; +; get reference bar from kgeom +if kgeom.refbar ge 0 and kgeom.refbar lt 120 then $ + refb = kgeom.refbar $ +else refb = 57 ; default +; +; log +kcwi_print_info,ppar,pre,'cross-correlating with reference bar',refb +kcwi_print_info,ppar,pre,'using cc window (px)',ccwn +; +; cross-correlate to reference line +for i=0,ns-1 do begin + off = 0. + if i ne refb then $ + off = ccpeak(spec[*,i],spec[*,refb],ccwn,offset=kgeom.ccoff[i/5]) + ; + ; record offsets + boff[i] = off + ; + ; log + kcwi_print_info,ppar,pre,'Bar#, yoffset',i,off,format='(a,i4,f9.3)',info=2 + ; + ; plot if requested + if do_plots ge 2 then begin + plot,w,spec[*,refb],thick=th,xthick=th,ythick=th, $ + charsi=si,charthi=th, $ + xran=xrng,xtitle='Pixel',xstyle=1, $ + ytitle='Int', /nodata, $ + title='Image: '+strn(kgeom.arcimgnum)+ $ + ', Bar: '+strn(i)+', Slice: '+strn(fix(i/5))+$ + ', Y Offset: '+strtrim(string(off,form='(f9.3)'),2)+' px' + oplot,w,spec[*,refb],color=colordex('R'),thick=th + oplot,w+off,spec[*,i],color=colordex('G') + kcwi_legend,['RefBar: '+strn(refb),'Bar: '+strn(i)], $ + linesty=[0,0], thick=[th,0], $ + color=[colordex('R'),colordex('G')],box=0 + read,'Next? (Q-quit plotting, -next): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then $ + do_plots = 0 + endif +endfor +; +; update Kgeom +kgeom.baroff = boff +; +; Kgeom timestamp +kgeom.progid = pre +kgeom.timestamp = systime(1) +; +return +end diff --git a/pro/kcwi_fit_center.pro b/pro/kcwi_fit_center.pro new file mode 100644 index 0000000..6543e31 --- /dev/null +++ b/pro/kcwi_fit_center.pro @@ -0,0 +1,504 @@ +; +; Copyright (c) 2014, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_FIT_CENTER +; +; PURPOSE: +; Solve the wavelength solutions for the central third of +; each bar of the arc spectrum using a fourth degree +; polynomial approximation +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_FIT_CENTER, Specs, Kgeom, Ppar, Centcoeff +; +; INPUTS: +; Specs - a array of arc spectra produced by KCWI_EXTRACT_ARCS +; Kgeom - KCWI_GEOM struct from KCWI_TRACE_CBARS and KCWI_EXTRACT_ARCS +; Ppar - KCWI_PPAR pipeline parameter struct +; +; OUTPUTS: +; Centcoeff - coefficients of central fit (4th degree) +; +; INPUT KEYWORDS: +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Matt Matuszewski +; 2014-SEP-18 Initial Revision +; 2015-APR-23 JDN: added cosine bell taper to minimize edge effects +;- +; +pro kcwi_fit_center, specs, kgeom, ppar, centcoeff + +pre = 'KCWI_FIT_CENTER' +q='' +; +; log info +kcwi_print_info,ppar,pre,systime(0) +; +; do we want to display stuff? +display = (ppar.display ge 2) +ddisplay = (ppar.display ge 3) + +if ppar.display ge 2 then begin + if !d.window lt 0 then $ + window,0,title=kcwi_drp_version() $ + else wset,0 + deepcolor + !p.multi=0 + !p.background=colordex('white') + !p.color=colordex('black') + th=2.0 + si=2.0 +endif +; +; set some system parameters -- these may need to be populated +; differently later +pix = 0.0150d ; pixel size in mm +ybin = kgeom.ybinsize ; binning in spectral direction +fcam = 305.0d ; focal length of camera in mm +gamma = 4.0d ; mean out-of-plane angle for diffraction. +; +; which image number +imgnum = kgeom.arcimgnum +imglab = 'Img # '+strn(imgnum) +; +; which slicer? +ifunum = kgeom.ifunum +ifunam = kgeom.ifunam +; +; which grating? +grating = kgeom.gratid +; +; which filter? +filter = kgeom.filter +; +; is the N&S mask in? +nasmask = kgeom.nasmask +; +; central wavelength? +cwvl = kgeom.cwave +; +; any anomolous tilt in the grating? +gratanom = kgeom.gratanom +; +; set the grating specific parameters +rho = kgeom.rho +; +; which is the reference bar? +refbar = kgeom.refbar +; +; we will be using a third degree fit here +degree = 4 +kcwi_print_info,ppar,pre,'Using polynomial approximation of degree',degree, $ + format='(a,i5)' +; +; report taper fraction +kcwi_print_info,ppar,pre, $ + 'Using cross-correlation bell cosine taper fraction of',ppar.taperfrac,$ + format='(a,f9.3)' +; +; load the reference atlas spectrum. +kcwi_read_atlas,kgeom,ppar,refspec,refwvl,refdisp +; +; make sure spectrum is double precision +specs = double(specs) +; +; get dimensions +specsz = size(specs,/dim) +; +; coefficients for central region fit +centcoeff = dblarr(9,120) +; +; Next we refine the central dispersion estimate +; +; 0- compute alpha +prelim_alpha = kgeom.grangle - 13.0d - kgeom.adjang +; +; 1- compute the prelim angle of diffraction +;prelim_beta = asin(cwvl/10000.0 * rho/2.0)+slant/!radeg +prelim_beta = kgeom.camang - prelim_alpha +; +; 1b - add the grating tilt anomoly +;prelim_beta = prelim_beta + gratanom/!radeg +; +; 2- compute the preliminary dispersion +prelim_disp = cos(prelim_beta/!radeg)/rho/fcam*(pix*ybin)*1e4 +; the 1e4 is there to make the units Angstrom/binnedpixel +; +; need to correct for the out-of-band angle here... not much, but +; there is some... so +; +prelim_disp *= cos(gamma/!radeg) +; +; report results +kcwi_print_info,ppar,pre,'Initial alpha, beta',prelim_alpha, prelim_beta, $ + format='(a,2f9.2)' +kcwi_print_info,ppar,pre,'Initial calculated dispersion (A/binned pixel)', $ + prelim_disp,format='(a,f8.3)' +; +; 3- generate an index array with its 0 point at the center of the +; detector +x0 = specsz[0]/2 +xvals = dindgen(specsz[0])-x0 +; +; 4- Pick out the central third of the detector in the spectral +; direction. This is where the dispersion is linear and is the +; best place to try to cross-correlate the reference bar spectrum +; with the ThAr spectrum to refine the shift. +; +; the min and max row to use for the adjustment +minrow = (1*specsz[0])/3 +maxrow = (2*specsz[0])/3 +; +; Unless we are working on the BL grating. In this case pick the +; central 3 fifths +if strpos(grating,'BL') ge 0 then begin + minrow = (1*specsz[0])/5 + maxrow = (4*specsz[0])/5 +endif +; +; the corresponding preliminary wavelengths +prelim_wvl = cwvl + xvals*prelim_disp +prelim_minwvl = min( [prelim_wvl[minrow],prelim_wvl[maxrow]] ) +prelim_maxwvl = max( [prelim_wvl[minrow],prelim_wvl[maxrow]] ) +; +; now we have to interpolate the bar spectrum to the same scale as the +; atlas spectrum. +; +; subspectrum to interpolate and subindex to interpolate from +prelim_spec = reform(specs[minrow:maxrow,refbar]) +prelim_xvals = xvals[minrow:maxrow] +prelim_subwvl = prelim_wvl[minrow:maxrow] +; +; determine the wavelengths to interpolate to and extract the relevant +; atlas portion +qwvl = where(refwvl gt prelim_minwvl and refwvl lt prelim_maxwvl, nqwvl) +if nqwvl eq 0 then begin + kcwi_print_info,ppar,pre,'Did not find atlas data to match to',/error + kgeom.status=2 + return +endif +; +prelim_refspec = refspec[qwvl] +prelim_refwvl = refwvl[qwvl] +; +; and interpolate +prelim_intspec = interpol(prelim_spec,prelim_subwvl,prelim_refwvl,/spline) +; +; check for scattered light problems +mmm,prelim_intspec,skymod,skysig +if skymod gt 0. and skysig gt 0. and skymod-2.*skysig gt 0. then begin + prelim_intspec = prelim_intspec - (skymod-2.*skysig) + kcwi_print_info,ppar,pre,'subtracting scattered light offset of', $ + (skymod-2.*skysig),format='(a,f9.3)' +endif +; +if ppar.display ge 2 then begin + plot,prelim_subwvl,prelim_spec/max(prelim_spec),charsi=si,charthi=th,thick=th, $ + xthick=th, xtitle='Wave(A)', $ + ythick=th, ytitle='Rel. Flux',title=imglab+' ('+kgeom.refname+'), No Offset',/xs + oplot,prelim_refwvl,prelim_refspec/max(prelim_refspec),color=colordex('red'),thick=th + oplot,[cwvl,cwvl],!y.crange,color=colordex('green'),thick=th,linesty=2 + kcwi_legend,['Ref Bar ('+strn(refbar)+')','Atlas','CWAVE'],linesty=[0,0,2], $ + thick=[th,th,th],box=0,charsi=si,charthi=th, $ + color=[colordex('black'),colordex('red'),colordex('green')] + read,'next: ',q +endif +; +; let's apply cosine bell taper to both +prelim_intspec = prelim_intspec * tukeywgt(n_elements(prelim_intspec),ppar.taperfrac) +prelim_refspec = prelim_refspec * tukeywgt(n_elements(prelim_refspec),ppar.taperfrac) +; +; now we have two spectra we can try to cross-correlate +; (prelim_intspec and prelim_refspec), so let's do that: +if ddisplay then window,1,title='kcwi_xspec' +kcwi_xspec,prelim_intspec,prelim_refspec,ppar,prelim_offset,prelim_value, $ + /min,/shift,/plot,label='Obj(0) vs '+kgeom.refname+'(1)',central=5. +if ddisplay then wset,0 +; +; record initial offset +kcwi_print_info,ppar,pre,'Initial arc-atlas offset (px, Ang)',prelim_offset, $ + prelim_offset*refdisp,format='(a,1x,f9.2,1x,f9.2)' +if ppar.display ge 2 then begin + q='test' + while strlen(q) gt 0 do begin + plot,prelim_subwvl-prelim_offset*refdisp,prelim_spec/max(prelim_spec), $ + charsi=si,charthi=th,thick=th,xthick=th,xtitle='Wave(A)', $ + ythick=th,ytitle='Rel. Flux',title=imglab+' ('+kgeom.refname+'), Offset = ' + $ + strtrim(string(prelim_offset*refdisp,form='(f9.2)'),2)+' Ang ('+$ + strtrim(string(prelim_offset,form='(f9.3)'),2)+' px)',/xs + oplot,prelim_refwvl,prelim_refspec/max(prelim_refspec), $ + color=colordex('red'),thick=th + oplot,[cwvl,cwvl],!y.crange,color=colordex('green'),thick=th,linesty=2 + kcwi_legend,['Ref Bar ('+strn(refbar)+')','Atlas','CWAVE'],linesty=[0,0,2], $ + thick=[th,th,th],box=0,charsi=si,charthi=th, $ + color=[colordex('black'),colordex('red'),colordex('green')] + read,'Enter: - next, new offset (px): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then $ ; just in case user enters 'q' + q = '' + if strlen(q) gt 0 then $ + prelim_offset = float(q) + endwhile +endif +; +; record final offset +kcwi_print_info,ppar,pre,'Final arc-atlas offset (px, Ang)',prelim_offset, $ + prelim_offset*refdisp,format='(a,1x,f9.2,1x,f9.2)' +; +; At this point we have the offsets between bars and the approximate offset from +; the reference bar to the actual spectrum and the approximate +; dispersion. +; +; let's populate the 0 points array. +p0 = cwvl + kgeom.baroff*prelim_disp - prelim_offset * refdisp +; +; next we are going to brute-force scan around the preliminary +; dispersion for a better solution. We will wander 5% away from it. +; +;we will try nn values +max_ddisp = 0.05d ; fraction (0.05 equiv to 5%) +;nn = (fix((1+max_ddisp)*max_ddisp*abs(prelim_disp)/refdisp*(maxrow-minrow)/2.0))>10<25 +nn = (fix(max_ddisp*abs(prelim_disp)/refdisp*(maxrow-minrow)/3.0))>10<25 +delta = (nn/10)>2;<3 ; may want to adjust this more? +kcwi_print_info,ppar,pre,'N disp. samples: ',nn +; +; which are: +disps = prelim_disp * ( 1.0d + max_ddisp * (dindgen(nn+1)-double(nn)/2.0d)*2.0d/double(nn)) +; +;containers for output values +maxima = dblarr(nn+1) +shifts = maxima +maxidx = dindgen(nn+1) +dspstat = intarr(nn+1) +; +; containers for bar-specific values +bardisp = dblarr(120) +barshift = dblarr(120) +barstat = intarr(120) +; +; data coefficients +coeff = dblarr(9) +; +; x values for central fit +subxvals = xvals[minrow:maxrow] +; +; loop over bars +for b = 0,119 do begin + ; + ; get sub spectrum for this bar + subspec = reform(specs[minrow:maxrow,b]) + ; + ; now loop over the dispersions... + for d = 0, nn do begin + ; + ; populate the coefficients + coeff[0] = p0[b] + coeff[1] = disps[d] + beta = acos( (coeff[1]/(pix*ybin)*rho*fcam*1d-4) < 1.d) + coeff[2] = -(pix*ybin/fcam)^2*sin(beta)/2.0d/rho*1d4 + coeff[3] = -(pix*ybin/fcam)^3*cos(beta)/6.0d/rho*1d4 + coeff[4] = (pix*ybin/fcam)^4*sin(beta)/24.0d/rho*1d4 + ; + ; what are the minimum and maximum wavelengths to consider + ; for the bar? + ; + minwvl = min( [poly(xvals[minrow],coeff), poly(xvals[maxrow],coeff)] ) + maxwvl = max( [poly(xvals[minrow],coeff), poly(xvals[maxrow],coeff)] ) + ; + ; where will we need to interpolate to cross correlate? + qwvl = where(refwvl gt minwvl and refwvl lt maxwvl, nqwvl) + if nqwvl le 0 then begin + kcwi_print_info,ppar,pre,'Insufficient reference wavelengths',/err + kgeom.status=3 + return + endif + ; + subrefwvl = refwvl[qwvl] + subrefspec = refspec[qwvl] + ; + ; get bell cosine taper to avoid nasty edge effects + tkwgt = tukeywgt(n_elements(subrefspec), ppar.taperfrac) + ; + ; apply taper to atlas spectrum + subrefspec = subrefspec * tkwgt + ; + ; adjust the spectra + waves = poly(subxvals,coeff) + ; + ; interpolate the bar spectrum + intspec = interpol(subspec,waves,subrefwvl,/spline) * prelim_disp/disps[d] + ; + ; apply taper to bar spectrum + intspec = intspec * tkwgt + ; + ; get a label + xslab = 'Bar '+strn(b)+', '+strn(d)+'/'+strn(nn)+', Dsp = '+string(disps[d],form='(f6.3)') + ; + ; cross correlate the interpolated spectrum with the atlas spectrum + if ddisplay then wset,1 + kcwi_xspec,intspec,subrefspec,ppar,soffset,svalue,status=status, $ + /min,/shift,/pad,/central,plot=ddisplay,label=xslab + if ppar.display ge 3 then wset,0 + ; + maxima[d] = double(svalue)/total(subrefspec)/total(intspec) + shifts[d] = soffset + dspstat[d] = status + ; + endfor ; d + ; + ; now find the max of the + ; cross-correlation scan and determine + ; what the corresponding shift and + ; dispersion are + mx = max(maxima,mi) + ; + submax = maxima[mi-delta>0:mi+delta0:mi+delta - next): ',q $ + else read,'Next? (Q - quit plotting, - next): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then begin + display = (1 eq 0) + if ppar.display ge 3 then begin + ddisplay = (1 eq 0) + wdelete,1 + endif + endif + if strupcase(strmid(q,0,1)) eq 'D' then ddisplay = (1 eq 1) + endif +endfor ; b +; +; clean up +if ddisplay then wdelete,1 +; +; now clean each slice of outlying bars +if ppar.cleancoeffs then $ + kcwi_clean_coeffs,centcoeff,degree,ppar +; +; let's check the status of the fits +if total(barstat) gt n_elements(barstat)/2 then begin + kcwi_print_info,ppar,pre,'Too many bar fit failures', $ + fix(total(barstat)),/error + kgeom.status = 4 +endif else begin + kgeom.status = 0 + ; + ; plot results + if ppar.display ge 1 then begin + ; + ; image label + imglab = 'Img # '+strn(imgnum)+' ('+kgeom.refname+') Sl: '+ $ + strtrim(ifunam,2)+ ' Fl: '+strtrim(filter,2)+' Gr: '+ $ + strtrim(grating,2) + if !d.window lt 0 then $ + window,0,title=kcwi_drp_version() $ + else wset,0 + !p.multi=[0,1,2] + si = 1.5 + ys = reform(centcoeff[0,*]) + yrng = get_plotlims(ys) + plot,centcoeff[0,*],psym=4,charsi=si,charthi=th,thick=th, $ + title = imglab+' Central Fit Coef0', $ + xthick=th,xtitle='Bar #', xrange=[-1,120],/xs, $ + ythick=th,ytitle='Ang',yrange=yrng,/ys + kcwi_oplot_slices + ys = reform(centcoeff[1,*]) + yrng = get_plotlims(ys) + plot,centcoeff[1,*],psym=4,charsi=si,charthi=th,thick=th, $ + title = imglab+' Central Fit Coef1', $ + xthick=th,xtitle='Bar #', xrange=[-1,120],/xs, $ + ythick=th,ytitle='Ang/px',yrange=yrng,/ys + kcwi_oplot_slices + if ppar.display ge 2 then $ + read,'next: ',q + !p.multi=0 + endif +endelse +; +return +end ; kcwi_fit_center diff --git a/pro/kcwi_flatten_cube.pro b/pro/kcwi_flatten_cube.pro new file mode 100644 index 0000000..4b83d9e --- /dev/null +++ b/pro/kcwi_flatten_cube.pro @@ -0,0 +1,134 @@ +; +; Copyright (c) 2017, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_FLATTEN_CUBE +; +; PURPOSE: +; Take a 3-D cube and output a 2-D version +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_FLATTEN_CUBE, Cfile +; +; Returns: +; 2-d or 3-d image, or -1 if image file not found +; +; INPUTS: +; cfile - Image cube file name +; +; OUTPUTS: +; +; KEYWORDS: +; ISCALE - set to scale output into integers +; +; SIDE EFFECTS: +; Writes a 2-D image with name based on input. +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2017-APR-27 Initial version +;- +pro kcwi_flatten_cube,cfile,iscale=iscale + ; + ; setup + pre = 'KCWI_FLATTEN_CUBE' + ; + ; check inputs + if n_params(0) lt 1 then begin + cfile = '' + read,'Input cube image file name: ',cfile + endif + ; + ; check if it exists + if file_test(cfile) then begin + ; + ; read in cube + cub = mrdfits(cfile,0,hdr,/fscale,/silent) + ; + ; is it 3-D? + sz = size(cub,/dimen) + if n_elements(sz) eq 3 then begin + ; + ; output image + oim = fltarr(sz[0]*sz[1],sz[2]) + ; + ; pack slices + for i=0,sz[0]-1 do $ + oim[i*sz[1]:(i+1)*sz[1]-1,*] = cub[i,*,*] + ; + ; get output name + ofil = repstr(cfile,'.fits','_2d.fits') + ; + ; get wavelength values + w0 = sxpar(hdr,'CRVAL3') + dw = sxpar(hdr,'CD3_3') + crpixw = sxpar(hdr,'CRPIX3') + ; + ; set spatial scale + s0 = 0. + ds = 24.0 / (sz[0]*sz[1]) + crpixs = 1.0 + ; + ; update header + sxaddpar,hdr,'HISTORY',' '+pre+' '+systime(0) + sxaddpar,hdr,'INCUBEF',cfile,' Inpute cube filename' + ; + ; remove old WCS + sxdelpar,hdr,'NAXIS3' + sxdelpar,hdr,'CTYPE1' + sxdelpar,hdr,'CTYPE2' + sxdelpar,hdr,'CTYPE3' + sxdelpar,hdr,'CUNIT1' + sxdelpar,hdr,'CUNIT2' + sxdelpar,hdr,'CUNIT3' + sxdelpar,hdr,'CNAME1' + sxdelpar,hdr,'CNAME2' + sxdelpar,hdr,'CNAME3' + sxdelpar,hdr,'CRVAL1' + sxdelpar,hdr,'CRVAL2' + sxdelpar,hdr,'CRVAL3' + sxdelpar,hdr,'CRPIX1' + sxdelpar,hdr,'CRPIX2' + sxdelpar,hdr,'CRPIX3' + sxdelpar,hdr,'CD1_1' + sxdelpar,hdr,'CD1_2' + sxdelpar,hdr,'CD2_1' + sxdelpar,hdr,'CD2_2' + sxdelpar,hdr,'CD3_3' + ; + ; set wavelength axis WCS values + sxaddpar,hdr,'WCSDIM',2 + sxaddpar,hdr,'CTYPE1','SPATIAL',' SLICE' + sxaddpar,hdr,'CUNIT1','slu',' SLICE units' + sxaddpar,hdr,'CNAME1','KCWI SLICE',' SLICE name' + sxaddpar,hdr,'CRVAL1',s0,' SLICE zeropoint' + sxaddpar,hdr,'CRPIX1',crpixs,' SLICE reference pixel' + sxaddpar,hdr,'CDELT1',ds,' SLICE per pixel' + sxaddpar,hdr,'CTYPE2','AWAV',' Air Wavelengths' + sxaddpar,hdr,'CUNIT2','Angstrom',' Wavelength units' + sxaddpar,hdr,'CNAME2','KCWI 2D Wavelength',' Wavelength name' + sxaddpar,hdr,'CRVAL2',w0,' Wavelength zeropoint' + sxaddpar,hdr,'CRPIX2',crpixw,' Wavelength reference pixel' + sxaddpar,hdr,'CDELT2',dw,' Wavelength Angstroms per pixel' + ; + ; write it out + mwrfits,oim,ofil,hdr,/create,iscale=iscale + endif else begin + print,pre+': not a 3-D image cube - ', cfile + endelse + ; + ; report non-existence + endif else begin + print,pre+': file not found - ',cfile + endelse + ; + return +end diff --git a/pro/kcwi_geom__define.pro b/pro/kcwi_geom__define.pro new file mode 100644 index 0000000..2244a70 --- /dev/null +++ b/pro/kcwi_geom__define.pro @@ -0,0 +1,198 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_GEOM__DEFINE +; +; PURPOSE: +; This procedure defines the structure for all the information for +; KCWI geometric transformations. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; kgeom = {KCWI_GEOM} +; +; INPUTS: +; None +; +; OUTPUTS: +; None +; +; SIDE EFFECTS: +; defines the KCWI_GEOM data structure +; +; PROCEDURE: +; Provide the automatic structure definition for the KCWI geometric +; transformation structure KCWI_GEOM. +; +; EXAMPLE: +; Instantiate and initialize a single geom struct for KCWI: +; +; kgeom = {kcwi_geom} +; kgeom = struct_init(kgeom) +; +; NOTES: +; Keep structure name on a line by itself (see struct_init.pro). +; Keep tag names 15 chars or less. +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-MAY-03 Initial version +; 2013-AUG-12 Added padlo, padhi tags +; 2016-APR-22 Added rotoff tag +;- +pro kcwi_geom__define +; +tmp = { kcwi_geom, $ +; +; output filename + geomfile:'', $ ; output save file for geom struct +; +; calibration images + cbarsimgnum:0l, $ ; cbars image number (0 - 9999) + cbarsjd:0.d0, $ ; cbars image julian date + cbarsfname:'', $ ; cbars image file name + arcimgnum:0l, $ ; arc image number + arcjd:0.d0, $ ; arc image julian date + arcfname:'', $ ; arc image file name +; +; reference spectra + refname:'', $ ; name of reference spectrum ('ThAr', e.g.) + refspec:'', $ ; reference spectrum file + reflist:'', $ ; reference line list +; +; nod and shuffle? + nasmask:0, $ ; 0 - no, 1 - yes +; +; IFU properties + ifunum: 0, $ ; Slicer number (1-3, Large, Medium, Small) + ifunam: '', $ ; Slicer name ("Large", etc.) +; +; grating properties + gratid:'', $ ; grating id + gratnum:-1, $ ; grating number + rho:3.0d, $ ; grating lines/mm + adjang:0.0d, $ ; grating adjustment angle (0 or 180 deg) + lastdegree:4, $ ; highest order for full-ccd wavelength fit +; +; filter properties + filter:'', $ ; filter id + filtnum:-1, $ ; filter number +; +; encoder positions + campos:0L, $ ; Camera encoder position + camang:0.0d, $ ; Camera articulation angle in degrees + grenc:0L, $ ; grating rotator encoder position + grangle:0., $ ; grating rotator angle (degrees) + gratanom:0., $ ; grating angle anomoly (degrees) +; +; pixel scale + pxscl:0.00008096d0, $ ; degrees per unbinned spatial pixel +; +; slice scale + slscl:0.00075437d0, $ ; degrees per slice +; +; IFU rotator offset + rotoff:0.0, $ ; degrees from rotator PA +; +; CCD binning + xbinsize:1, $ ; binning in x (pixels) + ybinsize:1, $ ; binning in y (pixels) +; +; CCD size + nx:0, $ ; number of x pixels in image + ny:0, $ ; number of y pixels in image +; +; CCD readnoise + rdnoise:3., $ ; e-/pixel +; +; Y ranges + goody0:0, $ ; lowest good y pixel in spectrum + goody1:0, $ ; highest good y pixel in spectrum + trimy0:0, $ ; lowest y pixel to include in output + trimy1:0, $ ; highest y pixel to include in output + goodmy0:0, $ ; lowest good y pixel in masked spectrum + goodmy1:0, $ ; highest good y pixel in masked spectrum + trimmy0:0, $ ; lowest y pixel to include in masked output + trimmy1:0, $ ; highest y pixel to include in masked output +; +; cross correlation parameters + ccwn:100., $ ; cross correlation window for bar-to-bar offset + ccoff:fltarr(24), $ ; known offsets relative to reference slice in pixels +; +; wavelengths + cwave:-9., $ ; central wavelength (Angstroms) + wave0out:-9., $ ; output wavelength zeropoint + wave1out:-9., $ ; output ending wavelength + dwout:-9., $ ; output dispersion (Ang/pix) + wavran:-9., $ ; approximate wavelength range (Ang) + halfwidth:3, $ ; instrumental half-width of resolution (Ang) + resolution:1., $ ; Gaussian sigma of instrumental resol. (Ang) + atsig:1., $ ; Gaussian sigma for convolving with atlas (Ang) + pxwindow:0, $ ; window size for finding lines in pixels + waveall0:-9., $ ; low wavelength that includes all data + waveall1:-9., $ ; high wavelength that includes all data + wavegood0:-9., $ ; low wavelength that includes all good data + wavegood1:-9., $ ; high wavelength that includes all good data + waveallm0:-9., $ ; low wavelength that includes all masked data + waveallm1:-9., $ ; high wavelength that includes all masked data + wavegoodm0:-9., $ ; low wavelength that includes all good masked data + wavegoodm1:-9., $ ; high wavelength that includes all good masked data + wavemid:-9., $ ; wavelength in the middle of the ranges + wavefid:3000., $ ; fiducial wavelength for aligning all solns. +; +; control points + xi:fltarr(6120), $ ; x input + yi:fltarr(6120), $ ; y input + xo:fltarr(6120), $ ; x output + yo:fltarr(6120), $ ; y output + xw:fltarr(6120), $ ; x output in wavelength space + yw:fltarr(6120), $ ; y output in wavelength space + bar:intarr(6120), $ ; bar number for each point + slice:intarr(6120), $ ; slice number for each point +; +; Reference bar for wavelength solution + refbar:57, $ ; 0 - 119 (defaults to 57) + rbcoeffs:fltarr(9), $ ; polynomial wavelength solution up to 9 coeffs + refoutx:fltarr(5), $ ; ref output x control point positions (pixels) + refdelx:-9., $ ; ref delta x (pixels) +; +; X values for bars + barx:fltarr(120), $ ; measured in cbars image + refx:fltarr(120), $ ; reference x for extracted arc spectra + x0out:30, $ ; output spatial zeropoint (unbinned pix) +; +; Y offset for each bar + baroff:fltarr(120), $ ; pixel offset relative to refbar +; +; Wavelength fit for each bar + bfitcoeffs:fltarr(9,120),$ ; polynomial wavelength solution + bfitord:6, $ ; order of polynomial wavelength solution + bclean:0, $ ; clean coeffs of errant bars? 0 - no, 1 - yes +; +; Which image row was used for canonical x position? + midrow:0L, $ ; image pixel +; +; Coefficients + kx:dblarr(4,4), $ ; simple de-warping + ky:dblarr(4,4), $ ; + kwx:dblarr(6,6,24), $ ; wavelength solution included for each slice + kwy:dblarr(6,6,24), $ ; + ypad:600, $ ; number of pixels to pad image in y +; +; Status + xrsd:fltarr(24), $ ; transform x residuals (pixels) + yrsd:fltarr(24), $ ; transform y residuals (pixels) + avewavesig:0., $ ; average bar wavelength sigma in Angstroms + stdevwavesig:0., $ ; standard deviation of bar sigmas in Angstroms + initialized:0, $ ; 0 - no, 1 - yes + status:-1, $ ; 0 - good fit, else not good + progid:'', $ ; program that last modified the geom file +; +; timestamp + timestamp:0.d0 $ ; timestamp for struct (unix seconds) + } +end diff --git a/pro/kcwi_geom_resid.pro b/pro/kcwi_geom_resid.pro new file mode 100644 index 0000000..de71a25 --- /dev/null +++ b/pro/kcwi_geom_resid.pro @@ -0,0 +1,74 @@ +; +; Copyright (c) 2014, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_GEOM_RESID +; +; PURPOSE: +; Calculates the residuals of the geometric transformation +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_GEOM_RESID,Xi,Yi,Xw,Yw,Degree,Kwx,Kwy,Xrsd,Yrsd +; +; INPUTS: +; Xi,Yi - Initial control points +; Xw,Yw - Wavelength corrected control points +; Degree - Degree passed to POLYWARP +; Kwx,Kwy - Output coeffs from POLYWARP +; +; OUTPUTS: +; Xrsd,Yrsd - Differences between Xi,Yi inputs and calculated +; points from coeffs +; +; SIDE EFFECTS: +; None +; +; EXAMPLE: +; polywarp,xi,yi,xw,yw,3,kwx,kwy,/double,status=status +; kcwi_geom_resid,xi,yi,xw,yw,3,kwx,kwy,xrsd,yrsd +; !p.multi=[0,1,2] +; plot,xi,xrsd,psym=4 +; plot,yi,yrsd,psym=4 +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2014-OCT-01 Initial Revision +;- +; +pro kcwi_geom_resid,xi,yi,xw,yw,degree,kwx,kwy,xrsd,yrsd +; +; setup +pre = 'KCWI_GEOM_RESID' +xrsd = -1. +yrsd = -1. +; +np = n_elements(xi) +if n_elements(yi) ne np or n_elements(xw) ne np or $ + n_elements(yw) ne np then begin + print,pre,': ERROR - control point vectors must all have the same number of elements' + return +endif +; +kwxsz = size(kwx,/dim) +kwysz = size(kwy,/dim) +if kwxsz[0] ne degree[1]+1 or kwxsz[1] ne degree[0]+1 or $ + kwysz[0] ne degree[1]+1 or kwysz[1] ne degree[0]+1 then begin + print,pre,': ERROR - 2D coeff array dimension error' + return +endif +; +xp=0. +yp=0. +for i=0,degree[1] do for j=0,degree[0] do begin + xp = xp + kwx[i,j]*xw^j*yw^i + yp = yp + kwy[i,j]*xw^j*yw^i +endfor +xrsd = xi - xp +yrsd = yi - yp +; +return +end diff --git a/pro/kcwi_get_atlas.pro b/pro/kcwi_get_atlas.pro new file mode 100644 index 0000000..3fa8041 --- /dev/null +++ b/pro/kcwi_get_atlas.pro @@ -0,0 +1,36 @@ +;+ +; kcwi_get_atlas - derive which atlas to use based on kcfg +; +; INPUTS: +; kcfg - KCWI_CFG struct +; +; OUTPUTS: +; atlas - atlas spectrum fits file +; atname - atlas name +;- +pro kcwi_get_atlas,kcfg,atlas,atname + ; + ; setup + pre = 'KCWI_GET_ATLAS' + ; + ; verify Kcfg + if kcwi_verify_cfg(kcfg,/silent) ne 0 then begin + print,pre+': Error - malformed KCWI_CFG struct' + return + endif + ; + ; Init + atlas = 0 + atname = 0 + ; + ; check Lamp 0 + if kcfg.lmp0stat eq 1 and kcfg.lmp0shst eq 1 then begin + atname = strtrim(kcfg.lmp0nam,2) + atlas = strlowcase(atname)+'.fits' + endif else if kcfg.lmp1stat eq 1 and kcfg.lmp1shst eq 1 then begin + atname = strtrim(kcfg.lmp1nam,2) + atlas = strlowcase(atname)+'.fits' + endif + ; + return +end diff --git a/pro/kcwi_get_imname.pro b/pro/kcwi_get_imname.pro new file mode 100644 index 0000000..0ab3370 --- /dev/null +++ b/pro/kcwi_get_imname.pro @@ -0,0 +1,75 @@ +function kcwi_get_imname,ppar,imgnum,tail, $ + nodir = nodir, calib = calib, raw = raw, reduced = reduced, $ + master = master, exist=exist, help = help + ; + ; setup + pre = 'KCWI_GET_IMNAME' + ; + ; help request + if n_params(0) lt 2 or keyword_set(help) then begin + print,pre+': Info - Usage: full_filename = '+pre+'( KcwiPpar, ImgNum, [Tail], /CALIB, /INPUT, /OUTPUT, /MASTER)' + return,'' + endif + ; + ; verify ppar + if kcwi_verify_ppar(ppar,/init) ne 0 then begin + print,pre+': ERROR - invalid or uninitialized ppar struct' + return,'' + endif + ; + ; test imgnum + if imgnum lt 0 then begin + kcwi_print_info,ppar,pre,'invalid image number',/error + return,'' + endif + ; + ; check tail + if n_params(0) lt 3 then $ + tail = '' + ; + ; default is raw dir + dir = ppar.rawdir + if keyword_set(reduced) then $ + dir = ppar.reddir + if keyword_set(calib) then $ + dir = ppar.caldir + if keyword_set(nodir) then $ + dir = '' + ; + ; image number and root image name strings + ; + ; are we a master file? (mbias, mdark, mflat) + if keyword_set(master) then begin + if strlen(master) gt 0 then begin + rutstr = master + endif else begin + kcwi_print_info,ppar,pre,'invalid master file root',/error + return,'' + endelse + ; + ; regular image file + endif else begin + rutstr = ppar.froot + endelse + imgstr = string(imgnum,'(i0'+strn(ppar.fdigits)+')') + ; + ; construct file name + file = dir + rutstr + imgstr + tail + '.fits' + ; + ; does it exist? + if keyword_set(exist) then begin + if not file_test(file,/regular,/read) then begin + ofil = file + file = file + '.gz' + if not file_test(file,/regular,/read) then begin + kcwi_print_info,ppar,pre, $ + 'File not found: '+ofil+' nor '+file,/error + file = ofil + endif + endif + endif + ; + ; all done + return,file + +end diff --git a/pro/kcwi_group_dgeom.pro b/pro/kcwi_group_dgeom.pro new file mode 100644 index 0000000..5ae9440 --- /dev/null +++ b/pro/kcwi_group_dgeom.pro @@ -0,0 +1,123 @@ +; +; Copyright (c) 2014, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_GROUP_DGEOM +; +; PURPOSE: +; This procedure groups continuum bars (cbars) and arcs in the KCWI_CFG +; struct for a given night. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_GROUP_DGEOM, Kcfg, Ppar, Ccfg, Acfg +; +; INPUTS: +; Kcfg - array of struct KCWI_CFG for a given directory +; Ppar - KCWI_PPAR pipeline parameter struct +; +; OUTPUTS: +; Ccfg - a KCWI_CFG struct vector with one entry for each arcbars from +; a calibration set +; Acfg - a KCWI_CFG struct vector with one entry for each arc from +; a calibration set +; +; KEYWORDS: +; +; SIDE EFFECTS: +; Updates Ppar struct with number of geometry groups, number of +; arcs and number continuum bars (ndirect, ndarcs). +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2014-APR-01 Initial version +;- +pro kcwi_group_dgeom, kcfg, ppar, ccfg, acfg + ; + ; setup + pre = 'KCWI_GROUP_DGEOM' + ppar.ndggrps = 0 + ; + ; check input + if kcwi_verify_cfg(kcfg) ne 0 then return + if kcwi_verify_ppar(ppar) ne 0 then return + ; + ; required calibration image types + bg = where(kcfg.imgtype eq 'arcbars' and $ + strpos(kcfg.obstype,'direct') ge 0, nbar) + ag = where(kcfg.imgtype eq 'arc' and $ + strpos(kcfg.obstype,'direct') ge 0, narc) + ; + ; check if we have required types + if narc le 0 or nbar le 0 then begin + ; + ; record results + if narc le 0 then $ + kcwi_print_info,ppar,pre,'no direct arcs found' + if nbar le 0 then $ + kcwi_print_info,ppar,pre,'no direct arcbars found' + return + endif + ; + ; get config records + cfg = kcfg[bg] + afg = kcfg[ag] + ; + ; collect matches + bmatch = lonarr(nbar) - 1L + amatch = lonarr(nbar) - 1L + ; + ; loop over bars images and gather direct geom pairs + for i=0,nbar-1 do begin + ; + ; loop over arcs and find a good match + for j=0,narc-1 do begin + ; + ; check configuration + tcfg = kcwi_match_cfg(afg[j],cfg[i],ppar,count=nm,/silent) + ; + ; check for match + if nm eq 1 then begin + ; + ; are we close in image number? + if abs(cfg[i].imgnum - afg[j].imgnum) lt 2 then begin + bmatch[i] = i + amatch[i] = j + endif ; end check if we are close in image number + endif ; end check for match + endfor ; end loop over arcs + endfor ; end loop over cbars + ; + ; get the good ones + good = where(bmatch ge 0, ndirect) + ; + ; collect the good calibs + if ndirect gt 0 then begin + ; + ; cbars + ccfg = cfg[bmatch[good]] + ; + ; arcs + acfg = afg[amatch[good]] + for i=0,ndirect-1 do $ + kcwi_print_info,ppar,pre,'Direct config', $ + (kcwi_cfg_string(ccfg[i],/long,/delim))[0] + endif else begin + acfg = -1 + ccfg = -1 + kcwi_print_info,ppar,pre,'no direct geom image sets found',/warning + endelse + ; + ; report number of direct geom groups + ppar.ndggrps = ndirect + kcwi_print_info,ppar,pre,'Number of direct geom groups',ndirect + ; + return +end diff --git a/pro/kcwi_group_geom.pro b/pro/kcwi_group_geom.pro new file mode 100644 index 0000000..209f178 --- /dev/null +++ b/pro/kcwi_group_geom.pro @@ -0,0 +1,143 @@ +; +; Copyright (c) 2014, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_GROUP_GEOM +; +; PURPOSE: +; This procedure groups continuum bars (cbars) and arcs in the KCWI_CFG +; struct for a given night. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_GROUP_GEOM, Kcfg, Ppar, Ccfg, Acfg +; +; INPUTS: +; Kcfg - array of struct KCWI_CFG for a given directory +; Ppar - KCWI_PPAR pipeline parameter struct +; +; OUTPUTS: +; Ccfg - a KCWI_CFG struct vector with one entry for each cbars from +; a calibration set +; Acfg - a KCWI_CFG struct vector with one entry for each arc from +; a calibration set +; +; KEYWORDS: +; +; SIDE EFFECTS: +; Updates Ppar struct with number of geometry groups, number of +; arcs and number continuum bars (ngeom, narcs, nbars). +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2014-APR-01 Initial version +;- +pro kcwi_group_geom, kcfg, ppar, ccfg, acfg + ; + ; setup + pre = 'KCWI_GROUP_GEOM' + ppar.nggrps = 0 + ; + ; check input + if kcwi_verify_cfg(kcfg) ne 0 then return + if kcwi_verify_ppar(ppar) ne 0 then return + ; + ; required calibration image types + bg = where(kcfg.imgtype eq 'cbars' and $ + strpos(kcfg.obstype,'direct') lt 0, nbar) + ag = where(kcfg.imgtype eq 'arc' and $ + strpos(kcfg.obstype,'direct') lt 0, narc) + ; + ; check if we have required types + if narc le 0 or nbar le 0 then begin + ; + ; record results + if narc le 0 then $ + kcwi_print_info,ppar,pre,'no arcs found!',/error + if nbar le 0 then $ + kcwi_print_info,ppar,pre,'no bars found!',/error + return + endif + ; + ; get config records + cfg = kcfg[bg] + afg = kcfg[ag] + ; + ; collect matches + bmatch = lonarr(nbar) - 1L + amatch = lonarr(nbar) - 1L + lamp = strarr(nbar) + ; + ; loop over bars images and gather geom pairs + for i=0,nbar-1 do begin + ; + ; loop over arcs and find a good match + for j=0,narc-1 do begin + ; + ; check configuration + tcfg = kcwi_match_cfg(afg[j],cfg[i],ppar,count=nm,/silent) + ; + ; check for match + if nm eq 1 then begin + ; + ; are we close in image number? + if abs(cfg[i].imgnum - afg[j].imgnum) lt 4 then begin + bmatch[i] = i + if strpos(cfg[i].gratid,'BH') ge 0 or $ + strpos(cfg[i].gratid,'BL') ge 0 or $ + strpos(cfg[i].gratid,'BM') ge 0 then begin + if afg[j].lmp1stat eq 1 and afg[j].lmp1shst eq 1 then begin + amatch[i] = j + lamp[i] = afg[j].lmp1nam + endif else $ + if amatch[i] lt 0 then amatch[i] = j + endif else begin + if afg[j].lmp0stat eq 1 and afg[j].lmp0shst eq 1 then begin + amatch[i] = j + lamp[i] = afg[j].lmp0nam + endif else $ + if amatch[i] lt 0 then amatch[i] = j + endelse + endif ; end check if we are close in image number + endif ; end check for match + endfor ; end loop over arcs + endfor ; end loop over cbars + ; + ; get the good ones + good = where(bmatch ge 0, ngeom) + ; + ; collect the good calibs + if ngeom gt 0 then begin + kcwi_print_info,ppar,pre,'Number of geom groups',ngeom + ; + ; cbars + ccfg = cfg[bmatch[good]] + ; + ; arcs + acfg = afg[amatch[good]] + ; + ; lamps + lamp = lamp[good] + for i=0,ngeom-1 do $ + kcwi_print_info,ppar,pre,'BarImg, Geom config, lamp', $ + ccfg[i].imgnum, $ + ', '+(kcwi_cfg_string(ccfg[i],/delim,/long))[0] +', '+lamp[i], $ + format='(a,i7,a)' + endif else begin + acfg = -1 + ccfg = -1 + kcwi_print_info,ppar,pre,'no geom groups found',/warning + endelse + ; + ; report number of geom groups + ppar.nggrps = ngeom + ; + return +end diff --git a/pro/kcwi_plot_arcfits.pro b/pro/kcwi_plot_arcfits.pro new file mode 100644 index 0000000..37086d6 --- /dev/null +++ b/pro/kcwi_plot_arcfits.pro @@ -0,0 +1,356 @@ +; +; Copyright (c) 2014, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_PLOT_ARCFITS +; +; PURPOSE: +; plot arc wavelength solution diagnostics +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_PLOT_ARCFITS, Specs, Kgeom, Ppar, CntCoeff, FinCoeff, Sigmas, Fwaves, Dwaves +; +; INPUTS: +; Specs - a array of arc spectra produced by KCWI_EXTRACT_ARCS +; Kgeom - KCWI_GEOM struct from KCWI_TRACE_CBARS and KCWI_EXTRACT_ARCS +; Ppar - KCWI_PPAR pipeline parameter struct +; CntCoeff - Coefficients from fits to central region +; FinCoeff - Coefficients from fits to full region +; Sigmas - Line center residuals (Ang) of fits +; Fwaves - Final fit wavelengths of matched lines +; Dwaves - Final difference between fit and atlas wavelengths +; +; INPUT KEYWORDS: +; TWEAK - set to indicate the fits are from iteratively tweaked fits +; PLOT_FILE - if set to a string, will be used to produce postscript +; output of diagnostic plots +; Ftype - Label describing fit type: FullCCD, Central, etc. +; +; SIDE EFFECTS: +; outputs postscript plot file specified in PLOT_FILE +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Matt Matuszewski, Don Neill +; 2014-SEP-19 Initial Revision +;- +; +pro kcwi_plot_arcfits, specs, kgeom, ppar, cntcoeff, fincoeff, sigmas, $ + fwaves, dwaves, tweak=tweak, plot_file=plot_file, ftype=ftype + +pre = 'KCWI_PLOT_ARCFITS' +q='' +; +; do we want to display stuff? +display = (ppar.display ge 2 or keyword_set(plot_file)) +; +; what to do if we are not plotting or displaying +if not display then return +; +; which image number +imgnum = kgeom.arcimgnum +; +; which slicer +ifunam = kgeom.ifunam +imglab = 'Img # '+strn(imgnum)+' ('+kgeom.refname+') Sl: '+strtrim(ifunam,2)+ $ + ' Fl: '+strtrim(kgeom.filter,2)+' Gr: '+strtrim(kgeom.gratid,2) +; +; is the N+S mask in? +nasmask = kgeom.nasmask +; +; degree defaults to 3 for central region fit +degree = 3 +; +; check if we need the last degree for full-ccd tweaked fits +if keyword_set(tweak) and not nasmask then $ + degree = kgeom.lastdegree +; +; fit type +if keyword_set(ftype) then begin + fittype = ftype +endif else begin + if nasmask then $ + fittype = 'Central' $ + else fittype = 'FullCCD' +endelse +; +; load the reference atlas spectrum. +kcwi_read_atlas,kgeom,ppar,refspec,refwvl,refdisp +; +; input spectrum dimensions +specsz = size(specs,/dim) +; +; set up array with zero point in the center +x0 = specsz[0]/2 +fxvals = dindgen(specsz[0]) ; fullCCD xvals +cxvals = fxvals-x0 ; central xvals +; +; if not making hardcopy, open a new window +if not keyword_set(plot_file) then $ + window,1,title='kcwi_plot_arcfits' +; +; set up plots +!p.multi=0 +if keyword_set(plot_file) then begin + psfile,plot_file + !p.multi=[0,1,2] + kcwi_print_info,ppar,pre,'plotting diagnostics to: '+plot_file+'.ps' +endif +deepcolor +!p.background=colordex('white') +!p.color=colordex('black') +th=2.0 +si=2.0 +; +; get good bars +bargood = where(sigmas ge 0., nbargood) +; +; bad bars +barbad = where(sigmas lt 0., nbarbad) +; +; number of comparison lines +nlines = intarr(120) +; +; first plot sigmas +if keyword_set(tweak) then begin + ; + ; get fit rms status + mom = moment(sigmas[bargood]) + fitrmslab = strtrim(string(mom[0],format='(f9.3)'),2) + ' +- ' + $ + strtrim(string(sqrt(mom[1]),format='(f9.3)'),2) + ; + ; set title + tlab=imglab+', '+fittype+' Fit : ' + fitrmslab + ; + ; plot range + yrng = get_plotlims(sigmas[bargood]) + ; + plot,sigmas,psym=4,title=tlab, $ + xtitle='Bar #',xrange=[-1,120],/xs, $ + ytitle='RMS (Ang)',yrange=yrng,/ys + oplot,!x.crange,[mom[0],mom[0]],linesty=5,thick=th + oplot,!x.crange,[mom[0]+sqrt(mom[1]),mom[0]+sqrt(mom[1])],linesty=1,thick=th + oplot,!x.crange,[mom[0]-sqrt(mom[1]),mom[0]-sqrt(mom[1])],linesty=1,thick=th + if nbarbad gt 0 then $ + oplot,barbad,replicate(mom[0],nbarbad),psym=7,thick=th + kcwi_oplot_slices + ; + ; next + if not keyword_set(plot_file) then $ + read,'next: ',q + ; + ; plot number of lines + for b = 0,119 do begin + ffwaves = reform(fwaves[b,*]) + li = where(ffwaves gt 0, nli) + nlines[b] = nli + endfor + mom = moment(nlines) + fitnlslab = strtrim(string(mom[0],format='(f9.1)'),2) + ' +- ' + $ + strtrim(string(sqrt(mom[1]),format='(f9.1)'),2) + ; + ; set title + tlab=imglab+', '+fittype+' Fit : ' + fitnlslab + ; + ; plot range + yrng = get_plotlims(nlines) + ; + plot,nlines,psym=4,title=tlab, $ + xtitle='Bar #',xrange=[-1,120],/xs, $ + ytitle='N Lines',yrange=yrng,/ys + oplot,!x.crange,[mom[0],mom[0]],linesty=5,thick=th + oplot,!x.crange,[mom[0]+sqrt(mom[1]),mom[0]+sqrt(mom[1])],linesty=1,thick=th + oplot,!x.crange,[mom[0]-sqrt(mom[1]),mom[0]-sqrt(mom[1])],linesty=1,thick=th + if nbarbad gt 0 then $ + oplot,barbad,replicate(mom[0],nbarbad),psym=7,thick=th + kcwi_oplot_slices + ; + ; next + if not keyword_set(plot_file) then $ + read,'next: ',q +endif +; +; plot each coeff +for i=0,degree do begin + if display then begin + ; + ; set title + tlab = imglab + ', '+fittype+' '+strn(i) + ; + ; set y title + if i eq 0 then begin + ylab = 'Ang' + endif else if i eq 1 then begin + ylab = 'Ang/pix' + endif else begin + ylab = 'Ang/pix^'+strn(i) + endelse + ; + ; plot range + yrng = get_plotlims(fincoeff[i,*]) + ; + plot,fincoeff[i,*],psym=4,title=tlab, $ + xtitle='Bar #',xrange=[-1,120],/xs, $ + ytitle=ylab,yrange=yrng,/ys + if nbarbad gt 0 then $ + oplot,barbad,fincoeff[i,barbad],psym=7,thick=th + kcwi_oplot_slices + if not keyword_set(plot_file) then begin + read,'Next? (Q-quit plotting, - next): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then $ + display = (1 eq 0) + endif + endif ; display? +endfor +; +if not keyword_set(plot_file) then $ + !p.multi = 0 $ +else !p.multi = [0,1,2] +; +; inclusive, all, and trim wavelengths +wavall0 = kgeom.waveall0 +wavall1 = kgeom.waveall1 +wavgood0 = kgeom.wavegood0 +wavgood1 = kgeom.wavegood1 +wavmid = kgeom.wavemid +cwave = kgeom.cwave +minwav = kgeom.wave0out +maxwav = kgeom.wave1out +; +; atlas name +atnam = kgeom.refname +; +; loop over each bar +for b=0,119 do begin + ; + ; are we displaying or plotting? + if display or keyword_set(plot_file) then begin + ; + ; fullCCD or Central? + if not nasmask then begin + cf = fincoeff[*,b] + pcolor = colordex('blue') + rcolor = colordex('s') + xvals = fxvals + endif else begin + cf = cntcoeff[*,b] + pcolor = colordex('red') + rcolor = colordex('pink') + xvals = cxvals + endelse + ; + ; check coeffs + gcf = where(finite(cf) eq 1, ngcf) + ; + ; generate wavelengths + if ngcf gt 0 then $ + waves = poly(xvals,cf[gcf]) $ + else waves = replicate(0.,n_elements(xvals)) + ; + ; mark off good regions for plot limits + ta = where(refwvl ge wavgood0 and refwvl le wavgood1) + ts = where(waves ge wavgood0 and waves le wavgood1,nts) + ; + ; are there good object wavelengths? + if nts gt 0 then begin + yrng = get_plotlims(specs[ts,b],/minzero) + fac = max(refspec[ta])/max(specs[ts,b]) + endif else begin + yrng = get_plotlims(refspec[ta],/minzero) + fac = 1. + endelse + tlab = imglab + ', Bar = '+string(b,"(i03)") + $ + ', Slice = '+string(fix(b/5),"(i02)") + if keyword_set(tweak) and sigmas[b] gt 0 then begin + ddwaves = reform(dwaves[b,*]) + ffwaves = reform(fwaves[b,*]) + li = where(ffwaves gt 0, nli) + if nli gt 0 then begin + ffwaves = ffwaves[li] + ddwaves = ddwaves[li] + endif else nli = 0 + tlab = tlab + ', RMS = '+string(sigmas[b],"(f5.3)") + $ + ', N = '+strn(nli) + endif else nli = 0 + plot,refwvl,refspec/fac,thick=th,charthi=th, title=tlab, $ + xthick=th,xtitle='Wavelength (A)',xrange=[minwav,maxwav],/xs, $ + ythick=th,ytitle='Flux',yrange=yrng,/ys + oplot,waves,specs[*,b],color=pcolor,thick=1.0 + oplot,[wavgood0,wavgood0],!y.crange,color=colordex('green'),linesty=1,thick=th + oplot,[wavgood1,wavgood1],!y.crange,color=colordex('green'),linesty=1,thick=th + oplot,[wavall0,wavall0],!y.crange,color=colordex('orange'),linesty=1,thick=th + oplot,[wavall1,wavall1],!y.crange,color=colordex('orange'),linesty=1,thick=th + oplot,[wavmid,wavmid],!y.crange,color=colordex('green'),linesty=1,thick=th + oplot,[cwave,cwave],!y.crange,color=colordex('black'),linesty=1,thick=th + ; + ; plot lines used for RMS calc + if keyword_set(tweak) and nli gt 0 then begin + ffwaves = ffwaves[li] + ddwaves = ddwaves[li] + for i=0,nli-1 do $ + oplot,[ffwaves[i],ffwaves[i]],!y.crange,color=rcolor + endif + kcwi_legend,[atnam+' Atlas',fittype+' Arc Fit'],linesty=[0,0],charthi=th, $ + color=[colordex('black'),pcolor], $ + thick=[th,th],/clear,clr_color=!p.background + ; + if not keyword_set(plot_file) then begin + read,'Next? (Q-quit plotting, - next): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then $ + display = (1 eq 0) + endif + ; + ; plot residuals + if keyword_set(tweak) and display then begin + if nli gt 0 then begin + residrng = get_plotlims(ddwaves) + resrng = [residrng[0]<(-0.2),residrng[1]>0.2] + endif else begin + ffwaves = fltarr(2) + ddwaves = fltarr(2) + resrng = [-0.2, 0.2] + endelse + plot,ffwaves,ddwaves,psym=4,charthi=th,thick=th, $ + xthick=th,xtitle='Wave(Ang)',xrange=[minwav,maxwav],/xs, $ + ythick=th,ytitle='Resid(Ang)',yrange=resrng,/ys, $ + title=tlab,/nodata + oplot,!x.crange,[0,0] + oplot,[wavgood0,wavgood0],!y.crange,color=colordex('green'),linesty=1,thick=th + oplot,[wavgood1,wavgood1],!y.crange,color=colordex('green'),linesty=1,thick=th + oplot,[wavall0,wavall0],!y.crange,color=colordex('orange'),linesty=1,thick=th + oplot,[wavall1,wavall1],!y.crange,color=colordex('orange'),linesty=1,thick=th + oplot,[wavmid,wavmid],!y.crange,color=colordex('green'),linesty=1,thick=th + oplot,[cwave,cwave],!y.crange,color=colordex('black'),linesty=1,thick=th + if nli gt 0 then begin + oplot,ffwaves,ddwaves,psym=4,thick=th + oplot,!x.crange,[sigmas[b],sigmas[b]],linesty=2 + oplot,!x.crange,-[sigmas[b],sigmas[b]],linesty=2 + endif + kcwi_legend,['All','Good','CWAVE'],linesty=[1,1,1],charthi=th, $ + color=[colordex('orange'),colordex('green'),colordex('black')], $ + thick=[th,th,th],/clear,clr_color=!p.background + ; + if not keyword_set(plot_file) then begin + read,'Next? (Q-quit plotting, - next): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then $ + display = (1 eq 0) + endif + endif ; are we plotting residuals? + endif ; display or keyword_set(plot_file) +endfor ;b + +if keyword_set(plot_file) then $ + psclose $ +else wdelete,1 +; +!p.multi=0 +; +return +end ; kcwi_plot_arcfits diff --git a/pro/kcwi_print_cfgs.pro b/pro/kcwi_print_cfgs.pro new file mode 100644 index 0000000..a447cdd --- /dev/null +++ b/pro/kcwi_print_cfgs.pro @@ -0,0 +1,125 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_PRINT_CFGS +; +; PURPOSE: +; This function prints a summary of the configurations passed, one +; line per image. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_PRINT_CFGS,Kcfg +; +; INPUTS: +; Kcfg - An array of struct KCWI_CFG. +; +; OUTPUTS: +; imsum - image summary (string) +; +; KEYWORDS: +; header - set to print headings for the columns +; silent - just return string, do not print +; outfile - filename to print to +; help - print help info for this routine +; +; PROCEDURE: +; Prints a summary allowing comparison of configurations of each image. +; +; EXAMPLE: +; Read in the stage one processed image data headers in directory +; 'redux' and return an array of struct KCWI_CFG. Find all the +; continuum flats and print their configuration summary. +; +; KCFG = KCWI_READ_CFGS('redux',filespec='*_int.fits') +; KCWI_PRINT_CFG, KCFG +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-JUL-08 Initial version +; 2013-NOV-13 Added outfile keyword +;- +pro kcwi_print_cfgs,kcfg,imsum, $ + header=header,silent=silent,outfile=outfile,help=help + ; + ; setup + pre = 'KCWI_PRINT_CFGS' + imsum = '' + ; + ; help request + if keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', Kcfg, Imsum' + print,pre+': Info Keywords: /header, /silent, outfile=' + return + endif + ; + ; check inputs + if kcwi_verify_cfg(kcfg) ne 0 then return + ; + ; outfile? + if keyword_set(outfile) then begin + filestamp,outfile,/arch + openw,ol,outfile,/get_lun + printf,ol,'# '+pre+' '+systime(0) + printf,ol,'# R = CCD Readout Speed : 0 - slow, 1 - fast' + printf,ol,'# G = Gain multiplier : 10, 5, 2, 1' + printf,ol,'# SSM = Sky, Shuffle, Mask: 0 - no, 1 - yes' + printf,ol,'# #/ N Imno Bin AMPS R G SSM IFU GRAT FILT Cwave JDobs Expt Type Imno RA Dec PA Air Object' + endif + ; + ; header? + if keyword_set(header) and not keyword_set(silent) then begin + print,' R = CCD Readout Speed : 0 - slow, 1 - fast, G = Gain multiplier' + print,' G = Gain multiplier : 10, 5, 2, 1' + print,' SSM = Sky, Shuffle, Mask: 0 - no, 1 - yes' + print,' #/ N Imno Bin AMPS R G SSM IFU GRAT FILT Cwave JDobs Expt Type Imno RA Dec PA Air Object' + endif + ; + ; current date + cdate = 0.d0 + ; + ; loop over elements + n = n_elements(kcfg) + for i=0,n-1l do begin + ; + ; prepare summary + imsum = string(i+1,'/',n,kcfg[i].imgnum, $ + kcfg[i].xbinsize,kcfg[i].ybinsize, $ + strtrim(kcfg[i].ampmode,2),kcfg[i].ccdmode, $ + kcfg[i].gainmul, $ + kcfg[i].skyobs,kcfg[i].shuffmod,kcfg[i].nasmask, $ + strtrim(kcfg[i].ifunam,2), $ + strtrim(kcfg[i].bgratnam,2), $ + strtrim(kcfg[i].bfiltnam,2), $ + kcfg[i].cwave,kcfg[i].juliandate, $ + kcfg[i].exptime,strtrim(kcfg[i].imgtype,2)+strtrim(kcfg[i].lampname,2), $ + kcfg[i].imgnum,kcfg[i].ra,kcfg[i].dec,kcfg[i].rotposn, $ + kcfg[i].airmass, $ + format='(i4,a1,i4,i7,2i2,1x,a-5,i1,1x,i2,1x,3i1,1x,a-3,1x,a-4,1x,a-4,1x,f8.1,f12.3,f7.1,1x,a-11,i7,2f13.8,2x,f7.2,f7.3)') + ; + ; add object info + if strpos(kcfg[i].imgtype,'object') ge 0 then begin + imsum = imsum + string(strtrim(kcfg[i].targname,2),form='(2x,a)') + endif + if not keyword_set(silent) then print,imsum + if keyword_set(outfile) then begin + if i gt 0 then $ + deljd = kcfg[i].juliandate - kcfg[i-1].juliandate $ + else deljd = 1.0 + if deljd gt 0.25 and kcfg[i].juliandate-cdate gt 0.75 then begin + cdate = kcfg[i].juliandate + caldat,long(cdate),month,day,year + printf,ol,'# Run: ',year-2000.,month,day, $ + format='(a,i02,i02,i02)' + endif + printf,ol,imsum + endif + endfor + if keyword_set(outfile) then free_lun,ol + ; + return +end diff --git a/pro/kcwi_print_info.pro b/pro/kcwi_print_info.pro new file mode 100644 index 0000000..43907ba --- /dev/null +++ b/pro/kcwi_print_info.pro @@ -0,0 +1,106 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_PRINT_INFO +; +; PURPOSE: +; This function prints info to the screen and to logfile +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_PRINT_INFO,Ppar,Pre,Msg,Var1,Var2... +; +; INPUTS: +; Ppar - array of struct KCWI_PPAR pipeline parameters +; Pre - a prefix string containing the name of the calling procedure +; Msg - the message to print +; Var1 - the variables to print +; +; KEYWORDS: +; ERROR - Set if message is an error message +; WARNING - Set if message is a warning message +; INFO - Set if message is only for information, can also set to +; level: will output if ppar.verbose ge level +; FORMAT - Set to format string for variables +; +; PROCEDURE: +; Constructs an output string using passed variables and format +; string. All warnings and errors are printed. Info messages are +; not printed unless their priority is less than or equal to the +; verbose level in ppar.verbose. Thus, if ppar.verbose = 0, a +; standard info message is not printed. By raising the priority +; rank, this requires the verbosity to be higher before it is printed. +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-SEP-16 Initial version +; 2013-NOV-21 Refinement of print rank logic +;- +pro kcwi_print_info,ppar,ipre,msg,v1,v2,v3,v4,v5,v6,v7,v8,v9, $ + error=error, warning=warning, info=info, format=format + ; + ; check inputs + ; ppar struct + if kcwi_verify_ppar(ppar) ne 0 then $ + ppar = { kcwi_ppar } + ; pre string + if size(ipre,/type) ne 7 then begin + print,'KCWI_PRINT_INFO: Error - input pre string error, returning' + return + endif + ; + ; make pre string standard length (20) + pre = ipre + while strlen(pre) lt 20 do $ + pre = pre + ' ' + ; + ; check keywords + if keyword_set(info) then $ + prank = info $ + else prank = 1 ; default rank + ; + ; default label (info) plus rank + ostr = pre+' [INFO '+strtrim(string(prank),2)+' ] '+msg + if keyword_set(warning) then begin + prank = 0 + ostr = pre+' [WARNING] '+msg + endif + if keyword_set(error) then begin + prank = 0 + ostr = pre+' [ERROR ] '+msg + endif + ; + ; check input variables + np = n_params(0) - 3 + if np gt 0 then begin + ostr = ostr + ': ' + case np of + 1: ostr = string(ostr,v1,format=format) + 2: ostr = string(ostr,v1,v2,format=format) + 3: ostr = string(ostr,v1,v2,v3,format=format) + 4: ostr = string(ostr,v1,v2,v3,v4,format=format) + 5: ostr = string(ostr,v1,v2,v3,v4,v5,format=format) + 6: ostr = string(ostr,v1,v2,v3,v4,v5,v6,format=format) + 7: ostr = string(ostr,v1,v2,v3,v4,v5,v6,v7,format=format) + 8: ostr = string(ostr,v1,v2,v3,v4,v5,v6,v7,v8,format=format) + 9: ostr = string(ostr,v1,v2,v3,v4,v5,v6,v7,v8,v9,format=format) + else: print,'KCWI_PRINT_INFO [ERROR] variable overflow' + endcase + endif + ; + ; check logging first + if ppar.loglun gt 0 then begin + printf,ppar.loglun,ostr + flush,ppar.loglun + endif + ; + ; check verbosity + if prank le ppar.verbose then $ + print,ostr + ; + return +end diff --git a/pro/kcwi_read_atlas.pro b/pro/kcwi_read_atlas.pro new file mode 100644 index 0000000..cf7871e --- /dev/null +++ b/pro/kcwi_read_atlas.pro @@ -0,0 +1,79 @@ +; +; Copyright (c) 2014, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_READ_ATLAS +; +; PURPOSE: +; Read the atlas spectrum and convolve to nominal KCWI resolution +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_READ_ATLAS, Kgeom, Ppar, Refspec, Refwave, Refdisp +; +; INPUTS: +; Kgeom - KCWI_GEOM struct from KCWI_TRACE_CBARS and KCWI_EXTRACT_ARCS +; Ppar - KCWI_PPAR pipeline parameter struct +; +; OUTPUTS: +; Refspec - Atlas reference spectrum +; Refwave - Atlas reference spectrum wavelengths +; Refdisp - Atlas reference spectrum dispersion in Ang/px +; +; INPUT KEYWORDS: +; +; SIDE EFFECTS: +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill +; 2014-SEP-18 Initial Revision +;- +; +pro kcwi_read_atlas, kgeom, ppar, refspec, refwave, refdisp + +pre = 'KCWI_READ_ATLAS' +q='' +; +; init +refspec = -1. +refwave = -1. +refdisp = -1. +; +; check inputs +if kcwi_verify_geom(kgeom,/init) ne 0 then return +if kcwi_verify_ppar(ppar,/init) ne 0 then return +; +; canonical resolution? +sigma = kgeom.atsig +; +; check if file is available +if not file_test(kgeom.refspec,/read,/regular) then begin + kcwi_print_info,ppar,pre,'Atlas spectrum file not found',kgeom.refspec,$ + format='(a,a)',/error + return +endif +; +; report the read +kcwi_print_info,ppar,pre,'Reading atlas spectrum in',kgeom.refspec, $ + format='(a,1x,a)' +; +; load the reference atlas spectrum. +rdfits1dspec,kgeom.refspec,refwave,atlas, $ + wavezero=refw0, deltawave=refdisp, refpix=refpix +refspec = atlas>0 +; +; we want to degrade this spectrum to the instrument resolution +xx = findgen(99)-50.0d +gaus = gaussian(xx,[1.0,0.0,sigma]) +gaus /= total(gaus) +refspec = convolve(refspec,gaus) +; +return +end ; kcwi_read_atlas diff --git a/pro/kcwi_read_cfg.pro b/pro/kcwi_read_cfg.pro new file mode 100644 index 0000000..c276afc --- /dev/null +++ b/pro/kcwi_read_cfg.pro @@ -0,0 +1,193 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_READ_CFG +; +; PURPOSE: +; This function reads the header from a KCWI image file and returns +; a kcwi_cfg structure. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; Result = KCWI_READ_CFG( OBSFNAME ) +; +; INPUTS: +; obsfname- filename of KCWI image file +; +; KEYWORDS: +; VERBOSE - set this to get extra screen output +; +; RETURNS: +; KCWI CFG struct (as defined in kcwi_cfg__define.pro) +; +; PROCEDURE: +; Reads in fits header and populates the KCWI_CFG struct. +; +; EXAMPLE: +; +; This will read in the file 'm82.fits' and return a kcwi_cfg struct: +; +; m82cfg = KCWI_READ_CFG('m82.fits') +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-MAY-03 Initial version +; 2014-SEP-12 Put in code to verify header +;- +function kcwi_read_cfg,obsfname,verbose=verbose +; +; initialize + pre = 'KCWI_READ_CFG' + A = {kcwi_cfg} ; get blank parameter struct + cfg = struct_init(A) ; initialize it +; +; get filename if not passed as a parameter + if n_params(0) lt 1 then begin + obsfname = '' + read,'Enter KCWI image FITS file name: ',obsfname + endif +; +; check file + fi = file_info(obsfname) + if not fi.exists or not fi.read or not fi.regular then begin + print,pre+': Error - file not accessible: ',obsfname + return,cfg + endif +; +; read header of KCWI image FITS file + hdr = headfits(obsfname) +; +; verify header + if size(hdr,/type) ne 7 or size(hdr,/dim) le 10 then begin + print,pre+': Error - bad header: ',obsfname + return,cfg + endif +; +; get parameter struct tags + keys = tag_names(cfg) + nkeys = n_elements(keys) + stopky = where(strcmp(keys,'JULIANDATE') eq 1)>0 change how things are loaded? +; > add some headers to the data. + +pro KCWI_REVERSE_GEOM, kgeom, ppar, degree=degree + +pre = "KCWI_REVERSE_GEOM" + +; Check structs +if kcwi_verify_geom(kgeom,/init) ne 0 then return +if kcwi_verify_ppar(ppar) ne 0 then begin + ppar = {kcwi_ppar} +endif + +; read in quantities from the geometry. +; trying to avoid using these inside of the code proper + imno = kgeom.cbarsimgnum + arcno = kgeom.arcimgnum + nx = kgeom.nx + ny = kgeom.ny + ypad = kgeom.ypad + nasmask = kgeom.nasmask + trimy0 = kgeom.trimy0 + trimy1 = kgeom.trimy1 + kwx = kgeom.kwx + kwy = kgeom.kwy + xi = kgeom.xi + yi = kgeom.yi + xw = kgeom.xw + yw = kgeom.yw + slice = kgeom.slice + refoutx = kgeom.refoutx + x0out = kgeom.x0out + dwout = kgeom.dwout + wave0out = kgeom.wave0out + ldeg = kgeom.lastdegree + + xpad = 3 + +; process the degree, default to 3. + if n_elements(degree) eq 0 then degree = ldeg + if degree lt 2 or degree gt 4 then degree = 3 + + outfile = kcwi_get_imname(ppar,kgeom.cbarsimgnum,"_wavemap",/reduced) + + x = dindgen(nx) + y = dindgen(ny); ypad + onex = x-x+1.00000d + oney = y-y+1.000000d + + xx = x#oney + yy = onex#y + + reverse_image = xx-xx-10 + + ; loop over slices + for s=0, 23 do begin + qs = where(slice eq s and xi gt 0 and finite(xw) and finite(yw), nqs) + if nqs eq 0 then begin + kcwi_print_info,ppar,pre, $ + "Sorry, no points in this slice. Confused. Exiting.",/error + return + endif + x0 = xi[qs] + y0 = yi[qs]+ypad + x1 = xw[qs] + y1 = yw[qs] + x0min = min(x0) + x0temp = x0-x0min+x0out +; stop + ; generate a mapping + polywarp,x1,y1,x0temp,y0,3,kx,ky,/double + ; reset the "in" values + xmm = minmax(x0) + qin = where(xx gt xmm[0]-xpad and xx lt xmm[1]+xpad and yy ge trimy0 and yy le trimy1, nqxin) + if nqxin eq 0 then begin + kcwi_print_info,ppar,pre, $ + "Sorry, no suitable points found. Confused. Exiting.",/error + return + endif + xin = xx[qin]-x0min+x0out + yin = yy[qin]+ypad + kcwi_poly_map,xin,yin,kx,ky,xout,yout + + ; set the pixel values to the wavelengths. + reverse_image[xin-x0out+x0min,yin-ypad] = yout*dwout+wave0out + + endfor + + ; + ; get header + hdr = headfits(kgeom.cbarsfname) + ; + ; update header + ; + ; spatial scale and zero point + sxaddpar,hdr,'BARSEP',kgeom.refdelx,' separation of bars (binned pix)' + sxaddpar,hdr,'BAR0',kgeom.x0out,' first bar pixel position' + ; + ; wavelength ranges + sxaddpar,hdr, 'WAVALL0', kgeom.waveall0, ' Low inclusive wavelength' + sxaddpar,hdr, 'WAVALL1', kgeom.waveall1, ' High inclusive wavelength' + sxaddpar,hdr, 'WAVGOOD0',kgeom.wavegood0, ' Low good wavelength' + sxaddpar,hdr, 'WAVGOOD1',kgeom.wavegood1, ' High good wavelength' + sxaddpar,hdr, 'WAVMID',kgeom.wavemid, ' middle wavelength' + ; + ; wavelength solution RMS + sxaddpar,hdr,'AVWVSIG',kgeom.avewavesig,' Avg. bar wave sigma (Ang)' + sxaddpar,hdr,'SDWVSIG',kgeom.stdevwavesig,' Stdev. bar wave sigma (Ang)' + ; + ; geometry solution RMS + xmo = moment(kgeom.xrsd,/nan) + ymo = moment(kgeom.yrsd,/nan) + sxaddpar,hdr, 'GEOXGSG', xmo[0], ' Global geometry X sigma (pix)' + sxaddpar,hdr, 'GEOYGSG', ymo[0], ' Global geometry Y sigma (pix)' + ; + ; pixel scales + sxaddpar,hdr,'PXSCL', kgeom.pxscl*kgeom.xbinsize,' Pixel scale along slice' + sxaddpar,hdr,'SLSCL', kgeom.slscl,' Pixel scale purpendicular to slices' + ; + ; geometry origins + sxaddpar,hdr, 'CBARSFL', kgeom.cbarsfname,' Continuum bars image' + sxaddpar,hdr, 'ARCFL', kgeom.arcfname, ' Arc image' + sxaddpar,hdr, 'CBARSNO', kgeom.cbarsimgnum,' Continuum bars image number' + sxaddpar,hdr, 'ARCNO', kgeom.arcimgnum, ' Arc image number' + sxaddpar,hdr, 'GEOMFL', kgeom.geomfile,' Geometry file' + + ; write the file + kcwi_print_info,ppar,pre,"Writing",outfile,/info,format='(a,1x,a)' + mwrfits, reverse_image, outfile, hdr,/create,/iscale + kcwi_print_info,ppar,pre,"Generated reverse map.",/info +end diff --git a/pro/kcwi_set_dgeom.pro b/pro/kcwi_set_dgeom.pro new file mode 100644 index 0000000..a112c71 --- /dev/null +++ b/pro/kcwi_set_dgeom.pro @@ -0,0 +1,154 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_SET_DGEOM +; +; PURPOSE: +; This procedure uses the input KCWI_CFG struct to set the basic +; parameters in the KCWI_DGEOM struct. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_SET_DGEOM, Kdgeom, iKcfg +; +; INPUTS: +; Kdgeom - Input KCWI_DGEOM struct. +; iKcfg - Input KCWI_CFG struct for a given observation. +; Ppar - Input KCWI_PPAR struct. +; +; KEYWORDS: +; None +; +; OUTPUTS: +; None +; +; SIDE EFFECTS: +; Sets the following tags in the KCWI_DGEOM struct according to the +; configuration settings in KCWI_CFG. +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2016-NOV-09 Initial version +;- +pro kcwi_set_dgeom,kdgeom,ikcfg,ppar,help=help + ; + ; setup + pre = 'KCWI_SET_DGEOM' + ; + ; help request + if keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', Kdgeom, Kcfg, Ppar' + return + endif + ; + ; verify Kdgeom + if kcwi_verify_geom(kdgeom,/init) ne 0 then return + ; + ; verify Kcfg + if kcwi_verify_cfg(ikcfg,/silent) ne 0 then begin + print,pre+': Error - malformed KCWI_CFG struct' + return + endif + ; + ; verify Ppar + psz = size(ppar) + if psz[2] eq 8 then begin + if ppar.initialized ne 1 then begin + print,pre+': Error - KCWI_PPAR struct not initialized.' + return + endif + endif else begin + print,pre+': Error - malformed KCWI_PPAR struct' + return + endelse + ; + ; take singleton of KCWI_CFG + kcfg = ikcfg[0] + ; + ; check image type + if strtrim(strupcase(kcfg.imgtype),2) ne 'ARC' then begin + kcwi_print_info,ppar,pre,'arc images are the direct geom reference files, this file is of type',kcfg.imgtype,/error + return + endif + ; + ; get output dgeom file name + odir = ppar.reddir + kdgeom.geomfile = ppar.reddir + $ + strmid(kcfg.obsfname,0,strpos(kcfg.obsfname,'_int')) + '_dgeom.fits' + ; + ; set basic configuration parameters + kdgeom.arcimgnum = kcfg.imgnum + kdgeom.ifunum = kcfg.ifunum + kdgeom.ifunam = kcfg.ifunam + kdgeom.filter = kcfg.filter + kdgeom.filtnum = kcfg.filtnum + kdgeom.campos = kcfg.campos + kdgeom.camang = kcfg.camang + kdgeom.xbinsize = kcfg.xbinsize + kdgeom.ybinsize = kcfg.ybinsize + kdgeom.nx = kcfg.naxis1 + kdgeom.ny = kcfg.naxis2 + ; + ; get noise model + rdnoise = 0. + ; + ; sum over amp inputs + switch kcfg.nvidinp of + 4: rdnoise = rdnoise + kcfg.biasrn4 + 3: rdnoise = rdnoise + kcfg.biasrn3 + 2: rdnoise = rdnoise + kcfg.biasrn2 + 1: rdnoise = rdnoise + kcfg.biasrn1 + endswitch + ; + ; take average + rdnoise /= float(kcfg.nvidinp) + kdgeom.rdnoise = rdnoise + ; + ; spatial scales + kdgeom.pxscl = 0.00004048d0 ; deg/unbinned pixel + kdgeom.slscl = 0.00037718d0 ; deg/slice, Large slicer + ; + ; set IFU-specific params + case kdgeom.ifunum of + 1: begin ; Large slicer + kdgeom.do_gauss = 0 + kdgeom.minpix = 6 / kdgeom.ybinsize + end + 2: begin ; Medium slicer + if kdgeom.ybinsize eq 2 then begin ; 2x2 binning + kdgeom.do_gauss = 1 + kdgeom.minpix = 3 + endif else begin ; 1x1 binning + kdgeom.do_gauss = 0 + kdgeom.minpix = 6 + endelse + kdgeom.slscl = kdgeom.slscl/2.d0 + end + 3: begin ; Small slicer + kdgeom.do_gauss = 1 + if kdgeom.ybinsize eq 2 then $ + kdgeom.minpix = 1 $ + else kdgeom.minpix = 2 + kdgeom.slscl = kdgeom.slscl/4.d0 + end + else: begin + kcwi_print_info,ppar,pre,'Unknown slicer ID', $ + kdgeom.ifunum,/error + return + end + endcase + ; + ; log our update of the geom struct + kdgeom.progid = pre + kdgeom.timestamp = systime(1) + ; + return +end diff --git a/pro/kcwi_set_geom.pro b/pro/kcwi_set_geom.pro new file mode 100644 index 0000000..f995c4f --- /dev/null +++ b/pro/kcwi_set_geom.pro @@ -0,0 +1,302 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_SET_GEOM +; +; PURPOSE: +; This procedure uses the input KCWI_CFG struct to set the basic +; parameters in the KCWI_GEOM struct. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_SET_GEOM, Kgeom, iKcfg +; +; INPUTS: +; Kgeom - Input KCWI_GEOM struct. +; iKcfg - Input KCWI_CFG struct for a given observation. +; Ppar - Input KCWI_PPAR struct. +; +; KEYWORDS: +; Atlas - Arc atlas fits file (string, eg: 'fear.fits') +; Atname - Arc atlas name (string, eg: 'FeAr') +; +; OUTPUTS: +; None +; +; SIDE EFFECTS: +; Sets the following tags in the KCWI_GEOM struct according to the +; configuration settings in KCWI_CFG. +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-AUG-13 Initial version +; 2016-JUN-12 Added KCWI gratings BH2,BH3, BM, BL +; 2016-JUL-01 Added atlas, atname keywords +;- +pro kcwi_set_geom,kgeom,ikcfg,ppar,atlas=atlas,atname=atname, help=help + ; + ; setup + pre = 'KCWI_SET_GEOM' + ; + ; help request + if keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', Kgeom, Kcfg, Ppar' + return + endif + ; + ; verify Kgeom + ksz = size(kgeom) + if ksz[2] eq 8 then begin + if kgeom.initialized ne 1 then begin + print,pre+': Error - KCWI_GEOM struct not initialized.' + return + endif + endif else begin + print,pre+': Error - malformed KCWI_GEOM struct' + return + endelse + ; + ; verify Kcfg + if kcwi_verify_cfg(ikcfg,/silent) ne 0 then begin + print,pre+': Error - malformed KCWI_CFG struct' + return + endif + ; + ; verify Ppar + psz = size(ppar) + if psz[2] eq 8 then begin + if ppar.initialized ne 1 then begin + print,pre+': Error - KCWI_PPAR struct not initialized.' + return + endif + endif else begin + print,pre+': Error - malformed KCWI_PPAR struct' + return + endelse + ; + ; take singleton of KCWI_CFG + kcfg = ikcfg[0] + ; + ; check image type + if strtrim(strupcase(kcfg.imgtype),2) ne 'CBARS' then begin + kcwi_print_info,ppar,pre,'cbars images are the geom reference files, this file is of type',kcfg.imgtype,/error + return + endif + ; + ; set up bar offsets + baroffs = [[[-220., 20., -220., 0., -210., -14., -210., -30., -210., -40., -200., 0., $ ; BH3 Lar + -220., -20., -170., -10., -180., 0., -190., 20., -190., 40., -190., 70.], $ + [-260., -20., -250., -30., -240., -40., -230., -50., -220., -50., -200., 0., $ ; BH3 Med + -215., -10., -160., 5., -160., 25., -160., 50., -150., 80., -140., 118.], $ + [-286., -40., -270., -50., -260., -50., -240., -50., -225., -60., -210., 0., $ ; BH3 Sma + -215., -7., -151., 15., -150., 40., -140., 70., -130., 100., -120., 140.]], $ + [[ -80., 150., -120., 90., -160., 40., -190., -4., -200., -40., -220., 0., $ ; BL Lar + -240., -20., -180., 8., -170., 40., -160., 90., -120., 150., -80., 240.], $ + [-140., 106., -175., 50., -200., 10., -215., -20., -220., -50., -225., 0., $ ; BL Med + -240., -8., -160., 25., -150., 70., -120., 130., -80., 200., -30., 290.], $ + [-160., 80., -195., 35., -215., -4., -220., -30., -230., -50., -230., 0., $ ; BL Sma + -240., -3., -160., 35., -140., 85., -105., 150., -60., 220., -4., 315.]], $ + [[-220., 20., -220., 0., -210., -16., -210., 30., -200., -40., -200., 0., $ ; BH2 Lar + -220., -20., -170., -10., -180., 0., -190., 20., -190., 40., -190., 60.], $ + [-260., -25., -250., -30., -240., -40., -230., -40., -215., -50., -200., 0., $ ; BH2 Med + -210., -10., -160., 5., -160., 25., -160., 50., -150., 75., -140., 110.], $ + [-280., -50., -270., -50., -250., -50., -240., -50., -220., -60., -205., 0., $ ; BH2 Sma + -210., -7., -150., 15., -150., 40., -140., 60., -130., 95., -120., 135.]], $ + [[-140., 100., -160., 60., -180., 20., -200., -15., -200., -40., -210., 0., $ ; BM Lar + -240., -20., -180., 0., -180., 20., -160., 60., -150., 100., -120., 170.], $ + [-180., 60., -200., 20., -210., -8., -220., -30., -220., -50., -220., 0., $ ; BM Med + -230., -10., -160., 20., -150., 50., -135., 100., -110., 150., -70., 220.], $ + [-210., 40., -220., 5., -230., -20., -230., -40., -230., -50., -220., 0., $ ; BM Sma + -230., -5., -150., 30., -140., 70., -120., 120., -90., 175., -50., 250.]]] + + ; + ; get output geom file name + odir = ppar.reddir + kgeom.geomfile = ppar.reddir + $ + strmid(kcfg.obsfname,0,strpos(kcfg.obsfname,'_int')) + '_geom.fits' + ; + ; set basic configuration parameters + kgeom.ifunum = kcfg.ifunum + kgeom.ifunam = kcfg.ifunam + kgeom.gratid = kcfg.gratid + kgeom.gratnum = kcfg.gratnum + kgeom.filter = kcfg.filter + kgeom.filtnum = kcfg.filtnum + kgeom.campos = kcfg.campos + kgeom.camang = kcfg.camang + kgeom.grenc = kcfg.grenc + kgeom.grangle = kcfg.grangle + kgeom.gratanom = kcfg.gratanom + kgeom.xbinsize = kcfg.xbinsize + kgeom.ybinsize = kcfg.ybinsize + kgeom.nx = kcfg.naxis1 + kgeom.ny = kcfg.naxis2 + kgeom.x0out = 30 / kgeom.xbinsize ; re-set in kcwi_trace_cbars + kgeom.goody0 = 10 + kgeom.goody1 = kgeom.ny - 10 + kgeom.trimy0 = 0 + kgeom.trimy1 = kgeom.ny + kgeom.ypad = 1400 / kgeom.ybinsize + kgeom.nasmask = kcfg.nasmask + kgeom.goodmy0 = kcfg.nsobjr0 + kgeom.goodmy1 = kcfg.nsobjr1 - 18 / kgeom.ybinsize + kgeom.trimmy0 = kcfg.nsobjr0 + kgeom.trimmy1 = kcfg.nsobjr1 - 18 / kgeom.ybinsize + ; + ; get noise model + rdnoise = 0. + ; + ; sum over amp inputs + switch kcfg.nvidinp of + 4: rdnoise = rdnoise + kcfg.biasrn4 + 3: rdnoise = rdnoise + kcfg.biasrn3 + 2: rdnoise = rdnoise + kcfg.biasrn2 + 1: rdnoise = rdnoise + kcfg.biasrn1 + endswitch + ; + ; take average + rdnoise /= float(kcfg.nvidinp) + kgeom.rdnoise = rdnoise + ; + ; wavelength numbers default from header + kgeom.cwave = kcfg.cwave + kgeom.wave0out = kcfg.wave0 + kgeom.wave1out = kcfg.wave1 + kgeom.dwout = kcfg.dwav + ; + ; reference spectrum: ppar value has top priority + if strlen(strtrim(ppar.atlas,2)) gt 0 then begin + kgeom.refspec = ppar.datdir+ppar.atlas + kgeom.refname = ppar.atlasname + endif else if keyword_set(atlas) then begin + kgeom.refspec = ppar.datdir+atlas + kgeom.refname = atname + endif else begin + kgeom.refspec = ppar.datdir+'fear.fits' + kgeom.refname = 'FeAr' + endelse + ; + ; default to no cc offsets + kgeom.ccoff = reform(baroffs[*,kgeom.ifunum-1,kgeom.gratnum-1]) + kgeom.ccoff = kgeom.ccoff/kgeom.ybinsize + ; + ; grating parameters BH1 + if strtrim(kcfg.gratid,2) eq 'BH1' then begin + kgeom.atsig = 2.5 + if kgeom.ifunum gt 2 then $ + kgeom.atsig = 1.5 + kgeom.ccwn = 360./kgeom.ybinsize + kgeom.rho = 3.751d + kgeom.adjang = 180.d + kgeom.lastdegree = 4 + kgeom.bclean = 1 + ; + ; output disperison + kgeom.dwout = 0.125 * float(kcfg.ybinsize) + endif + ; + ; grating parameters BH2 + if strtrim(kcfg.gratid,2) eq 'BH2' then begin + kgeom.atsig = 2.5 + if kgeom.ifunum gt 2 then $ + kgeom.atsig = 1.5 + kgeom.ccwn = 100. ;360./kgeom.ybinsize + kgeom.rho = 3.255d + kgeom.adjang = 180.d + kgeom.lastdegree = 4 + kgeom.bclean = 0 + ; + ; output disperison + kgeom.dwout = 0.125 * float(kcfg.ybinsize) + endif + ; + ; grating parameters BH3 + if strtrim(kcfg.gratid,2) eq 'BH3' then begin + kgeom.atsig = 2.5 + if kgeom.ifunum gt 2 then $ + kgeom.atsig = 1.5 + kgeom.ccwn = 100. ;360./kgeom.ybinsize + kgeom.rho = 2.80d + kgeom.adjang = 180.d + kgeom.lastdegree = 4 + kgeom.bclean = 0 + ; + ; output disperison + kgeom.dwout = 0.125 * float(kcfg.ybinsize) + endif + ; + ; grating parameters BM + if strtrim(kcfg.gratid,2) eq 'BM' then begin + kgeom.atsig = 4.0 + if kgeom.ifunum ge 2 then $ + kgeom.atsig = 2.0 + kgeom.ccwn = 75. ;260./kgeom.ybinsize + kgeom.rho = 1.900d + kgeom.adjang = 0.d + kgeom.lastdegree = 4 + kgeom.bclean = 0 + ; + ; output disperison + kgeom.dwout = 0.25 * float(kcfg.ybinsize) + endif + ; + ; grating parameters BL + if strtrim(kcfg.gratid,2) eq 'BL' then begin + kgeom.atsig = 14.0 + if kgeom.ifunum eq 2 then $ + kgeom.atsig = 10.0 + if kgeom.ifunum eq 3 then $ + kgeom.atsig = 7.0 + kgeom.ccwn = 75. ;320./kgeom.ybinsize + kgeom.rho = 0.870d + kgeom.adjang = 0.d + kgeom.lastdegree = 4 + kgeom.bclean = 1 + ; + ; output disperison + kgeom.dwout = 0.5 * float(kcfg.ybinsize) + endif + ; + ; spatial scales + kgeom.pxscl = 0.00004048d0 ; deg/unbinned pixel + kgeom.slscl = 0.00037718d0 ; deg/slice, Large slicer + if kcfg.ifunum eq 2 then begin + kgeom.slscl = kgeom.slscl/2.d0 + endif else if kcfg.ifunum eq 3 then begin + kgeom.slscl = kgeom.slscl/4.d0 + endif + ; + ; check central wavelength + if kgeom.cwave le 0. then begin + kcwi_print_info,ppar,pre,'No central wavelength found',/error + return + endif + ; + ; now check ppar values which override defaults + if ppar.dw gt 0. then $ + kgeom.dwout = ppar.dw + if ppar.wave0 gt 0. then $ + kgeom.wave0out = ppar.wave0 + if ppar.wave1 gt 0. then $ + kgeom.wave1out = ppar.wave1 + if ppar.cleancoeffs gt -1. then $ + kgeom.bclean = ppar.cleancoeffs + ; + ; print log of values + kcwi_print_info,ppar,pre,'Data cube output Disp (A/px), Wave0 (A): ', $ + kgeom.dwout,kgeom.wave0out,format='(a,f8.3,f9.2)' + ; + ; log our change of the kgeom struct + kgeom.progid = pre + kgeom.timestamp = systime(1) + ; + return +end diff --git a/pro/kcwi_solve_arcs.pro b/pro/kcwi_solve_arcs.pro new file mode 100644 index 0000000..c86d46f --- /dev/null +++ b/pro/kcwi_solve_arcs.pro @@ -0,0 +1,692 @@ +; +; Copyright (c) 2013-2017, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_SOLVE_ARCS +; +; PURPOSE: +; Solve the wavelength solutions for each bar of the arc spectrum +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_SOLVE_ARCS, Specs, Cntcoeff, Kgeom, Ppar +; +; INPUTS: +; Specs - a array of arc spectra produced by KCWI_EXTRACT_ARCS +; Cntcoeff - Central fit coefficients (from kcwi_fit_center.pro) +; Kgeom - KCWI_GEOM struct from KCWI_TRACE_CBARS and KCWI_EXTRACT_ARCS +; Ppar - KCWI_PPAR pipeline parameter struct +; +; INPUT KEYWORDS: +; PLOT_FILE - if set to a string, will be used to produce +; postscript output of diagnostic plots +; +; SIDE EFFECTS: +; Modifies KCWI_GEOM struct by updating the bar coeffiecients and +; wavelength fit order and calculating new control points that take +; into account the wavelength solution. +; NOTE: sets KGEOM.STATUS to 0 if fitting succeeded, otherwise sets to +; 1 or greater depending on reason for failure (see code below). +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Matt Matuszewski +; 2017-MAY-10 Initial Revision +;- +; +pro kcwi_solve_arcs, specs, cntcoeff, kgeom, ppar, plot_file=plot_file + +pre = 'KCWI_SOLVE_ARCS' +q='' +; +; check status of central fit +if kgeom.status gt 0 then begin + kcwi_print_info,ppar,pre,'Bad central solutions, cannot solve arcs',/error + return +endif +; +; do we want to display stuff? +display = (ppar.display gt 2) +; +; set up plots +if ppar.display ge 1 then begin + deepcolor + !p.background=colordex('white') + !p.color=colordex('black') + th=2.0 + si=1.5 +endif +; +; which image number +imgnum = kgeom.arcimgnum +; +; which slicer? +ifunum = kgeom.ifunum +ifunam = kgeom.ifunam +; +; which grating? +grating = kgeom.gratid +; +; which filter? +filter = kgeom.filter +; +; image label +imglab = 'Img # '+strn(imgnum)+' ('+kgeom.refname+') Sl: '+strtrim(ifunam,2)+ $ + ' Fl: '+strtrim(filter,2)+' Gr: '+strtrim(grating,2) +; +; is this N+S mask in? +nasmask = kgeom.nasmask +; +; central wavelength? +cwvl = kgeom.cwave +; +; canonical resolution in Angstroms? +resolution = kgeom.atsig +; +; log info +kcwi_print_info,ppar,pre,systime(0) +kcwi_print_info,ppar,pre,'img, grat, filt, nasmsk, cwave', $ + imgnum,grating,filter,nasmask,cwvl, $ + format='(a,i6,2x,a-8,2x,a-8,i4,f12.3)' +; +; last degree is for full ccd tweaked fits +degree = kgeom.lastdegree +; +; which is the reference bar? +refbar = kgeom.refbar +; +; input spectrum dimensions +specsz = size(specs,/dim) +; +; set up array with zero point in the center +x0 = specsz[0]/2 +;xvals = dindgen(specsz[0])-x0 +xvals = dindgen(specsz[0]) +; +; start with central fit coefficients +twkcoeff = cntcoeff +; +; keep track of final sigmas +sigmas = dblarr(120) +; +; keep track of individual bar fits +barstat = intarr(120) +barrej = intarr(120) +; +; keep track of observed versus atlas comparison +fwaves = fltarr(120,1000) +dwaves = fltarr(120,1000) +; +; keep track of reference wavelengths and matched pixel positions +rwaves = fltarr(120,1000) +xcents = fltarr(120,1000) +; +; x range of spectra to consider +if nasmask then begin + minrow = (1*specsz[0])/3 + maxrow = (2*specsz[0])/3 +endif else begin + lastrow = 50 + minrow = lastrow + maxrow = specsz[0]-lastrow-1 +endelse ; no nasmask +ftype = 'Std' +; +; find wavelength range +; and pascal shift coeffs +mnwvs = fltarr(120) +mxwvs = fltarr(120) +for b=0,119 do begin + twkcoeff[*,b] = pascal_shift(cntcoeff[*,b],x0,/silent) + waves = poly(xvals,twkcoeff[*,b]) + mnwvs[b] = min(waves) + mxwvs[b] = max(waves) +endfor +; +minwav = min(mnwvs) +maxwav = max(mxwvs) +; +; do a continuum subtraction +if nasmask then sectors = 6 else sectors = 16 +; +div = (maxrow-minrow)/sectors +; +xv = fltarr(sectors) +yv = fltarr(sectors) +for b = 0, 119 do begin + ; + ; get continuum point in each sector + for sec = 0,sectors-1 do begin + ; + ; We may need this for low-res gratings! + mn = min(specs[minrow+sec*div:minrow+(sec+1)*div,b],mi) + xv[sec] = mi+minrow+sec*div + ;mn = median(specs[minrow+sec*div:minrow+(sec+1)*div,b]) + ;xv[sec] = minrow+sec*div + div/2 + yv[sec] = mn + endfor + ; + ; fit sector sky points + res = poly_fit(xv,yv,3,yfit=yf) + ; + xvs = findgen(specsz[0]) + specs[*,b] -= res[0]+res[1]*xvs+res[2]*xvs^2+res[3]*xvs^3 + ; +endfor ; b +; +; use reference bar +b = refbar +; +; fill out the arrays we are working with. +subxvals = xvals[minrow:maxrow] +subyvals = smooth(reform(specs[minrow:maxrow,b]),3) +subwvals = poly(subxvals,twkcoeff[*,b]) +; +; find good peaks in object spectrum +;smooth_width = fix(resolution/abs(twkcoeff[1,b]))>4 ; in pixels +smooth_width = 4 +;peak_width = fix(smooth_width*2.5) ; for fitting peaks +peak_width = fix(resolution/abs(twkcoeff[1,b]))>4 ; in pixels +;slope_thresh = 0.7*smooth_width/resolution/100.0 ; more severe for object +slope_thresh = 0.7*smooth_width/2./100.0 ; more severe for object +kcwi_print_info,ppar,pre,'using a peak_width (px) of', $ + peak_width,format='(a,1x,i5)' +; +; set to get the most lines +ampl_thresh = 0. +spec_cent = findpeaks(subwvals,subyvals,smooth_width,slope_thresh, $ + ampl_thresh,peak_width,count=spec_npks,avsg=avwsg) +avwfwhm = avwsg * 2.355 +kcwi_print_info,ppar,pre,'starting with ',spec_npks,' lines.' +kcwi_print_info,ppar,pre,'avg line width (A)',avwsg,form='(a,1x,f9.3)' +kcwi_print_info,ppar,pre,'avg line FWHM (A)',avwfwhm,form='(a,1x,f9.3)' +; +; load the reference atlas spectrum. +kcwi_read_atlas,kgeom,ppar,refspec,refwvl,refdisp +; +; first trim to the grating bandpass (and a buffer of 100 Ang). +qwvs = where(refwvl gt minwav-100. and $ + refwvl lt maxwav+100.,nqwvs) +if nqwvs eq 0 then begin + kcwi_print_info,ppar,pre,'No wavelength overlap with atlas spectrum',/error + kgeom.status=5 + return +endif +refspec = refspec[qwvs] +refwvl = refwvl[qwvs] +kcwi_print_info,ppar,pre,'atlas disp (A/pix)',refdisp,form='(a,1x,f9.3)' +; +; set spectra fitting width +; BM, BL are well sampled, but other +; gratings need a bigger window +if strpos(grating,'BH') ge 0 or strpos(grating,'BM') ge 0 then $ + fwid = avwfwhm $ +else fwid = avwsg +; +; generate an Atlas line list +; +; setup wavelength fitting arrays +refxs = fltarr(spec_npks)-1. +refws = fltarr(spec_npks)-1. +refwd = fltarr(spec_npks)-1. +refof = fltarr(spec_npks)-1. +fi = 0 ; fit line index +; +; examine lines +for pp = 0,spec_npks-1 do begin + is_good = (1 eq 1) + !p.multi=[0,2,1] + ; + ; fit Atlas peak in wavelength + ffs = get_line_window(refwvl,refspec,spec_cent[pp],count=nffs) + if nffs gt 5 then begin + !quiet = 1 ; supress curvefit error messages + yf = gaussfit(refwvl[ffs],refspec[ffs],a,nterms=3) + !quiet = 0 ; turn error messages back on + pka = (where(refspec[ffs] eq max(refspec[ffs])))[0] + xoff = abs(refwvl[ffs[pka]] - a[1])/refdisp + woff = abs(spec_cent[pp]-a[1]) + wrat = a[2]/fwid + if woff gt 5. or xoff gt 1.5 or wrat gt 1.1 then $ + is_good = (1 eq 0) + at_pk_wl = a[1] + if is_good and display then begin + wvs = where(refwvl gt spec_cent[pp]-30. and $ + refwvl lt spec_cent[pp]+30.) + plot,refwvl[wvs],refspec[wvs],psym=10,title='Atlas', $ + xtitle='Wave(A)',/xs,ytitle='Flux' + oplot,refwvl[ffs],yf,color=colordex('G') + endif + endif else begin + is_good = (1 eq 0) + xoff = -1. + woff = -1. + wrat = 0. + endelse + ;print,'Atlas: ',is_good,nffs,xoff,woff,wrat + ; + ; fit Spec peak in X pixels + ffs = get_line_window(subwvals,subyvals,spec_cent[pp],count=nffs) + if nffs gt 5 then begin + !quiet = 1 ; supress curvefit error messages + yf = gaussfit(subxvals[ffs],subyvals[ffs],a,nterms=3) + !quiet = 0 ; turn error messages back on + pka = (where(subyvals[ffs] eq max(subyvals[ffs])))[0] + xoff = abs(subxvals[ffs[pka]] - a[1]) + xrat = a[2]/(avwsg/abs(twkcoeff[1,b])) + if xrat gt 1.25 or not finite(a[1]) or a[2] gt nffs then $ + is_good = (1 eq 0) + sp_pk_x = a[1] + if is_good and display then begin + wvs = where(subwvals gt spec_cent[pp]-30. and $ + subwvals lt spec_cent[pp]+30.) + plot,subxvals[wvs],subyvals[wvs],psym=10,title='Spec', $ + xtitle='Wave(A)',/xs,ytitle='Flux' + oplot,subxvals[ffs],yf,color=colordex('G') + endif + endif else begin + is_good = (1 eq 0) + xoff = -1. + xrat = 0. + endelse + ;print,'Spec : ',is_good,nffs,xoff,xrat + ;print,' ' + ; + ; are we good? + if is_good then begin + if display then begin + q = '' + read,' - skip, f - fit: ',q + if strupcase(strtrim(q,2)) eq 'F' then begin + refxs[fi] = sp_pk_x + refws[fi] = at_pk_wl + refwd[fi] = wrat + refof[fi] = woff + fi += 1 + endif else if strupcase(strtrim(q,2)) eq 'Q' then begin + display = (1 eq 0) + endif + endif else begin + refxs[fi] = sp_pk_x + refws[fi] = at_pk_wl + refwd[fi] = wrat + refof[fi] = woff + fi += 1 + endelse + endif ; are we good? +endfor ; examine lines +; +; get fit lines +good = where(refws gt 0., ngood) +refws = refws[good] +refxs = refxs[good] +refwd = refwd[good] +refof = refof[good] +; +; preliminary fit +newcoeff = poly_fit(refxs,refws,degree,yfit=fitw) +; +; get residual stats +diff = refws - fitw +ims, diff, rmn, rsg, rwgt +bad = where(rwgt le 0, nbad) +; +; clean outliers +if nbad gt 0 then begin + good = where(rwgt gt 0, ngood) + refws = refws[good] + refxs = refxs[good] + refwd = refwd[good] + refof = refof[good] + newcoeff = poly_fit(refxs,refws,degree,yfit=fitw) +endif else ngood = n_elements(refws) +diff = refws - fitw +; +; account for cleaning zero-point wavelength +mo = moment(diff) +diff = diff - mo[0] +; +; final RMS +mo = moment(diff) +sigmas[b] = sqrt(mo[1]) +; +; store for later plots +fwaves[b,0:(ngood-1)] = fitw +dwaves[b,0:(ngood-1)] = diff +rwaves[b,0:(ngood-1)] = refws +xcents[b,0:(ngood-1)] = refxs +barstat[b] = 0 +barrej[b] = nbad +; +; log stats +kcwi_print_info,ppar,pre,ftype+' Fit: Bar#,Npks,RMS,Coefs', $ + b,ngood,sigmas[b],newcoeff[0:degree], $ + format='(a,2i5,2x,f7.3,f9.2,f8.4,'+strn((degree-1)>1)+'g13.5)' +; +; store resulting coeffs +twkcoeff[*,b] = 0. +twkcoeff[0:degree,b] = newcoeff +; +; plot residuals +if ppar.display ge 1 then begin + !p.multi=0 + yran = get_plotlims(diff) + if abs(newcoeff[1]) gt 0.4 then $ + yrng = [yran[0]<(-0.4),yran[1]>0.4] $ + else yrng = [yran[0]<(-0.2),yran[1]>0.2] + plot,refws,diff,psym=4,charsi=si,charthi=th,thick=th, $ + title=imglab+' (bar '+strn(b)+' REF)', $ + xthick=th, xtitle='Wavelength (A)', $ + ythick=th, ytitle='Atlas - Spec (A)', yrange=yrng, /ys + oplot,!x.crange,[0,0], linesty=2, color=colordex('blue') + oplot,!x.crange,[sigmas[b],sigmas[b]],linesty=1,color=colordex('blue') + oplot,!x.crange,-[sigmas[b],sigmas[b]],linesty=1,color=colordex('blue') + kcwi_legend,["RMS = "+string(sigmas[b],form='(f7.3)')+' Ang', $ + 'NMatch = '+strn(ngood), $ + 'NRej = '+strn(nbad)],charsi=si,charthi=th +endif +if display then $ + read,'next: ',q +; +; write out atlas line list +atllist = ppar.reddir+ppar.froot+string(imgnum,'(i0'+strn(ppar.fdigits)+')')+'_atlas.txt' +openw,al,atllist,/get_lun +printf,al,'# '+pre+': Atlas lines (Ang) printed on '+systime(0) +printf,al,'# Atlas filename: '+kgeom.refspec +printf,al,'# Cont Bars filename: '+kgeom.cbarsfname +printf,al,'# Arc filename: '+kgeom.arcfname +printf,al,'# Reference bar: '+string(refbar) +printf,al,'# RMS of fit: '+string(rsg,form='(f9.3)') +printf,al,'# Avg Spec Width (A): '+string(avwsg,form='(f9.3)') +printf,al,'# Atlas disp (A/px): '+string(refdisp,form='(f9.3)') +printf,al,'# RefWave SpecPix FitResid LWidRat LOffRat' +for j=0,ngood-1 do printf,al,refws[j],refxs[j],diff[j], $ + refwd[j],refof[j],format='(5f12.3)' +free_lun,al +kcwi_print_info,ppar,pre,'Atlas lines written to',atllist,format='(a,a)' +; +; now loop over other bars +for b = 0,119 do begin + ; + ; setup wavelength fitting arrays + fitxs = refxs - refxs + ; + ; skip refbar (already done) + if b ne refbar then begin + ; + ; fill out the arrays we are working with. + subxvals = xvals[minrow:maxrow] + subyvals = smooth(reform(specs[minrow:maxrow,b]),3) + subwvals = poly(subxvals,twkcoeff[*,b]) + ; + ; loop over reference lines + for pp = 0,n_elements(refws)-1 do begin + ; + ; fit Spec peak + ffs = get_line_window(subwvals,subyvals,refws[pp], $ + count=nffs) + if nffs gt 5 then begin + !quiet = 1 ; supress curvefit error msgs + yf = gaussfit(subxvals[ffs],subyvals[ffs],a, $ + nterms=3) + !quiet = 0 ; turn error messages back on + if finite(a[1]) and a[2] lt nffs then $ + fitxs[pp] = a[1] + endif + endfor ; loop over reference lines + ; + ; get fit lines + good = where(fitxs gt 0., ngood) + fitws = refws[good] + fitxs = fitxs[good] + ; + ; preliminary fit + newcoeff = poly_fit(fitxs,fitws,degree,yfit=fitw) + ; + ; get residual stats + diff = fitws - fitw + ims, diff, rmn, rsg, rwgt + bad = where(rwgt le 0, nbad) + ; + ; clean outliers + if nbad gt 0 then begin + good = where(rwgt gt 0, ngood) + fitws = fitws[good] + fitxs = fitxs[good] + newcoeff = poly_fit(fitxs,fitws,degree,yfit=fitw) + endif else ngood = n_elements(fitws) + diff = fitws - fitw + ; + ; account for cleaning zero-point wavelength + mo = moment(diff) + diff = diff - mo[0] + ; + ; final RMS + mo = moment(diff) + sigmas[b] = sqrt(mo[1]) + ; + ; store for later plots + fwaves[b,0:(ngood-1)] = fitw + dwaves[b,0:(ngood-1)] = diff + rwaves[b,0:(ngood-1)] = fitws + xcents[b,0:(ngood-1)] = fitxs + ; + ; log stats + kcwi_print_info,ppar,pre,ftype+' Fit: Bar#,Npks,RMS,Coefs', $ + b,ngood,sigmas[b],newcoeff[0:degree], $ + format='(a,2i5,2x,f7.3,f9.2,f8.4,'+strn((degree-1)>1)+'g13.5)' + ; + ; store resulting coeffs + twkcoeff[*,b] = 0. + twkcoeff[0:degree,b] = newcoeff + ; + ; plot residuals + if display then begin + yran = get_plotlims(diff) + if abs(newcoeff[1]) gt 0.4 then $ + yrng = [yran[0]<(-0.4),yran[1]>0.4] $ + else yrng = [yran[0]<(-0.2),yran[1]>0.2] + plot,fitws,diff,psym=4,charsi=si,charthi=th,thick=th, $ + title=imglab+' (bar '+strn(b)+')', $ + xthick=th, xtitle='Wavelength (A)', $ + ythick=th, ytitle='Atlas - Spec (A)',yrange=yrng,/ys + oplot,!x.crange,[0,0],linestyle=2,color=colordex('blue') + oplot,!x.crange,[sigmas[b],sigmas[b]],linesty=1, $ + color=colordex('blue') + oplot,!x.crange,-[sigmas[b],sigmas[b]],linesty=1, $ + color=colordex('blue') + kcwi_legend,["RMS = "+string(sigmas[b],form='(f7.3)')+' Ang',$ + 'NMatch = '+strn(ngood), $ + 'NRej = '+strn(nbad)],charsi=si,charthi=th + read,'Next? (Q - quit plotting, - next): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then display = (1 eq 0) + endif + endif ; b ne refbar +endfor +; +; our final coefficients +fincoeff = twkcoeff +; +; get wavelength ranges +y0wvs = fltarr(120) ; good wavelengths +y1wvs = fltarr(120) +t0wvs = fltarr(120) ; trim wavelengths +t1wvs = fltarr(120) +ym0wvs = fltarr(120) ; masked good wavelengths +ym1wvs = fltarr(120) +tm0wvs = fltarr(120) ; masked trim wavelengths +tm1wvs = fltarr(120) +for b=0,119 do begin + y0wvs[b] = poly(kgeom.goody0,fincoeff[*,b]) + y1wvs[b] = poly(kgeom.goody1,fincoeff[*,b]) + t0wvs[b] = poly(kgeom.trimy0,fincoeff[*,b]) + t1wvs[b] = poly(kgeom.trimy1,fincoeff[*,b]) + ym0wvs[b] = poly(kgeom.goodmy0,fincoeff[*,b]) + ym1wvs[b] = poly(kgeom.goodmy1,fincoeff[*,b]) + tm0wvs[b] = poly(kgeom.trimmy0,fincoeff[*,b]) + tm1wvs[b] = poly(kgeom.trimmy1,fincoeff[*,b]) +endfor +y0max = max(y0wvs) +y0min = min(y0wvs) +y1max = max(y1wvs) +y1min = min(y1wvs) +trim0 = min(t0wvs) +trim1 = max(t1wvs) +ym0max = max(ym0wvs) +ym0min = min(ym0wvs) +ym1max = max(ym1wvs) +ym1min = min(ym1wvs) +trimm0 = min(tm0wvs) +trimm1 = max(tm1wvs) +; +; check for negative dispersion +if trim0 gt trim1 then begin + trim0 = min(t1wvs) + trim1 = max(t0wvs) + trimm0 = min(tm1wvs) + trimm1 = max(tm0wvs) +endif +; +wavgood0 = min([y0max,y1max]) +wavgood1 = max([y0min,y1min]) +wavall0 = min([y0min,y1min]) +wavall1 = max([y0max,y1max]) +wavgoodm0 = min([ym0max,ym1max]) +wavgoodm1 = max([ym0min,ym1min]) +wavallm0 = min([ym0min,ym1min]) +wavallm1 = max([ym0max,ym1max]) +; +; mid-wavelength +wavmid = total([wavgood0,wavgood1,wavall0,wavall1])/4. +; +; now let's update kgeom struct +; +; check bar statuses +kgeom.status = max(barstat) +if kgeom.status gt 0 then begin + kcwi_print_info,ppar,pre,'Some bars had issues, check log for details',/warning + ; + ; if N&S, then this was just for info, carry on + if kgeom.nasmask then $ + kgeom.status=0 +endif +; +; order of the fits +kgeom.bfitord = degree+1 +; +; reference bar coeff's +kgeom.rbcoeffs = fincoeff[*,refbar] +; +; inclusive, all, and trim wavelengths +kgeom.waveall0 = wavall0 +kgeom.waveall1 = wavall1 +kgeom.wavegood0 = wavgood0 +kgeom.wavegood1 = wavgood1 +kgeom.waveallm0 = wavallm0 +kgeom.waveallm1 = wavallm1 +kgeom.wavegoodm0 = wavgoodm0 +kgeom.wavegoodm1 = wavgoodm1 +kgeom.wavemid = wavmid +; +; make sure wavelengths are on fiducial scale and +; zeropoint is integer number of bins from fiducial wl +ndels = long((trim0 - kgeom.wavefid)/kgeom.dwout) +kgeom.wave0out = kgeom.wavefid + float(ndels) * kgeom.dwout +ndels = long((trim1 - kgeom.wavefid)/kgeom.dwout) +kgeom.wave1out = kgeom.wavefid + float(ndels) * kgeom.dwout +; +; log values +kcwi_print_info,ppar,pre,'min,max good wavelengths',kgeom.wavegood0,kgeom.wavegood1, $ + format='(a,2f9.2)' +kcwi_print_info,ppar,pre,'min,max inclusive wavelengths',kgeom.waveall0,kgeom.waveall1, $ + format='(a,2f9.2)' +kcwi_print_info,ppar,pre,'min,max trim wavelengths',kgeom.wave0out,kgeom.wave1out, $ + format='(a,2f9.2)' +kcwi_print_info,ppar,pre,'middle wavelength',kgeom.wavemid,format='(a,f9.2)' +; +; set fit rms values +mom = moment(sigmas) +kgeom.avewavesig = mom[0] +kgeom.stdevwavesig = sqrt(mom[1]) +; +; log results +kcwi_print_info,ppar,pre,'ave/stddev bar wavelength sigma (Ang)', $ + mom[0],sqrt(mom[1]),format='(a,f7.3,2x,f7.3)' +kcwi_print_info,ppar,pre,'min/max bar wavelength sigma (Ang)', $ + minmax(sigmas), format='(a,f7.3,2x,f7.3)' +outliers = where(abs(sigmas-mom[0]) gt 3.*sqrt(mom[1]), noutliers) +if noutliers gt 0 then $ + kcwi_print_info,ppar,pre,'> 3sig outlier bars present: Imgnum, Bars', $ + imgnum,outliers,format='(a,i7,2x,'+strn(noutliers)+'i5)', $ + /warning +; +; central dispersion +disp = dblarr(120) +; +; data coefficients +dcoeff = dblarr(9) +; +; now apply all the coeffs and calculate central dispersion +for b=0,119 do begin + ; + ; update coeffs in kgeom struct + dcoeff = fincoeff[*,b] + kcwi_apply_coeffs,kgeom,b,dcoeff + ; + ; calculate central dispersion + subxvals = xvals[minrow:maxrow] + subwvals = poly(subxvals,fincoeff[*,b]) + ; + for d=0,degree do dcoeff[d] = d*fincoeff[d,b] + dcoeff = shift(dcoeff,-1) + xv = interpol(subxvals,subwvals,cwvl) + disp[b] = poly(xv,dcoeff) +endfor ; b +; +; stamp the KCWI_GEOM struct +kgeom.progid = pre +kgeom.timestamp = systime(1) +; +; now let's make some plots! +; +; plot diagnostics of arc fits +kcwi_plot_arcfits,specs,kgeom,ppar,cntcoeff,fincoeff,sigmas,fwaves,dwaves, $ + /tweak,plot_file=plot_file,ftype=ftype +; +; finally, let's make a diagnostic summary plot +if ppar.display ge 1 then begin + deepcolor + !p.background=colordex('white') + !p.color=colordex('black') + tlab = imglab+', '+ftype + ; + ; get fit rms status + mom = moment(sigmas) + fitrmslab = strtrim(string(mom[0],format='(f9.3)'),2) + ' +- ' + $ + strtrim(string(sqrt(mom[1]),format='(f9.3)'),2) + yrngs = get_plotlims(sigmas) + !p.multi=[0,1,2] + plot,sigmas,psym=4,charsi=si,charthi=th,thick=th,title=tlab+' Fit : '+fitrmslab, $ + xthick=th,xtitle='Bar #',xrange=[-1,120],/xs, $ + ythick=th,ytitle='RMS (Ang)',yrange=yrngs,/ys + oplot,!x.crange,[mom[0],mom[0]],linesty=5,thick=th + oplot,!x.crange,[mom[0]+sqrt(mom[1]),mom[0]+sqrt(mom[1])],linesty=1,thick=th + oplot,!x.crange,[mom[0]-sqrt(mom[1]),mom[0]-sqrt(mom[1])],linesty=1,thick=th + kcwi_oplot_slices + plot,disp,psym=4,charsi=si,charthi=th,thick=th, $ + title=tlab+' Dispersion @ '+string(cwvl,"(f8.2)")+' Ang', $ + xthick=th,xtitle='Bar #',xrange=[-1,120],/xs, $ + ythick=th,ytitle='Ang/pix',yrange=yrngd,/ys + kcwi_oplot_slices + if ppar.display ge 2 then read,'next: ',q +endif +!p.multi=0 +; +return +end ; kcwi_solve_arcs diff --git a/pro/kcwi_solve_arcs_iter.pro b/pro/kcwi_solve_arcs_iter.pro new file mode 100644 index 0000000..70a24c5 --- /dev/null +++ b/pro/kcwi_solve_arcs_iter.pro @@ -0,0 +1,887 @@ +; +; Copyright (c) 2013-2016, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_SOLVE_ARCS_ITER +; +; PURPOSE: +; Solve the wavelength solutions for each bar of the arc spectrum +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_SOLVE_ARCS_ITER, Specs, Kgeom, Ppar +; +; INPUTS: +; Specs - a array of arc spectra produced by KCWI_EXTRACT_ARCS +; Cntcoeff - Central fit coefficients (from kcwi_fit_center.pro) +; Kgeom - KCWI_GEOM struct from KCWI_TRACE_CBARS and KCWI_EXTRACT_ARCS +; Ppar - KCWI_PPAR pipeline parameter struct +; +; INPUT KEYWORDS: +; TWEAK - set to iteratively adjust wavelength solution to improve +; fit in outer thirds of wavelength range: for +; nod-and-shuffle images only one iteration is performed +; and only to assess rms of fit (no actual tweak is done). +; PLOT_FILE - if set to a string, will be used to produce +; postscript output of diagnostic plots +; +; SIDE EFFECTS: +; Modifies KCWI_GEOM struct by updating the bar coeffiecients and +; wavelength fit order and calculating new control points that take +; into account the wavelength solution. +; NOTE: sets KGEOM.STATUS to 0 if fitting succeeded, otherwise sets to +; 1 or greater depending on reason for failure (see code below). +; +; PROCEDURE: +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Matt Matuszewski +; 2013-DEC-10 Initial Revision +; 2013-DEC-12 Changes to make the cross-correlation more robust +; 2014-MAR-19 Use KCWI_CLEAN_COEFFS to fix errant bars +; 2014-MAR-20 Adjusts peak match thresh if no peaks matched, regardless of iteration number +; 2014-MAR-28 Use central ccor peak for initial wave soln. +; 2014-MAR-28 Check for scattered light in initial offset calc. +; 2014-APR-07 Changed to a function returning status of fit +; 2014-APR-07 Uses preliminary solution to determine wave range instead of filter +; 2014-APR-07 Uses parabola fit to CC peak when no zero derivs are found +; 2014-SEP-11 Converted function to pro with status in kgeom.status +; 2014-SEP-16 Single iteration mode for assessing rms of fits (N&S) +; 2014-NOV-06 Put stats (rms) after all tweaking done +; 2015-APR-22 Rework of peak finding +;- +; +pro kcwi_solve_arcs_iter, specs, cntcoeff, kgeom, ppar, $ + tweak=tweak, plot_file=plot_file + +pre = 'KCWI_SOLVE_ARCS_ITER' +q='' +; +; check status of central fit +if kgeom.status gt 0 then begin + kcwi_print_info,ppar,pre,'Cannot solve arcs',/error + return +endif +; +; do we want to display stuff? +display = (ppar.display ge 2) +; +if ppar.display ge 1 then begin + deepcolor + !p.multi=0 + !p.background=colordex('white') + !p.color=colordex('black') + th=2.0 + si=1.5 +endif +; +; peak matching +pkdel = ppar.pkdel +; +; which image number +imgnum = kgeom.arcimgnum +; +; which slicer? +ifunum = kgeom.ifunum +ifunam = kgeom.ifunam +; +; which grating? +grating = kgeom.gratid +; +; which filter? +filter = kgeom.filter +; +; image label +imglab = 'Img # '+strn(imgnum)+' ('+kgeom.refname+') Sl: '+strtrim(ifunam,2)+ $ + ' Fl: '+strtrim(filter,2)+' Gr: '+strtrim(grating,2) +; +; is this N+S mask in? +nasmask = kgeom.nasmask +; +; central wavelength? +cwvl = kgeom.cwave +; +; canonical resolution in Angstroms? +resolution = kgeom.resolution * float(kgeom.ybinsize) +; +; log info +kcwi_print_info,ppar,pre,systime(0) +kcwi_print_info,ppar,pre,'img, grat, filt, nasmsk, cwave', $ + imgnum,grating,filter,nasmask,cwvl, $ + format='(a,i6,2x,a-8,2x,a-8,i4,f12.3)' +; +; we will be initially using a third degree fit to the peak positions: +degree = 3 +; +; last degree is for full ccd tweaked fits +lastdegree = kgeom.lastdegree +; +; which is the reference bar? +refbar = kgeom.refbar +; +; load the reference atlas spectrum. +kcwi_read_atlas,kgeom,ppar,refspec,refwvl,refdisp +; +; input spectrum dimensions +specsz = size(specs,/dim) +; +; set up array with zero point in the center +x0 = specsz[0]/2 +xvals = dindgen(specsz[0])-x0 +; +; keep track of final sigmas +sigmas = dblarr(120) +; +; keep track of individual bar fits +barstat = intarr(120) +barrej = intarr(120) +; +; keep track of observed versus atlas comparison +fwaves = fltarr(120,1000) +dwaves = fltarr(120,1000) +; +; keep track of reference wavelengths and matched pixel positions +rwaves = fltarr(120,1000) +xcents = fltarr(120,1000) +; +; for tweaking: array of border fractions to ignore. +if nasmask then begin + edgespace = [ 0.34 ] + minrow = (1*specsz[0])/3 + maxrow = (2*specsz[0])/3 + ftype = 'Central' +endif else begin + edgespace = [ 0.25, 0.18, 0.12, 0.06, 0.00] + lastrow = 50 + minrow = lastrow + maxrow = specsz[0]-lastrow-1 + ftype = 'Iter' + if strpos(grating,'BL') ge 0 then begin + edgespace = [0.00] + ftype = 'LowRes' + lastdegree = degree + pkdel = pkdel * 2.0 + endif +endelse ; no nasmask +; +; find wavelength range +mnwvs = fltarr(120) +mxwvs = fltarr(120) +for b=0,119 do begin + waves = poly(xvals,cntcoeff[*,b]) + mnwvs[b] = min(waves) + mxwvs[b] = max(waves) +endfor +; +minwav = min(mnwvs) +maxwav = max(mxwvs) +; +; adjust higher order terms, if requested +; do a continuum subtraction if tweaking +if keyword_set(tweak) then begin + if nasmask then sectors = 6 else sectors = 16 + ; + div = (maxrow-minrow)/sectors + ; + xv = fltarr(sectors) + yv = fltarr(sectors) + for b = 0, 119 do begin + ; + ; get continuum point in each sector + for sec = 0,sectors-1 do begin + ; + ; We may need this for low-res gratings! + mn = min(specs[minrow+sec*div:minrow+(sec+1)*div,b],mi) + xv[sec] = mi+minrow+sec*div + ;mn = median(specs[minrow+sec*div:minrow+(sec+1)*div,b]) + ;xv[sec] = minrow+sec*div + div/2 + yv[sec] = mn + endfor + ; + ; fit sector sky points + res = poly_fit(xv,yv,3,yfit=yf) + ; + xvs = findgen(specsz[0]) + specs[*,b] -= res[0]+res[1]*xvs+res[2]*xvs^2+res[3]*xvs^3 + ; + endfor ; b + +endif ; this is part of the tweak. +; +; --- this is where the tweak happens, if requested +; this is how we do it: +; currently, we are not going to be tweaking nod-and-shuffle +; We assume that the initial solution is pretty good and only needs to +; be refined toward the edges of the CCD. +; We step through successively larger windows starting with the inner +; most 1/3 of the CCD +if keyword_set(tweak) then begin + ; + ; get number of iterations + niter = n_elements(edgespace) + if niter eq 1 then $ + kcwi_print_info,ppar,pre,'using '+strn(niter)+' iteration to measure RMS of initial central fit' $ + else kcwi_print_info,ppar,pre,'using '+strn(niter)+' iterations to tweak higher order terms' + ; + kcwi_print_info,ppar,pre,'Spectral line finding/matching thresh (frac. of res.): PkDel', $ + pkdel,format='(a,2x,f7.3)' + ; + ; start with central fit coefficients + twkcoeff = cntcoeff + ; + ; reference spectrum and wavelengths + twk_reference_spectrum = refspec + twk_reference_wavelengths = refwvl + ; + ; first trim to the grating bandpass (and a buffer of 100 Ang). + qwvs = where(twk_reference_wavelengths gt minwav-100. and $ + twk_reference_wavelengths lt maxwav+100.,nqwvs) + if nqwvs eq 0 then begin + kcwi_print_info,ppar,pre,'No wavelength overlap with atlas spectrum',/error + kgeom.status=5 + return + endif + twk_reference_spectrum = twk_reference_spectrum[qwvs] + twk_reference_wavelengths = twk_reference_wavelengths[qwvs] + ; + ; plot if requested + if ppar.display ge 3 then begin + xarng = get_plotlims([minwav,maxwav],pad=0.2) + ; + ; if plotting diagnostics, just focus on first iteration (central third) + if ppar.display ge 3 then begin + dwav = (maxwav - minwav)/3. + xarng = get_plotlims([minwav+dwav,maxwav-dwav],pad=0.3) + atitle = 'Atlas ('+kgeom.refname+') Spectrum Lines (Central 3rd)' + endif else $ + atitle = 'Atlas ('+kgeom.refname+') Spectrum Lines' + plot,twk_reference_wavelengths,twk_reference_spectrum, title=atitle, $ + thick=th,charsi=si,charthi=th, $ + xthick=th,xtitle='Wave(A)',xrange=xarng,/xs, $ + ythick=th,ytitle='Flux', $ + yrange=[min(twk_reference_spectrum)>10., $ + max(twk_reference_spectrum)+0.05*max(twk_reference_spectrum)],/ys,/ylog + endif + ; + ; let's find the peaks in the reference spectrum. + smooth_width = fix(resolution/refdisp)>4 ; in pixels + slope_thresh = 0.003 + fthr = 0.05 ; start at 5% of max + ampl_thresh = max(twk_reference_spectrum)*fthr + twk_ref_cent = findpeaks(twk_reference_wavelengths,twk_reference_spectrum, $ + smooth_width,slope_thresh,ampl_thresh,count=twk_ref_npks) + newlines = twk_ref_npks + oldlines = newlines + ; + ; find where noise starts to contaminate line list + while (newlines - oldlines) lt 10 and fthr gt 0.01 do begin + fthr -= 0.01 + oldlines = newlines + ampl_thresh = max(twk_reference_spectrum)*fthr + twk_ref_cent = findpeaks(twk_reference_wavelengths,twk_reference_spectrum, $ + smooth_width,slope_thresh,ampl_thresh,count=twk_ref_npks) + newlines = twk_ref_npks + endwhile + fthr += 0.01 + ampl_thresh = max(twk_reference_spectrum)*fthr + twk_ref_cent = findpeaks(twk_reference_wavelengths,twk_reference_spectrum, $ + smooth_width,slope_thresh,ampl_thresh,count=twk_ref_npks) + ; + if twk_ref_npks le lastdegree then begin + kcwi_print_info,ppar,pre,'Not enough good atlas points found',twk_ref_npks,/error + kgeom.status=6 + return + endif + ; + if twk_ref_npks lt 50 then $ + kcwi_print_info,ppar,pre,'Atlas wavelength coverage may be limited: N lines < 50',/warn + + kcwi_print_info,ppar,pre,'Atlas threshhold in percent of max',fix(fthr*100.) + kcwi_print_info,ppar,pre,'Number of clean atlas lines found',twk_ref_npks,$ + format='(a,i4)' + kcwi_print_info,ppar,pre,'Matching width in Angstroms',pkdel*refdisp,$ + format='(a,f6.4)' + ; + if ppar.display ge 3 then begin + for j=0,twk_ref_npks-1 do begin + oplot,[twk_ref_cent[j],twk_ref_cent[j]],10.^!y.crange,color=colordex('blue'),thick=1.0 + endfor + oplot,[minwav,minwav],10.^!y.crange,color=colordex('green'),thick=th,linesty=5 + oplot,[maxwav,maxwav],10.^!y.crange,color=colordex('green'),thick=th,linesty=5 + kcwi_legend,['PeakFit','WavRng'],linesty=[0,5],thick=[1.0,th], $ + color=[colordex('blue'),colordex('green')],/clear,charsi=si,charthi=th + kcwi_legend,['Thr='+string(fix(fthr*100.),form='(i02)')+'%'], $ + /clear,/right, charthi=th,charsi=si + endif + ; + ; at this point we have the peaks we want in the atlas spectrum. + ; + ; write out atlas line list + atllist = ppar.reddir+ppar.froot+string(imgnum,'(i0'+strn(ppar.fdigits)+')')+'_atlas.txt' + openw,al,atllist,/get_lun + printf,al,'# '+pre+': Atlas lines (Ang) printed on '+systime(0) + for j=0,twk_ref_npks-1 do printf,al,twk_ref_cent[j],format='(f12.3)' + free_lun,al + kcwi_print_info,ppar,pre,'Atlas lines written to',atllist,format='(a,a)' + ; + ; open file for output of object line list + objlist = ppar.reddir+ppar.froot+string(imgnum,'(i0'+strn(ppar.fdigits)+')')+'_object.txt' + openw,al,objlist,/get_lun + printf,al,'# '+pre+': Object lines (Ang) printed on '+systime(0) + printf,al,'# AtlWl Ypx Bar' + ; + ; now pop up a diagnostic window if requested + ddisplay = (ppar.display ge 3) + if ddisplay then $ + window,1,title='kcwi_solve_arcs_iter(2)' + ; + ; now loop over the iterations and over the bars + ; + if nasmask and niter gt 1 then $ + kcwi_print_info,ppar,pre,'tweaking not recommended when the NS mask is in',/warning + ; + centerpoints = intarr(niter) + leftpoints = centerpoints + rightpoints = centerpoints + ; + for iter = 0,niter-1 do begin + ; if this is the last iteration then we should see the full CCD + ; at this point, we want switch to last degree. + if iter eq niter-1 and not nasmask then degree = lastdegree + ; what region of the ccd are we fitting in this iteration? + twk_minrow = fix(specsz[0]*edgespace[iter]) + twk_maxrow = fix(specsz[0]*(1-edgespace[iter]))1., title='Object Spectrum Lines (Central 3rd): Bar '+strn(b), $ + thick=th,charsi=si,charthi=th, $ + xthick=th,xtitle='Wave(A)',xrange=xarng,/xs, $ + ythick=th,ytitle='Flux', $ + yrange=[min(subyvals)>1.,max(subyvals)+0.05*max(subyvals)],/ys,/ylog + endif + ; + ; find good peaks in object spectrum + smooth_width = fix(resolution/abs(twkcoeff[1,b]))>4 ; in pixels + peak_width = fix(smooth_width*1.5) ; for fitting peaks + slope_thresh = 0.7*smooth_width/2./100.0 ; more severe for object + fthr = 0.15 ; start at 15% of max + ampl_thresh = max(subyvals) * fthr + twk_spec_cent = findpeaks(subwvals,subyvals,smooth_width,slope_thresh, $ + ampl_thresh,peak_width,count=twk_spec_npks) + newlines = twk_spec_npks + oldlines = newlines + ; + ; find where the noise starts to contaminate the line list + while (newlines - oldlines) lt 10 and fthr gt 0.01 do begin + oldlines = newlines + fthr -= 0.01 + ampl_thresh = max(subyvals) * fthr + twk_spec_cent = findpeaks(subwvals,subyvals,smooth_width,slope_thresh, $ + ampl_thresh,peak_width,count=twk_spec_npks) + newlines = twk_spec_npks + endwhile + ; + ; put threshhold back above the noise + fthr += 0.02 + ampl_thresh = max(subyvals) * fthr + twk_spec_cent = findpeaks(subwvals,subyvals,smooth_width,slope_thresh, $ + ampl_thresh,peak_width,count=twk_spec_npks) + ; + ; test final line list + if twk_spec_npks gt 0 then begin + ; + ; at this point we have the catalog of good reference points (from + ; before) and list of good points in this bar. We have to associate + ; them. + one_ref = fltarr(twk_ref_npks)+1.0 + one_spec = fltarr(twk_spec_npks)+1.0 + + diff = abs((twk_spec_cent##one_ref) - (one_spec##twk_ref_cent)) + ; + mn = min(diff,dim=1,mi) + ; + ; here we match the peaks to one another. + pkd = pkdel > 1.0 + ; + ; never let this get smaller than a single atlas pixel + pkm = pkd*refdisp ; match thresh in Angstroms + matchedpeaks = where(mn lt pkm, nmatchedpeaks) + ; + ; for first iteration, make sure we have enough peaks + ; also handle if we have no matched peaks + if iter eq 0 or nmatchedpeaks eq 0 then begin + orig_nmp = nmatchedpeaks + while nmatchedpeaks lt 5 and pkd lt 2. do begin + ; + ; open up the match criterion + pkd += 0.5 + ; + ; try again + pkm = pkd*refdisp + matchedpeaks = where(mn lt pkm, nmatchedpeaks) + endwhile + ; + ; report any adjustments + if pkd ne pkdel then begin + print,'' + print,'Bar: ',b,', pkdel updated to ',pkd, '; ',orig_nmp, $ + ' --> ',nmatchedpeaks,' peaks', $ + format='(a,i3,a,f5.2,a,i2,a,i3,a)' + endif + endif + ; + ; no matched peaks + endif else begin + nmatchedpeaks = 0 + kcwi_print_info,ppar,pre,'No good peaks in arc for Bar # '+strn(b),/warning + endelse + ; + if nmatchedpeaks le 0 then begin + kcwi_print_info,ppar,pre,'No peaks matched in Bar # '+strn(b),/warning + barstat[b]=9 + goto, errbar + endif + ; + ; plot matches, if requested + if ddisplay and iter eq 0 then begin + for pp = 0,n_elements(twk_spec_cent)-1 do $ + oplot,[twk_spec_cent[pp],twk_spec_cent[pp]], $ + 10.^!y.crange,color=colordex('R') + for pp = 0,nmatchedpeaks-1 do $ + oplot,[twk_spec_cent[matchedpeaks[pp]],twk_spec_cent[matchedpeaks[pp]]], $ + 10.^!y.crange,color=colordex('G') + endif + ; + if nmatchedpeaks le degree then begin + print,' ' + kcwi_print_info,ppar,pre,'Not enough peaks matched in Bar # '+strn(b)+': '+strn(nmatchedpeaks),/warning + barstat[b]=9 + goto, errbar + endif + ; + twk_ref_idx = mi[matchedpeaks] mod twk_ref_npks + twk_spec_idx = matchedpeaks + ; + targetw = twk_ref_cent[twk_ref_idx] + ; + ; check if subwvals is monotonic (cspline requires this) + testw = subwvals(sort(subwvals)) + diff = subwvals - testw + junk = where(diff ne 0., ndiff) + if ndiff gt 0 then begin + kcwi_print_info,ppar,pre,'wavelengths not monotonic',/warn + initx = twk_spec_cent[twk_spec_idx] + endif else initx = cspline(subwvals,subxvals,twk_spec_cent[twk_spec_idx]) + ; + ; we need reverse coefficients. + fitdegree = degree < (nmatchedpeaks-1) + newcoeff = poly_fit(initx,targetw,fitdegree,yfit=fitw) + ; + ; get residual stats + ims,(targetw - fitw),rmn,rsg,rwgt + bad = where(rwgt le 0, nbad) + ; + ; reject baddies + nrej = 0 + while nbad gt 0 do begin + ; + ; count the bad guys + nrej += nbad + ; + ; overplot them if we are displaying + if ddisplay and iter eq 0 then begin + for pp = 0,nbad-1 do $ + oplot,[targetw[bad[pp]],targetw[bad[pp]]], $ + 10.^!y.crange,color=colordex('B') + endif + ; + ; get the good guys + good = where(rwgt eq 1, nmatchedpeaks) + ; + ; check to be sure we have some left + if nmatchedpeaks le 0 then begin + kcwi_print_info,ppar,pre,'All points rejected',/error + barstat[b]=10 + goto, errbar + endif + ; + ; clean fit data + initx = initx[good] + targetw = targetw[good] + ; + ; re-fit + fitdegree = degree < (nmatchedpeaks-1) + newcoeff = poly_fit(initx,targetw,fitdegree,yfit=fitw) + ; + ; check for more baddies + ims,(targetw - fitw),rmn,rsg,rwgt + bad = where(rwgt le 0, nbad) + endwhile + ; + ; record results + print,string(13B),iter+1,niter,b,nmatchedpeaks,nrej, $ + format='($,a1,"Iteration: ",i2," of ",i2," Bar:",i3," Valid peaks:",i3," Rejected peaks:",i3)' + ; + ; plot legend and query user + if ddisplay and iter eq 0 then begin + kcwi_legend,['Good','NoAtlas','Reject'],linesty=[0,0,0], $ + /bottom,/clear, charthi=th,charsi=si, $ + color=[colordex('G'),colordex('R'),colordex('B')] + kcwi_legend,['Thr='+string(fix(fthr*100.),form='(i02)')+'%'], $ + /clear,/right, charthi=th,charsi=si + print,'' + read,'Next? (Q - quit plotting, - next): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then ddisplay = (1 eq 0) + endif + ; + ; store for later + rwaves[b,0:(nmatchedpeaks-1)] = targetw + xcents[b,0:(nmatchedpeaks-1)] = initx + ; + ; write out + if iter eq niter-1 then begin + for j=0,nmatchedpeaks-1 do $ + printf,al,targetw[j],initx[j],b,format='(f12.3,f15.3,i7)' + endif + ; + if nmatchedpeaks ge degree+2 then begin + twkcoeff[*,b] = 0 + twkcoeff[0:fitdegree,b] = newcoeff + endif + ; + ; if we get here, bar fit was good + barstat[b] = 0 + barrej[b] = nrej + ; + ; target for bars with problems + errbar: + endfor ; b + ; + ; report iteration + print,string(13B),format='($,a1)' + kcwi_print_info,ppar,pre,"Iteration "+string(iter+1,"(i2)")+" of "+ $ + string(niter,"(i2)")+" done. " + ; + ; now clean each slice of outlying bars + ; don't bother for a single iteration + if kgeom.bclean and niter gt 1 then begin + ; + ; only go interactive on last iteration + if iter lt niter-1 then $ + clnplot = 0 $ + else clnplot = 1 + if ppar.display ge 3 and clnplot then $ + window,2,title='kcwi_clean_coeffs' + kcwi_clean_coeffs,twkcoeff,degree,ppar,plot=clnplot + if ppar.display ge 3 and clnplot then $ + wdelete,2 + endif + endfor; iter + ; + ; done tweaking + print,'' + print,'Done tweaking: now get final stats.' +endif; tweak +; endtweak +; +; close object line output list +free_lun,al +kcwi_print_info,ppar,pre,'Object lines written to',objlist,format='(a,a)' +; +; use the tweaked coefficients, if asked to. +if keyword_set(tweak) and niter gt 1 then $ + usecoeff = twkcoeff $ +else usecoeff = cntcoeff +; +; if we have tweaked at all let's get the final status +if keyword_set(tweak) then begin + ; + ; loop over the bars + for b=0,119 do begin + ; + ; get output coeffs + scoeff = pascal_shift(reform(usecoeff[*,b]),x0,/silent) + ; + ; get ref waves and corresponding object x values + targetw = rwaves[b,*] + initx = xcents[b,*] + ; + ; matched peaks are where we have a good ref wave + matchedpeaks = where(targetw gt 0., nmatchedpeaks) + if nmatchedpeaks gt 2 then begin + ; + ; get the good ones + targetw = targetw[matchedpeaks] + initx = initx[matchedpeaks] + ; + ; ref - observed wavelengths + finwaves = poly(initx,usecoeff[*,b]) + cmpdiff = targetw-finwaves + ; + ; account for cleaning zero-point wavelength + mom = moment(cmpdiff) + cmpdiff = cmpdiff - mom[0] + ; + ; final RMS + mom = moment(cmpdiff) + sigmas[b] = sqrt(mom[1]) + ; + ; store for later plots + fwaves[b,0:(nmatchedpeaks-1)] = finwaves + dwaves[b,0:(nmatchedpeaks-1)] = cmpdiff + ; + ; log stats + kcwi_print_info,ppar,pre,ftype+' Fit: Bar#,Npks,RMS,Cdisp,Coefs', $ + b,nmatchedpeaks,sigmas[b],usecoeff[1,b],scoeff[0:degree], $ + format='(a,2i5,2x,f7.3,f8.4,f9.2,f8.4,'+strn((degree-1)>1)+'g13.5)' + ; + ; for your viewing pleasure + if display then begin + ; + ; back to main plot window + wset,0 + if nasmask then begin + dwav = (maxwav - minwav)/3. + xrng = get_plotlims([minwav+dwav,maxwav-dwav],pad=0.3) + endif else xrng = [minwav,maxwav] + yran = get_plotlims(cmpdiff) + if abs(usecoeff[1,refbar]) gt 0.4 then $ + yrng = [yran[0]<(-0.4),yran[1]>0.4] $ + else yrng = [yran[0]<(-0.2),yran[1]>0.2] + plot,finwaves,cmpdiff,psym=4,charsi=si,charthi=th,thick=th, $ + xthick=th,xtitle='Wave(Ang)',xrange=xrng,/xs, $ + ythick=th,ytitle='Resid(Ang)',yrange=yrng,/ys, $ + title=imglab+', Bar: '+string(b,"(i3)") + $ + ', Slice: '+string(fix(b/5),"(i2)") + oplot,!x.crange,[0,0],linestyle=2,color=colordex('blue') + oplot,!x.crange,[sigmas[b],sigmas[b]], $ + linesty=1,color=colordex('blue') + oplot,!x.crange,-[sigmas[b],sigmas[b]], $ + linesty=1,color=colordex('blue') + kcwi_legend,["RMS = "+string(sigmas[b],form='(f7.3)')+' Ang', $ + 'NMatch = '+strn(nmatchedpeaks), $ + 'NRej = '+strn(barrej[b])], $ + charsi=si,charthi=th + read,'Next? (Q - quit plotting, - next): ',q + if strupcase(strmid(q,0,1)) eq 'Q' then display = (1 eq 0) + endif ; display + endif else $ ; nmatchedpeaks gt 2 + kcwi_print_info,ppar,pre,'Not enough matched peaks for bar',b,nmatchedpeaks,/error + endfor ; loop over bars to get final stats +endif ; have we done any tweaking? +; +; clean up +if ppar.display ge 3 then wdelete,1 +; +; our final coefficients +fincoeff = dblarr(9,120) +; +; shift solution reference to pixel 0 instead of center +for b=0,119 do $ + fincoeff[*,b] = pascal_shift(usecoeff[*,b],x0,/silent) +; +; good bars +bargood = where(barstat eq 0, nbargood) +if nbargood le 0 then begin + kcwi_print_info,'No good bar fits',/error + kgeom.status = 10 + return +endif +; +; bad bars +barbad = where(barstat gt 0, nbarbad) +if nbarbad gt 0 then $ + sigmas[barbad] = -1. +; +; get wavelength ranges +y0wvs = fltarr(120) ; good wavelengths +y1wvs = fltarr(120) +t0wvs = fltarr(120) ; trim wavelengths +t1wvs = fltarr(120) +for b=0,119 do begin + wy0 = poly(kgeom.goody0,fincoeff[*,b]) + y0wvs[b] = wy0 + wy1 = poly(kgeom.goody1,fincoeff[*,b]) + y1wvs[b] = wy1 + wt0 = poly(kgeom.trimy0,fincoeff[*,b]) + t0wvs[b] = wt0 + wt1 = poly(kgeom.trimy1,fincoeff[*,b]) + t1wvs[b] = wt1 +endfor +; +; for central fit only, all bars must be used +if nasmask then begin + y0max = max(y0wvs) + y0min = min(y0wvs) + y1max = max(y1wvs) + y1min = min(y1wvs) + trim0 = min(t0wvs) + trim1 = max(t1wvs) +endif else begin + y0max = max(y0wvs[bargood]) + y0min = min(y0wvs[bargood]) + y1max = max(y1wvs[bargood]) + y1min = min(y1wvs[bargood]) + trim0 = min(t0wvs[bargood]) + trim1 = max(t1wvs[bargood]) +endelse +; +; check for negative dispersion +if trim0 gt trim1 then begin + trim0 = min(t1wvs[bargood]) + trim1 = max(t0wvs[bargood]) +endif +; +wavgood0 = min([y0max,y1max]) +wavgood1 = max([y0min,y1min]) +wavall0 = min([y0min,y1min]) +wavall1 = max([y0max,y1max]) +; +; mid-wavelength +wavmid = total([wavgood0,wavgood1,wavall0,wavall1])/4. +; +; now let's update kgeom struct +; +; check bar statuses +kgeom.status = max(barstat) +if kgeom.status gt 0 then begin + kcwi_print_info,ppar,pre,'Some bars had issues, check log for details',/warning + ; + ; if N&S, then this was just for info, carry on + if kgeom.nasmask then $ + kgeom.status=0 +endif +; +; order of the fits +kgeom.bfitord = degree+1 +; +; reference bar coeff's +kgeom.rbcoeffs = fincoeff[*,refbar] +; +; inclusive, all, and trim wavelengths +kgeom.waveall0 = wavall0 +kgeom.waveall1 = wavall1 +kgeom.wavegood0 = wavgood0 +kgeom.wavegood1 = wavgood1 +kgeom.wavemid = wavmid +; +; make sure wavelengths are on fiducial scale and +; zeropoint is integer number of bins from fiducial wl +ndels = long((trim0 - kgeom.wavefid)/kgeom.dwout) +kgeom.wave0out = kgeom.wavefid + float(ndels) * kgeom.dwout +ndels = long((trim1 - kgeom.wavefid)/kgeom.dwout) +kgeom.wave1out = kgeom.wavefid + float(ndels) * kgeom.dwout +; +; log values +kcwi_print_info,ppar,pre,'min,max good wavelengths',kgeom.wavegood0,kgeom.wavegood1, $ + format='(a,2f9.2)' +kcwi_print_info,ppar,pre,'min,max inclusive wavelengths',kgeom.waveall0,kgeom.waveall1, $ + format='(a,2f9.2)' +kcwi_print_info,ppar,pre,'min,max trim wavelengths',kgeom.wave0out,kgeom.wave1out, $ + format='(a,2f9.2)' +kcwi_print_info,ppar,pre,'middle wavelength',kgeom.wavemid,format='(a,f9.2)' +; +; set fit rms values +if keyword_set(tweak) then begin + mom = moment(sigmas[bargood]) + kgeom.avewavesig = mom[0] + kgeom.stdevwavesig = sqrt(mom[1]) + ; + ; log results + outliers = where(abs(sigmas[bargood]-mom[0]) gt 3.*sqrt(mom[1]), noutliers) + kcwi_print_info,ppar,pre,'average bar wavelength sigma (Ang)',mom[0],' +- ',sqrt(mom[1]),$ + format='(a,f7.3,a,f7.3)' + kcwi_print_info,ppar,pre,'min/max bar wavelength sigma (Ang)',minmax(sigmas[bargood]), $ + format='(a,f7.3,2x,f7.3)' + if noutliers gt 0 then $ + kcwi_print_info,ppar,pre,'> 3sig outlier bars present, may want to tweak ppar.pkdel: Imgnum, Bars', $ + imgnum,outliers,format='(a,i7,2x,'+strn(noutliers)+'i5)',/warning +endif +; +; central dispersion +disp = dblarr(120) +; +; data coefficients +dcoeff = dblarr(9) +; +; now apply all the coeffs and calculate central dispersion +for b=0,119 do begin + ; + ; update coeffs in kgeom struct + dcoeff = fincoeff[*,b] + kcwi_apply_coeffs,kgeom,b,dcoeff + ; + ; calculate central dispersion + subxvals = xvals[minrow:maxrow] + subwvals = poly(subxvals,fincoeff[*,b]) + ; + for d=0,degree do dcoeff[d] = d*fincoeff[d,b] + dcoeff = shift(dcoeff,-1) + xv = interpol(subxvals,subwvals,cwvl) + disp[b] = poly(xv,dcoeff) +endfor ; b +; +; stamp the KCWI_GEOM struct +kgeom.progid = pre +kgeom.timestamp = systime(1) +; +; now let's make some plots! +; +; plot diagnostics of arc fits +kcwi_plot_arcfits,specs,kgeom,ppar,cntcoeff,fincoeff,sigmas,fwaves,dwaves, $ + tweak=tweak,plot_file=plot_file,ftype=ftype +; +; finally, let's make a diagnostic summary plot +if ppar.display ge 1 then begin + deepcolor + !p.background=colordex('white') + !p.color=colordex('black') + tlab = imglab+', '+ftype + if keyword_set(tweak) then begin + ; + ; get fit rms status + mom = moment(sigmas[bargood]) + fitrmslab = strtrim(string(mom[0],format='(f9.3)'),2) + ' +- ' + $ + strtrim(string(sqrt(mom[1]),format='(f9.3)'),2) + yrngs = get_plotlims(sigmas[bargood]) + !p.multi=[0,1,2] + plot,sigmas,psym=4,charsi=si,charthi=th,thick=th,title=tlab+' Fit : '+fitrmslab, $ + xthick=th,xtitle='Bar #',xrange=[-1,120],/xs, $ + ythick=th,ytitle='RMS (Ang)',yrange=yrngs,/ys + oplot,!x.crange,[mom[0],mom[0]],linesty=5,thick=th + oplot,!x.crange,[mom[0]+sqrt(mom[1]),mom[0]+sqrt(mom[1])],linesty=1,thick=th + oplot,!x.crange,[mom[0]-sqrt(mom[1]),mom[0]-sqrt(mom[1])],linesty=1,thick=th + if nbarbad gt 0 then $ + oplot,barbad,replicate(mom[0],nbarbad),psym=7,thick=th + kcwi_oplot_slices + endif else begin + yrngd = get_plotlims(disp) + !p.multi=0 + endelse + plot,disp,psym=4,charsi=si,charthi=th,thick=th, $ + title=tlab+' Dispersion @ '+string(cwvl,"(f8.2)")+' Ang', $ + xthick=th,xtitle='Bar #',xrange=[-1,120],/xs, $ + ythick=th,ytitle='Ang/pix',yrange=yrngd,/ys + if nbarbad gt 0 then $ + oplot,barbad,disp[barbad],psym=7 + kcwi_oplot_slices + if ppar.display ge 2 then read,'next: ',q +endif +!p.multi=0 +; +return +end ; kcwi_solve_arcs_iter diff --git a/pro/kcwi_solve_dgeom.pro b/pro/kcwi_solve_dgeom.pro new file mode 100644 index 0000000..903f83e --- /dev/null +++ b/pro/kcwi_solve_dgeom.pro @@ -0,0 +1,285 @@ +; +; Copyright (c) 2016, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_SOLVE_DGEOM +; +; PURPOSE: +; Solve the direct mode image geometry using and arc and cbars image set +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_SOLVE_DGEOM, Kdgeom, Ppar +; +; INPUTS: +; Kdgeom - KCWI_DGEOM struct from KCWI_TRACE_CBARS and KCWI_EXTRACT_ARCS +; Ppar - KCWI_PPAR pipeline parameter struct +; +; INPUT KEYWORDS: +; +; SIDE EFFECTS: +; Modifies KCWI_DGEOM struct with results of geometry solution. +; NOTE: sets KDGEOM.STATUS to 0 if fitting succeeded, otherwise sets to +; 1 or greater depending on reason for failure. +; +; PROCEDURE: +; Find the arc flat images and fit their shape, then extract the +; arcbars and find the relative offsets +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2016-NOV-09 Initial Revision +; 2017-APR-06 Updates for real data, more diagnostic output +;- +; +pro kcwi_solve_dgeom,kdgeom,ppar, help=help +; +; startup +pre = 'KCWI_SOLVE_DGEOM' +q = '' +; +; check inputs +if n_params(0) lt 2 or keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', ArcSpec, Kdgeom' + return +endif +; +; Check structs +if kcwi_verify_geom(kdgeom,/init) ne 0 then return +if kcwi_verify_ppar(ppar,/init) ne 0 then return +; +; read in the arc image +arf = kdgeom.arcfname +arc = mrdfits(arf,0,ahdr,/silent) +sz = size(arc,/dim) +nx = sz[0] +ny = sz[1] +; +; read in bars image +cbf = kdgeom.cbarsfname +bar = mrdfits(cbf,0,bhdr,/silent) +sz = size(bar,/dim) +if sz[0] ne nx or sz[1] ne ny then begin + kcwi_print_info,ppar,pre,'Arc and Bars size mis-match, returning',/error + kdgeom.status=1 + return +endif +; +; useful params +xbin = kdgeom.xbinsize +ybin = kdgeom.ybinsize +; +; set up data arrays +data = fltarr(3,24,1000) - 1. +errs = fltarr(24,1000) - 1. +wids = fltarr(24,1000) - 1. +sl = -1 +last = 0 +maximx = 0 +maximy = 0 +minimy = ny +dx = 2 / xbin +bspec = fltarr(1000,24) +; +; loop over image columns +for i=4,nx-4 do begin + ; + ; extract column + col = reform(median(arc[(i-dx):(i+dx),*],dimen=1)) + ; + ; find start of slice + t = where(col gt kdgeom.rdnoise * 10., nt) + if nt gt kdgeom.minpix then begin + ; + ; are we starting a new slice? + if last le kdgeom.minpix then begin + sl += 1 + nxp = 0 + print,'Slice: ',sl + if sl gt 23 then begin + kcwi_print_info,ppar,pre,'Slice overflow',sl, $ + /warn + break + endif + endif + ; + ; get ranges for bright pixels + rangepar,t,tstr + ran = strsplit(tstr,',',/extract,count=n) + ; + ; do we have more than one? + if n gt 1 then begin + for j=0,n-1 do begin + rangepar,ran[j],ys + ; + ; if large, this is the one we want + if n_elements(ys) gt 5 then begin + ran = ran[j] + break + endif + endfor + endif + ; + ; get trace of arc image + rangepar,ran[0],iy + cnt = cntrd1d(iy,col[iy]) + err = fltarr(n_elements(iy)) + kdgeom.rdnoise + if kdgeom.do_gauss then begin + while n_elements(iy) lt 9 do begin + iy = [min(iy)-1, iy, max(iy)+1] + err = [kdgeom.rdnoise, err, kdgeom.rdnoise] + endwhile + startp = [max(col[iy]), cnt, 2.] + yfit = gaussfit(iy, col[iy], pars, nterm=3, $ + estimates=startp, measure_err=err, sigma=perr) + if n_elements(perr) eq 3 then begin + fcnt = pars[1] + fwid = pars[2] + endif else begin + perr = [1.e9, 1.e9] + fcnt = 1.e9 + fwid = 1 + endelse + ferr = perr[1] + endif else begin + startp = [cnt, (n_elements(iy)-5)>5, 2., 2., $ + max(col[iy]), 0.] + pars = mpfitfun('erfcfit', iy, col[iy], err, startp, $ + yfit=yfit, perror=perr, /quiet) + if n_elements(perr) eq 6 then begin + fcnt = pars[0] + fwid = pars[1] + endif else begin + perr = 1.e9 + fcnt = 1.e9 + fwid = 1 + endelse + ferr = perr[0] + endelse + ; + ; get data + if ferr le 0.1 and abs(cnt-fcnt) le 1. then begin + data[0,sl,nxp] = float(i) + data[1,sl,nxp] = fcnt + data[2,sl,nxp] = fwid + errs[sl,nxp] = ferr + wids[sl,nxp] = fwid + if fcnt lt minimy then minimy = fcnt + if fcnt gt maximy then maximy = fcnt + nxp += 1 + endif + endif + last = nt +endfor +sl += 1 +; +; log number of slices found +kcwi_print_info,ppar,pre,'Found this many slices',sl,/info +; +; report on statistics +terrs = [-1.] +twids = [-1.] +for isl=0,23 do begin + es = reform(errs[isl,*]) + es = es[where(es ge 0.)] + ws = reform(wids[isl,*]) + ws = ws[where(ws ge 0.)] + emo = moment(es) + wmo = moment(ws) + terrs = [terrs, es] + twids = [twids, ws] + kcwi_print_info,ppar,pre, $ + 'Slice, , STD, min, max, , STD, min, max',isl, $ + emo[0],sqrt(emo[1]),minmax(es), $ + wmo[0],sqrt(wmo[1]),minmax(ws), format='(a,i4,2x,8f9.4)' +endfor +emo = moment(terrs[where(terrs ge 0.)]) +wmo = moment(twids[where(twids ge 0.)]) +kcwi_print_info,ppar,pre,'All: , STD, min, max, , STD, min, max', $ + emo[0],sqrt(emo[1]),minmax(es), $ + wmo[0],sqrt(wmo[1]),minmax(ws), format='(a,6x,8f9.4)' +; +; plot the traces +deepcolor +!p.background=colordex('white') +!p.color=colordex('black') +plot,[0,1],[0,1],xtitle='X pixel',xran=[0,nx],/xs, $ + ytitle='Y pixel',yran=[minimy,maximy], $ + /nodata,title='Direct '+strtrim(kdgeom.ifunam,2)+ $ + ' Slicer Traces for Image # '+strn(kdgeom.arcimgnum) +; +; now get spatial extent of arc images +; +; loop over the slices +for isl = 0, sl-1 do begin + xf = reform(data[0,isl,*]) + yf = reform(data[1,isl,*]) + wf = reform(data[2,isl,*]) + good = where(yf gt 0.) + xf = xf[good] + yf = yf[good] + wf = wf[good] + ; + ; trim outliers + ims,yf,ymn,ysg,ywt + good = where(ywt eq 1, ngood) + if ngood gt 10 then begin + xf = xf[good] + yf = yf[good] + wf = wf[good] + endif + oplot,xf,yf,psym=4 + c = poly_fit(xf,yf,1) + angle = atan(c[1]) / !DTOR + kdgeom.angles[isl] = angle + ims,wf,w,wstd + kdgeom.wid[isl] = fix(w/2.)>1 + ; + ; extract trace from bar image + xx = lindgen(fix(max(xf)-min(xf))) + long(min(xf)) + nxx = n_elements(xx) + x0 = xx[0] + x1 = xx[nxx-1] + y0 = fix(poly(x0,c)+0.5) + y1 = fix(poly(x1,c)+0.5) + if y0 lt y1 then begin + y0 = y0 - w*2. + y1 = y1 + w*2. + endif else begin + yt = y0 + y0 = y1 - w*2. + y1 = yt + w*2. + endelse + kdgeom.x0[isl] = x0 + kdgeom.x1[isl] = x1 + kdgeom.y0[isl] = y0 + kdgeom.y1[isl] = y1 + ; + ; extract subimg + sub = bar[x0:x1,y0:y1] + subr = rot(sub,angle,cubic=-0.5) + yc = fix( ( (y1-y0)+1 ) / 2 ) + bs = median(subr[0:nxx-1,yc-w:yc+w],dimen=2) + bspec[0:nxx-1,isl] = bs + if nxx gt maximx then $ + maximx = nxx +endfor +; +; now get x offsets w.r.t. ref slice +for isl = 0, sl-1 do $ + kdgeom.xoff[isl] = ccpeak(bspec[*,isl],bspec[*,kdgeom.refslice],20) +; +; we get here, all is good +kdgeom.status = 0 +; +; log our change to geom struct +kdgeom.progid = pre +kdgeom.timestamp = systime(1) +; +return +end diff --git a/pro/kcwi_solve_geom.pro b/pro/kcwi_solve_geom.pro new file mode 100644 index 0000000..fee49e4 --- /dev/null +++ b/pro/kcwi_solve_geom.pro @@ -0,0 +1,97 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_SOLVE_GEOM +; +; PURPOSE: +; Solve the wavelength solutions for each arc spectrum +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_SOLVE_GEOM,Spec, Kgeom, Ppar +; +; INPUTS: +; Spec - a array of arc spectra produced by KCWI_EXTRACT_ARCS +; Kgeom - KCWI_GEOM struct from KCWI_TRACE_CBARS and KCWI_EXTRACT_ARCS +; Ppar - KCWI_PPAR pipeline parameter struct +; +; INPUT KEYWORDS: +; Interactive - set to solve arcs interactively +; +; SIDE EFFECTS: +; Modifies KCWI_GEOM struct by calculating new control points that +; take into account the wavelength solution. +; NOTE: sets KGEOM.STATUS to 0 if fitting succeeded, otherwise sets to +; 1 or greater depending on reason for failure (see kcwi_solve_arcs.pro). +; +; PROCEDURE: +; Find the wavelength solution of the reference bar arc and then +; propogate it to the other bars. Record the wavelength solution +; in the wavelength control points in Kgeom. +; +; EXAMPLE: +; Define the geometry from a 'cbars' image and use it to extract and +; display the spectra from an 'arc' image from the same calibration +; sequence. +; +; cbars = mrdfits('image7142_int.fits',0,chdr) +; kcwi_trace_cbars,cbars,Kgeom,/centroid +; arc = mrdfits('image7140_int.fits',0,ahdr) +; kcwi_extract_arcs,arc,kgeom,arcspec +; kcwi_solve_geom,arcspec,kgeom,ppar +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2014-SEP-18 Initial Revision +; 2017-MAY-02 Added interactive option +;- +; +pro kcwi_solve_geom,spec,kgeom,ppar, interactive=interactive, help=help +; +; startup +pre = 'KCWI_SOLVE_GEOM' +q = '' +; +; check inputs +if n_params(0) lt 2 or keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', ArcSpec, Kgeom, Ppar, /INTERACTIVE, /HELP' + return +endif +; +; Check structs +if kcwi_verify_geom(kgeom,/init) ne 0 then return +if kcwi_verify_ppar(ppar,/init) ne 0 then return +; +; check spec +ssz = size(spec) +if ssz[0] ne 2 or ssz[2] ne 120 then begin + kcwi_print_info,ppar,pre,'Input spec array malformed, run KCWI_EXTRACT_ARCS first.',/error + return +endif +; +; plot file +p_fmt = '(i0'+strn(ppar.fdigits)+')' +plfil = ppar.reddir+'wave_cb' + string(kgeom.cbarsimgnum,p_fmt) + $ + '_arc' + string(kgeom.arcimgnum,p_fmt) +; +; do initial fit of central third of ccd +kcwi_fit_center,spec,kgeom,ppar,cntcoeff +; +; solve arc spectra with iterative method +if ppar.waveiter eq 1 then begin + kcwi_solve_arcs_iter,spec,cntcoeff,kgeom,ppar,/tweak,plot_file=plfil +; +; solve arc spectra automagically +endif else begin + kcwi_solve_arcs,spec,cntcoeff,kgeom,ppar,plot_file=plfil +endelse +; +; solve transformation on slice-by-slice basis +kcwi_solve_slices,ppar,kgeom +; +return +end diff --git a/pro/kcwi_stage4geom.pro b/pro/kcwi_stage4geom.pro new file mode 100644 index 0000000..078450e --- /dev/null +++ b/pro/kcwi_stage4geom.pro @@ -0,0 +1,442 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_STAGE4GEOM +; +; PURPOSE: +; This procedure takes the data from basic CCD reduction through the +; geometric correction, which includes solving for wavelength and +; spatial geometries. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_STAGE4GEOM, Procfname, Pparfname +; +; OPTIONAL INPUTS: +; Procfname - input proc filename generated by KCWI_PREP +; defaults to './redux/kcwi.proc' +; Pparfname - input ppar filename generated by KCWI_STAGE2_PREP +; defaults to './redux/kcwi.ppar' +; +; KEYWORDS: +; VERBOSE - set to verbosity level to override value in ppar file +; DISPLAY - set to display level to override value in ppar file +; +; OUTPUTS: +; None +; +; SIDE EFFECTS: +; Outputs processed files in output directory specified by the +; KCWI_PPAR struct read in from Pparfname. +; +; PROCEDURE: +; Reads Pparfname to derive input/output directories and reads the +; corresponding '*.proc' file in output directory to derive the list +; of input files and their associated geometric calibration files. +; +; EXAMPLE: +; Perform stage4geom reductions on the images in 'night1/redux' directory: +; +; KCWI_STAGE4GEOM,'night1/redux/kcwi.ppar' +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-AUG-01 Initial version +; 2014-JAN-30 Handles failure to trace cbars +; 2014-APR-03 Use master ppar and link files +; 2014-MAY-13 Include calibration image numbers in headers +; 2014-SEP-29 Added infrastructure to handle selected processing +; 2015-APR-25 Added CWI flexure hooks (MM) +; 2016-OCT-05 Removed CWI flexure routines +; 2017-MAY-24 Changed to proc control file and removed link file +;- +pro kcwi_stage4geom,procfname,ppfname,help=help,verbose=verbose, display=display + ; + ; setup + pre = 'KCWI_STAGE4GEOM' + startime=systime(1) + q = '' ; for queries + ; + ; help request + if keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', Proc_filespec, Ppar_filespec' + print,pre+': Info - default filespecs usually work (i.e., leave them off)' + return + endif + ; + ; get ppar struct + ppar = kcwi_read_ppar(ppfname) + ; + ; verify ppar + if kcwi_verify_ppar(ppar,/init) ne 0 then begin + print,pre+': Error - pipeline parameter file not initialized: ',ppfname + return + endif + ; + ; directories + if kcwi_verify_dirs(ppar,rawdir,reddir,cdir,ddir) ne 0 then begin + kcwi_print_info,ppar,pre,'Directory error, returning',/error + return + endif + ; + ; check keyword overrides + if n_elements(verbose) eq 1 then $ + ppar.verbose = verbose + if n_elements(display) eq 1 then $ + ppar.display = display + ; + ; log file + lgfil = reddir + 'kcwi_stage4geom.log' + filestamp,lgfil,/arch + openw,ll,lgfil,/get_lun + ppar.loglun = ll + printf,ll,'Log file for run of '+pre+' on '+systime(0) + printf,ll,'DRP Ver: '+kcwi_drp_version() + printf,ll,'Raw dir: '+rawdir + printf,ll,'Reduced dir: '+reddir + printf,ll,'Calib dir: '+cdir + printf,ll,'Data dir: '+ddir + printf,ll,'Filespec: '+ppar.filespec + printf,ll,'Ppar file: '+ppfname + if ppar.clobber then $ + printf,ll,'Clobbering existing images' + printf,ll,'Verbosity level : ',ppar.verbose + printf,ll,'Plot display level: ',ppar.display + ; + ; read proc file + kpars = kcwi_read_proc(ppar,procfname,imgnum,count=nproc) + ; + ; gather configuration data on each observation in reddir + kcwi_print_info,ppar,pre,'Number of input images',nproc + ; + ; loop over images + for i=0,nproc-1 do begin + ; + ; image to process + ; + ; check for flat fielded image first + obfil = kcwi_get_imname(kpars[i],imgnum[i],'_intf',/reduced) + ; + ; if not check for dark subtracted image + if not file_test(obfil) then $ + obfil = kcwi_get_imname(kpars[i],imgnum[i],'_intd',/reduced) + ; + ; if not just get stage1 output image + if not file_test(obfil) then $ + obfil = kcwi_get_imname(kpars[i],imgnum[i],'_int',/reduced) + ; + ; check if input file exists + if file_test(obfil) then begin + ; + ; read configuration + kcfg = kcwi_read_cfg(obfil) + ; + ; direct or dispersed? + if strpos(kcfg.obstype,'direct') ge 0 then $ + do_direct = (1 eq 1) $ + else do_direct = (1 eq 0) + ; + ; final output file + if do_direct then $ + ofil = kcwi_get_imname(kpars[i],imgnum[i],'_img',/reduced) $ + else ofil = kcwi_get_imname(kpars[i],imgnum[i],'_icube',/reduced) + ; + ; get image type + kcfg.imgtype = strtrim(kcfg.imgtype,2) + ; + ; check if output file exists already + if kpars[i].clobber eq 1 or not file_test(ofil) then begin + ; + ; print image summary + kcwi_print_cfgs,kcfg,imsum,/silent + if strlen(imsum) gt 0 then begin + for k=0,1 do junk = gettok(imsum,' ') + imsum = string(i+1,'/',nproc,format='(i3,a1,i3)')+' '+imsum + endif + print,"" + print,imsum + printf,ll,"" + printf,ll,imsum + flush,ll + ; + ; record input file + kcwi_print_info,ppar,pre,'input 2-D image',obfil,format='(a,a)' + ; + ; do we have the geom files? + do_geom = (1 eq 0) ; assume no to begin with + if (strtrim(kpars[i].geomcbar,2) ne '' and $ + strtrim(kpars[i].geomarc,2) ne '') or $ + strtrim(kpars[i].geom,2) ne '' then begin + ; + ; do we have a specified geom file? + if strtrim(kpars[i].geom,2) ne '' then begin + gfile = strtrim(kpars[i].geom,2) + endif + ; + ; do we have specified cbar and arc files (these take precedence) + if strtrim(kpars[i].geomcbar,2) ne '' and $ + strtrim(kpars[i].geomarc,2) ne '' then begin + ; + ; get filenames + cbf = kpars[i].geomcbar + arf = kpars[i].geomarc + ; + ; get corresponding kgeom file + if do_direct then $ + gfile = repstr(arf,'_int.fits','_dgeom.fits') $ + else gfile = repstr(cbf,'_int.fits','_geom.fits') + endif + ; + ; if it exists, read it + if file_test(gfile,/read) then begin + if do_direct then begin + kdgeom = mrdfits(gfile,1,ghdr,/silent) + do_geom = (kdgeom.status eq 0) + endif else begin + kgeom = mrdfits(gfile,1,ghdr,/silent) + do_geom = (kgeom.status eq 0) + endelse + ; + ; log it + kcwi_print_info,ppar,pre,'Using geometry from',gfile,format='(a,a)' + ; + ; if not, derive it + endif else begin + ; + ; time geometry generation + gstartime = systime(1) + ; + ; log + kcwi_print_info,ppar,pre,'Generating geometry solution' + ; + ; read configs + ccfg = kcwi_read_cfg(cbf) + acfg = kcwi_read_cfg(arf) + ; + ; direct image geometry + if do_direct then begin + ; + ; create a new Kdgeom + kdgeom = {kcwi_dgeom} + kdgeom = struct_init(kdgeom) + kdgeom.initialized = 1 + ; + ; populate it with goodness + kcwi_set_dgeom,kdgeom,acfg,kpars[i] + kdgeom.arcfname = arf + kdgeom.arcimgnum = acfg.imgnum + kdgeom.cbarsfname = cbf + kdgeom.cbarsimgnum = ccfg.imgnum + ; + ; get direct geometry + kcwi_solve_dgeom,kdgeom,kpars[i] + ; + ; log bad solution + if kdgeom.status ne 0 then begin + kcwi_print_info,ppar,pre,'bad direct geometry solution',/error + do_geom = (1 eq 0) + endif else $ + do_geom = (1 eq 1) + ; + ; write out result + kcwi_write_geom,kpars[i],kdgeom + ; + ; dispersed image geometry + endif else begin + ; + ; get arc atlas + kcwi_get_atlas,acfg,atlas,atname + ; + ; create a new Kgeom + kgeom = {kcwi_geom} + kgeom = struct_init(kgeom) + kgeom.initialized = 1 + ; + ; populate it with goodness + kcwi_set_geom,kgeom,ccfg,kpars[i],atlas=atlas,atname=atname + kgeom.cbarsfname = cbf + kgeom.cbarsimgnum = ccfg.imgnum + kgeom.arcfname = arf + kgeom.arcimgnum = acfg.imgnum + ; + ; read in cbars image + cbars = mrdfits(cbf,0,chdr,/fscale,/silent) + ; + ; trace the bars + kcwi_trace_cbars,cbars,kgeom,kpars[i],status=stat + ; + ; check status, if < 0 don't proceed + if stat ge 0 then begin + ; + ; log + kcwi_print_info,ppar,pre,'traced continuum bars in cbars image',cbf,format='(a,a)' + ; + ; read in arcs + arc = mrdfits(arf,0,ahdr,/fscale,/silent) + ; + ; extract along bars + kcwi_extract_arcs,arc,kgeom,spec,kpars[i] + ; + ; log + kcwi_print_info,ppar,pre,'extracted arc spectra from arc image',arf,format='(a,a)' + ; + ; do the solution + kcwi_solve_geom,spec,kgeom,kpars[i] + ; + ; log bad solution + if kgeom.status ne 0 then begin + kcwi_print_info,ppar,pre,'bad dispersed geometry solution',/error + do_geom = (1 eq 0) + endif else $ + do_geom = (1 eq 1) + endif else begin + kcwi_print_info,ppar,pre,'unable to trace cont bars',/error + do_geom = (1 eq 0) + endelse + ; + ; write out result + kcwi_write_geom,kpars[i],kgeom + ; + ; output a wavemap image which gives the wavelength at each pixel + kcwi_reverse_geom,kgeom,kpars[i] + endelse ; dispersed image geometry + ; + ; time for geometry + eltime = systime(1) - gstartime + print,'' + printf,ll,'' + kcwi_print_info,ppar,pre,'geom time in seconds',eltime + endelse + ; + ; is our geometry good? + if do_geom then begin + ; + ; read in, update header, apply geometry, write out + ; + ; object image + img = mrdfits(obfil,0,hdr,/fscale,/silent) + ; + sxaddpar,hdr, 'HISTORY',' '+pre+' '+systime(0) + ; + ; apply direct geometry + if do_direct then begin + kcwi_apply_dgeom,img,hdr,kdgeom,kpars[i],dimg,dhdr + ; + ; write out intensity image + ofil = kcwi_get_imname(kpars[i],imgnum[i],'_img',/nodir) + kcwi_write_image,dimg,dhdr,ofil,kpars[i] + ; + ; apply dispersed geometry + endif else begin + kcwi_apply_geom,img,hdr,kgeom,kpars[i],cube,chdr + ; + ; write out intensity cube + ofil = kcwi_get_imname(kpars[i],imgnum[i],'_icube',/nodir) + kcwi_write_image,cube,chdr,ofil,kpars[i] + ; + ; check for arc and output diagnostic 2d image + if strpos(kcfg.imgtype,'arc') ge 0 then begin + rcube = kcwi_get_imname(kpars[i],imgnum[i],'_icube',/reduced) + kcwi_flatten_cube,rcube,/iscale + kcwi_print_info,ppar,pre,'wrote image file', $ + repstr(rcube,'.fits','_2d.fits'), $ + format='(a,a)' + endif + ; + ; variance cube + vfil = repstr(obfil,'_int','_var') + if file_test(vfil,/read) then begin + var = mrdfits(vfil,0,varhdr,/fscale,/silent) + ; + sxaddpar,varhdr,'HISTORY',' '+pre+' '+systime(0) + kcwi_apply_geom,var,varhdr,kgeom,kpars[i],vcub,vchdr + ; + ; write out variance cube + ofil = kcwi_get_imname(kpars[i],imgnum[i],'_vcube',/nodir) + kcwi_write_image,vcub,vchdr,ofil,kpars[i] + endif else $ + kcwi_print_info,ppar,pre,'no variance image found',/warning + ; + ; mask cube + mfil = repstr(obfil,'_int','_msk') + if file_test(mfil,/read) then begin + msk = float(mrdfits(mfil,0,mskhdr,/silent)) + ; + sxaddpar,mskhdr,'HISTORY',' '+pre+' '+systime(0) + kcwi_apply_geom,msk,mskhdr,kgeom,kpars[i],mcub,mchdr + ; + ; write out mask cube + ofil = kcwi_get_imname(kpars[i],imgnum[i],'_mcube',/nodir) + kcwi_write_image,mcub,mchdr,ofil,kpars[i] + endif else $ + kcwi_print_info,ppar,pre,'no mask image found',/warning + ; + ; check for nod-and-shuffle sky images + sfil = kcwi_get_imname(kpars[i],imgnum[i],'_sky',/reduced) + if file_test(sfil,/read) then begin + sky = mrdfits(sfil,0,skyhdr,/fscale,/silent) + ; + sxaddpar,skyhdr,'HISTORY',' '+pre+' '+systime(0) + kcwi_apply_geom,sky,skyhdr,kgeom,kpars[i],scub,schdr + ; + ; write out sky cube + ofil = kcwi_get_imname(kpars[i],imgnum[i],'_scube',/nodir) + kcwi_write_image,scub,schdr,ofil,kpars[i] + endif + ; + ; check for nod-and-shuffle obj images + nfil = kcwi_get_imname(kpars[i],imgnum[i],'_obj',/reduced) + if file_test(nfil,/read) then begin + obj = mrdfits(nfil,0,objhdr,/fscale,/silent) + ; + sxaddpar,objhdr,'HISTORY',' '+pre+' '+systime(0) + kcwi_apply_geom,obj,objhdr,kgeom,kpars[i],ocub,ochdr + ; + ; write out obj cube + ofil = kcwi_get_imname(kpars[i],imgnum[i],'_ocube',/nodir) + kcwi_write_image,ocub,ochdr,ofil,kpars[i] + endif + endelse ; end apply dispersed geometry + ; end if do_geom + endif else $ + kcwi_print_info,ppar,pre,'unusable geom for: '+obfil+' type: '+kcfg.imgtype,/error + ; + ; end check for geom files + endif else begin + ; + ; no problem skipping darks + if strpos(kcfg.imgtype,'dark') ge 0 then $ + kcwi_print_info,ppar,pre,'darks do not get geometry: '+ $ + obfil,/info $ + else kcwi_print_info,ppar,pre,'missing calibration file(s) for: '+ $ + obfil,/warning + endelse + ; + ; end check if output file exists already + endif else begin + kcwi_print_info,ppar,pre,'file not processed: '+obfil+' type: '+kcfg.imgtype,/warning + if kpars[i].clobber eq 0 and file_test(ofil) then $ + kcwi_print_info,ppar,pre,'processed file exists already',/warning + endelse + ; + ; end check if input file exists + endif else $ + kcwi_print_info,ppar,pre,'input file not found: '+obfil,/error + endfor ; loop over images + ; + ; report + eltime = systime(1) - startime + print,'' + printf,ll,'' + kcwi_print_info,ppar,pre,'run time in seconds',eltime + kcwi_print_info,ppar,pre,'finished on '+systime(0) + ; + ; close log file + free_lun,ll + ; + return +end ; kcwi_stage4geom diff --git a/pro/kcwi_trace_cbars.pro b/pro/kcwi_trace_cbars.pro new file mode 100644 index 0000000..131134c --- /dev/null +++ b/pro/kcwi_trace_cbars.pro @@ -0,0 +1,567 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_TRACE_CBARS +; +; PURPOSE: +; This procedure traces the bars in the continuum bars (cbars) +; type calibration images and generates a geometric solution. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_TRACE_CBARS,Img,Kgeom,WarpIm +; +; INPUTS: +; Img - 2-D image of type 'cbars' +; +; INPUT KEYWORDS: +; PWDEG - degree of POLYWARP geometric solution +; PKBUF - number of pixels on each side of peak to use +; STEPSIZE- number of pixel rows to step in tracing peaks in y direction +; AVGROWS - number of rows to average when tracing at each step +; CENTROID- use Centroid only to find peaks (default) +; GAUSS - use Gaussian to fit peaks +; MOFFAT - use Moffat function to fit peaks +; +; OUTPUT KEYWORDS: +; WARPIM - warped version of input image +; +; OUTPUTS: +; Kgeom - KCWI geometric transformation struct +; (see KCWI_GEOM__DEFINE.PRO) +; +; SIDE EFFECTS: +; None. +; +; PROCEDURE: +; Starts by defining peaks from middle row of image. It then +; traces the bars upward and then downward in steps to generate +; a grid of control points that are sent to POLYWARP, which +; generates the Kx, Ky coefficients that are sent to POLY_2D. +; +; NOTE: +; Assumes the input image has already been through stage1 (see +; KCWI_STAGE1.PRO). +; +; EXAMPLE: +; +; Read in a reduced image and generate a warped image and the +; corresponding coefficients using a Moffat function to fit the +; peaks in each row. +; +; img = mrdfits('image7060_int.fits',0,hdr,/fscale) +; kcwi_trace_cbars,img,kgeom,ppar,/moffat +; +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-JUL-12 Initial version +; 2013-JUL-17 Added keywords, condition on warp output +; 2013-JUL-18 Now finds canonical number of peaks (120) +; 2013-JUL-23 Renamed and updated arguments and keywords +; 2013-SEP-14 Use ppar to pass loglun +; 2013-NOV-08 Use median of navg pixels instead of total to reject CRs +; 2014-JAN-30 Handle dome flat scattered light +; 2014-APR-09 Robust check for single-pixel peaks +; 2014-OCT-30 Bar detection now starts at half peak and decriments +;- +forward_function trace_u + +function trace_u, x, p + wid = abs(p[2]) > 1e-20 + return, ((x-p[1])/wid) +end + +pro kcwi_trace_cbars, img, kgeom, ppar, $ + warpim=warpim, pwdeg=pwdeg, pkbuf=pkbuf, stepsize=stepsize, $ + avgrows=avgrows, gauss=gauss, moffat=moffat, centroid=centroid, $ + status=status, help=help +; +; init +pre = 'KCWI_TRACE_CBARS' +q='' +status = -1 +; +; check inputs +if n_params(0) lt 2 or keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', CBarsImg, Kgeom, [Ppar]' + return +endif +if n_params(0) lt 3 then begin + ppar = { kcwi_ppar } + ppar = struct_init(ppar) +endif +; +; check keywords +nter = 4 +if keyword_set(moffat) then nter=5 +pdeg = 3 +if keyword_set(pwdeg) then pdeg=pwdeg +pbuf = 5 +if keyword_set(pkbuf) then pbuf=pkbuf +step = 80/kgeom.ybinsize +if keyword_set(stepsize) then step=stepsize +navg = 6/kgeom.ybinsize +if keyword_set(avgrows) then navg=avgrows +do_plots = ppar.display +; +; prepare plots, if needed +if do_plots ge 1 then begin + if !d.window lt 0 then $ + window,0,title=kcwi_drp_version() $ + else wset,0 + deepcolor + !p.background=colordex('white') + !p.color=colordex('black') +endif +; +; check image size +sz=size(img,/dim) +if sz[0] eq 0 then begin + kcwi_print_info,ppar,pre, $ + 'input image must be 2-D, preferrably a cbars image.', $ + /error + return +endif +; +; check Kgeom +if kcwi_verify_geom(kgeom,/silent) ne 0 then begin +; not a struct so get a shiny, new one + kgeom = { kcwi_geom } + kgeom = struct_init(kgeom) +endif +; +; log +kcwi_print_info,ppar,pre,'Tracing bars in image number',kgeom.cbarsimgnum +; +; number of measurements per peak +npts = sz[1]/step +; +; which rows to measure? +peaky = indgen(npts)*step + fix(step/2) +; +; ensure we don't run off image +while max(peaky) ge sz[1]-navg do begin + npts = npts - 1 + peaky = indgen(npts)*step+fix(step/2) +endwhile +; +; get middle row +print,pre+': Info - Finding peaks in middle row...' +p0 = npts/2 +y0 = peaky[p0] +midrow = y0 +; +row = reform(median(img[*,(y0-navg):(y0+navg)],dimen=2)) +nrpts = n_elements(row) +; +; get rough background vector +win = 50 +backv = row-row +for i=win/2,nrpts-win/2-1 do $ + backv[i] = min(row[i-win/2:i+win/2-1]) +backv[0:win/2-1] = backv[win/2] +backv[nrpts-win/2:*] = backv[nrpts-win/2-1] +; +; subtract +row = row - backv +; +; row maximum +rowmax = max(row) +; +; get sky and skysig +mmm,row,sky,skysig,/silent +; +; check sky value, if fails use brute-force histogram peak +if sky lt min(row) then begin + rhist = histogram(row,min=0.) + t = where(rhist eq max(rhist)) + sky = float(t[0]) +endif +; +; check for zero skysig +if skysig le 0. then $ + skysig = max(row)*0.01 +; +kcwi_print_info,ppar,pre,'Sky, SkySig',sky,skysig +; +; make sure we are above sky + 10.*skysig +skylim = sky + skysig * 10. +; +; find 120 peaks +smul = 0.5 +npks = 0 +tries = 0 +; +; start at half peak, but above sky limit +barth = rowmax*smul>skylim +while npks ne 120 and tries lt 10 do begin + ; + ; find indices for data above threshhold + t=where(row gt barth, nt) + ; + ; make range list with commas splitting each peak + rangepar,t,tstr + sta=strsplit(tstr,',',/extract,count=npks) + ; + ; check for single-pixel lines + ;good = where(strpos(sta,'-') ge 0, ngood) + ;if ngood ne npks then npks = 0 + ; + ; keep decrementing threshhold until we reach 120 peaks + smul -= 0.05 > 0.01 + barth = rowmax*smul>skylim + ; + ; keep incrementing tries until we reach limit + tries += 1 +endwhile +kcwi_print_info,ppar,pre,'final bar thresh, ntries',barth,tries, $ + format='(a,f9.2,i5)' +; +; did we succeed? +if npks ne 120 then begin + kcwi_print_info,ppar,pre,'unable to find 120 peaks',npks,/error + kgeom.status=1 + return +endif +kcwi_print_info,ppar,pre,'Found 120 peaks' +; +; keep track of max values +pkmax = fltarr(npks) +; +; position arrays +peakx = fltarr(npts,npks) - 99.9 +xi = fltarr(npts*npks) ; input x values +yi = fltarr(npts*npks) ; input y values +xo = fltarr(npts*npks) ; output x values + ; output y = input y +barx = fltarr(npks) ; canonical bar x positions +bar = intarr(npts*npks) - 9 ; bar id for each point +slice = intarr(npts*npks) - 9 ; slice id for each point +pp = 0L +; +; get starting positions from middle row +for j=0,npks-1 do begin + ; + ; get index range for this peak + st = sta[j] + ; + ; convert back into index list + rangepar,sta[j],x + ; + ; get fluxes in range + y = row[x] + ; + ; centroid on this peak + cnt=cntrd1d(x,y) + ; + ; re-center window based on centroid + t0 = fix(cnt) - pbuf + t1 = fix(cnt) + pbuf + nx = (t1-t0)+1 + ; + ; get x values + x = t0 + indgen(nx) + ; + ; get fluxes in range + y = row[x] + ; + ; re-centroid before fitting + cnt=cntrd1d(x,y) + ; + ; choose a fit method (default is centroid) + if keyword_set(gauss) or keyword_set(moffat) then begin + ; + ; x-vector for plotting + xp = t0 + indgen(nx*100)/100. + ; + ; initial estimates to improve fitting success ratio + est=[max(y),cnt,2.,1.] + ; + ; do fit + yfit = mpfitpeak(x,y,a,nterms=nter,estimate=est,error=sqrt(y), $ + chisq=chi2,dof=dof,perror=sig,gauss=gauss,moffat=moffat) + ; + ; sometimes mpfitpeak just doesn't return this + if n_elements(sig) le 0 then sig = (a-a)-99. + ; + ; reduced chi^2 + redchi = chi2/float(dof) + ; + ; generate moffat plot + if keyword_set(moffat) then begin + u = trace_u(xp,a) + if nter ge 5 then f = a[4] else f = 0 + if nter ge 6 then f = f + a[5]*xp + denom0 = (u^2 + 1) + denom = denom0^(-a[3]) + yp = f + a[0] * denom + ; + ; else generate gaussian plot + endif else yp = gaussian(xp,a) + ; + ; derive yrange for plot + yrng = [0.,max([yp,y,yfit])]*1.10 + ; + ; store values + barx[j] = a[1] + pkmax[j] = a[0] + peakx[p0,j] = a[1] + xi[pp] = a[1] + xo[pp] = a[1] + ; + ; print results + kcwi_print_info,ppar,pre,'Pk#, chi^2/dof, GauPk, sig, Cntr',$ + j,redchi,a[1],sig[1],cnt, $ + format='(a,i3,4f8.2)',info=2 + ; + ; do centroid instead of fitting + endif else begin + ; + ; get y range for plotting + yrng = [0.,max(y)]*1.10 + ; + ; store values + barx[j] = cnt + pkmax[j] = max(y) + peakx[p0,j] = cnt + xi[pp] = cnt + xo[pp] = cnt + ; + ; print results + kcwi_print_info,ppar,pre,'Pk#, Cntr',j,cnt, $ + format='(a,i3,f8.2)',info=2 + endelse + ; + ; values independant of fitting method + yi[pp] = y0 + bar[pp] = j + slice[pp] = fix(j/5) + pp += 1 + ; + ; plot results + if do_plots ge 2 then begin + xrng = [min(x)-1,max(x)+1] + plot,x,y,psym=-4,xran=xrng,/xsty,xtitle='PIX', $ + yran=yrng,/ysty,ytitle='INT', $ + charsi=1.5,symsi=1.5, $ + title = 'Image: '+strn(kgeom.cbarsimgnum) + $ + ', Middle Row - Bar: '+strn(j)+ $ + ', Slice: '+strn(fix(j/5)) + ; + ; over plot fit results if fitting done + if keyword_set(gauss) or keyword_set(moffat) then begin + oplot,xp,yp,linesty=5 + oplot,[a[1],a[1]],[-100,1e9],linesty=5 + endif + ; + ; overplot centroid + oplot,[cnt,cnt],[-100,1e9] + ; + ; query user + read,'Next? (Q-quit plotting, -next): ',q + if strupcase(strmid(strtrim(q,2),0,1)) eq 'Q' then $ + do_plots = 0 + endif +endfor +print,pre+': Info - Done.' +; +; recover plot state +do_plots = ppar.display +; +; calculate reference output x control point positions +; assumes kgeom.refbar is middle bar of slice (out of 5 bars per slice) +if kgeom.refbar ge 0 and kgeom.refbar lt 120 then $ + refb = kgeom.refbar $ +else refb = 57 ; default +delx = 0. +; +; maybe use i=2,4 to avoid vignetted first bar? +for i=1,4 do begin + bb = (refb-2) + i + delx = delx + barx[bb] - barx[bb-1] +endfor +delx = delx / 4.0 +kgeom.refdelx = delx +kgeom.refoutx = barx[refb] + (findgen(5)-2.) * delx +kgeom.x0out = fix(delx/2.) + 1 +; +; log results +kcwi_print_info,ppar,pre,'reference delta x',kgeom.refdelx +kcwi_print_info,ppar,pre,'reference x positions',kgeom.refoutx, $ + format='(a,5f9.2)' +; +; now trace each peak up and down +if do_plots ge 1 then begin + ; + ; plot control points over entire image + xrng = [0,sz[0]] + yrng = [0,sz[1]] + ; + ; start with middle row + plot,peakx[p0,*],replicate(peaky[p0],npks), $ + xran=xrng,/xsty,xtitle='X PIX', $ + yran=yrng,/ysty,ytitle='Y PIX', $ + title='Image: '+strn(kgeom.cbarsimgnum) + $ + ', CBARS Input Geometry Control Points', $ + charsi=1.5,symsi=1.5,psym=7 +endif +; +; trace both ways +msg = [ 'Tracing peaks in upper half of image...', $ + 'Tracing peaks in lower half of image...' ] +inc = [1,-1] +off = [-1,1] +lim = [npts-1,0] +; +; loop over up,down +for k=0,1 do begin + kcwi_print_info,ppar,pre,msg[k] + ; + ; loop over peaks + for j=0,npks-1 do begin + ; + ; loop over samples + for i=p0+inc[k],lim[k],inc[k] do begin + x0c = peakx[i+off[k],j] + ; + ; allow for skipping some bad fits + nskip = 1 + while x0c lt 0 and nskip lt 3 do begin ; can only skip 3 + nskip = nskip + 1 + x0c = peakx[i+nskip*off[k],j] + endwhile + ; + ; make sure we have a good starting point + if x0c gt 0 then begin + ; + ; working row + iy = peaky[i] + row=reform(total(img[*,(iy-navg):(iy+navg)],2)/float(2*navg+1.)) + ; + ; initial x range + t0 = fix(x0c) - pbuf + t1 = fix(x0c) + pbuf + nx = (t1-t0)+1 + x = t0 + indgen(nx) + ; + ; get fluxes in range + y = row[x] + ; + ; get centroid + cnt=cntrd1d(x,y) + ; + ; re-do x range + t0 = fix(cnt) - pbuf + t1 = fix(cnt) + pbuf + nx = (t1-t0)+1 + x = t0 + indgen(nx) + ; + ; update fluxes in range + y = row[x] + ; + ; get updated centroid + cnt=cntrd1d(x,y) + ; + ; do fitting if we are not centroiding + if keyword_set(gauss) or keyword_set(moffat) then begin + ; + ; initial estimates to improve fitting success ratio + est=[max(y),cnt,2.,1.] + ; + ; do fit + yfit = mpfitpeak(x,y,a,nterms=nter,estimate=est, $ + error=sqrt(y),chisq=chi2,dof=dof,sigma=sig, $ + gauss=gauss,moffat=moffat) + ; + ; sometimes mfitpeak just doesn't return this + if n_elements(sig) le 0 then sig = (a-a)-99. + ; + ; reduced chi^2 for testing + redchi = chi2/float(dof) + ; + ; check fit + if finite(redchi) eq 1 and redchi lt 199. and $ + redchi gt 0. and a[0] gt pkmax[j]*0.05 and $ + sig[1] lt 99. and sig[1] gt 0. and $ + abs(cnt-a[1]) lt 1.5 then begin + if do_plots ge 1 then plots,a[1],iy,psym=1 + ; + ; store values + peakx[i,j] = a[1] + xi[pp] = a[1] + xo[pp] = peakx[p0,j] + endif + ; + ; print results + ;kcwi_print_info,ppar,pre, $ + ; 'Bar#,Sam#,chi^2/DOF,GauCnt,GauSig,Cnt,Pk', $ + ; j,i,redchi,a[1],sig[1],cnt,a[0], $ + ; format='(a,2i5,5f9.2)',info=2 + ; + ; here we are just centroiding so no fitting + endif else begin + ; + ; check centroid + if abs(x0c - cnt) le 2.0 and $ + max(y) gt pkmax[j]*0.05 then begin + if do_plots ge 1 then plots,cnt,iy,psym=1 + ; + ; store values + peakx[i,j] = cnt + xi[pp] = cnt + xo[pp] = peakx[p0,j] + endif + ; + ; print results + ;kcwi_print_info,ppar,pre,'Bar#,Sam#,Cnt,Max[y],Pk', $ + ; j,i,cnt,max(y),pkmax[j], $ + ; format='(a,2i5,3f9.2)',info=2 + endelse ; centroiding + ; + ; values independant of fitting method + yi[pp] = iy + bar[pp] = j + slice[pp] = fix(j/5) + pp += 1 + endif ; x0c gt 0 + endfor ; looping over samples + endfor ; looping over peaks +endfor ; loop over directions (up,down) +print,pre+': Info - Done.' +; +; continue? +if do_plots ge 2 then $ + read,'Continue (): ',q +; +; yo is the same as yi +yo = yi +; +; call polywarp +polywarp,xi,yi,xo,yo,pdeg,fkx,fky,/double +; +; fill KCWI_GEOM struct +kgeom.kx = fkx +kgeom.ky = fky +kgeom.xi = xi +kgeom.yi = yi +kgeom.xo = xo +kgeom.yo = yo +kgeom.bar = bar +kgeom.slice = slice +kgeom.barx = barx +kgeom.midrow = midrow +kgeom.initialized = 1 +kgeom.progid = pre +kgeom.timestamp = systime(1) +; +; make warped version, if requested +warpim = poly_2d(img,fkx,fky,2,cubic=-0.5) +; +status = 0 ; we're good! +return +end diff --git a/pro/kcwi_verify_cfg.pro b/pro/kcwi_verify_cfg.pro new file mode 100644 index 0000000..b231a6a --- /dev/null +++ b/pro/kcwi_verify_cfg.pro @@ -0,0 +1,57 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_VERIFY_CFG +; +; PURPOSE: +; This function verifies the input KCWI_CFG struct. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; Result = KCWI_VERIFY_CFG(Kcfg) +; +; INPUTS: +; Kcfg - array of struct KCWI_CFG +; +; RETURNS: +; The status of the input KCWI_CFG struct: +; 0 - verified without problems +; 1 - a malformed KCWI_CFG struct was passed +; +; KEYWORDS: +; INITIALIZED - set to check if KCWI_CFG is initialized +; SILENT - set to silence output +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-JUL-08 Initial version +; 2013-SEP-14 Added initialized keyword +;- +function kcwi_verify_cfg,kcfg,initialized=initialized,silent=silent + ; + ; setup + pre = 'KCWI_VERIFY_CFG' + ; + ; check input + stat = 0 + sz = size(kcfg) + if sz[0] ne 1 or sz[1] lt 1 or sz[2] ne 8 then begin + if not keyword_set(silent) then $ + print,pre+': Error - malformed KCWI_CFG struct array' + stat = 1 + endif else begin + if keyword_set(initialized) then begin + test = n_elements(kcfg) + if total(kcfg.initialized) ne test then begin + print,pre+': Error - KCWI_CFG struct not initialized' + stat = 1 + endif + endif + endelse + ; + return,stat +end diff --git a/pro/kcwi_verify_dirs.pro b/pro/kcwi_verify_dirs.pro new file mode 100644 index 0000000..c841dfd --- /dev/null +++ b/pro/kcwi_verify_dirs.pro @@ -0,0 +1,120 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_VERIFY_DIRS +; +; PURPOSE: +; This function verifies the input KCWI_PPAR struct. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; Result = KCWI_VERIFY_DIRS(Ppar,Rawdir,Reddir,Cdir,Ddir) +; +; INPUTS: +; Ppar - array of struct KCWI_PPAR +; +; RETURNS: +; -1 - error +; 0 - all dirs exist and accessible +; 1 - raw dir not accessible +; 2 - reduced dir not accessible +; 3 - calib dir not accessible +; 4 - data dir not accessible +; +; OPTIONAL OUTPUTS: +; Rawdir - raw image directory +; Reddir - reduced image directory +; Cdir - calib directory +; Ddir - data directory +; +; KEYWORDS: +; NOCREATE - set to prevent creation of output directory +; READONLY - set to skip test for writablity +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-SEP-16 Initial version +; 2014-APR-03 Added calib dir check +; 2015-FEB-20 Renamed in/out dirs to raw/red +;- +function kcwi_verify_dirs,ppar,rawdir,reddir,cdir,ddir, $ + nocreate=nocreate, readonly=readonly + ; + ; setup + pre = 'KCWI_VERIFY_DIRS' + q = '' + ; + ; check input + if kcwi_verify_ppar(ppar,/init) ne 0 then return,-1 + ; + ; directories + rawdir = kcwi_expand_dir(ppar.rawdir) + reddir = kcwi_expand_dir(ppar.reddir) + cdir = kcwi_expand_dir(ppar.caldir) + ddir = kcwi_expand_dir(ppar.datdir) + ; + ; check if rawdir exists and is readable + if not file_test(rawdir,/directory,/executable,/read) then begin + kcwi_print_info,ppar,pre,'cannot access raw image dir',rawdir,/error + return,1 + endif + ; + ; check if reddir exists + if not file_test(reddir,/directory) then begin + if not keyword_set(nocreate) then begin + print,pre+': Warning - reduced image dir does not exist: ',reddir + read,'Create? (Y/n): ',q + q = strupcase(strtrim(q,2)) + if strmid(q,0,1) ne 'N' then begin + file_mkdir,reddir,/noexpand + kcwi_print_info,ppar,pre,'created directory',reddir + endif + endif else begin + kcwi_print_info,ppar,pre,'no reduced image dir',/error + return,2 + endelse + endif + ; + ; check if reddir accessible + if keyword_set(readonly) then begin + if not file_test(reddir,/directory,/executable,/read) then begin + kcwi_print_info,ppar,pre, $ + 'reduced image dir not readable',/error + return,2 + endif + endif else begin + if not file_test(reddir,/directory,/executable,/read,/write) then begin + kcwi_print_info,ppar,pre, $ + 'reduced image dir not read/writable',/error + return,2 + endif + endelse + ; + ; check if cdir accessible + if keyword_set(readonly) then begin + if not file_test(cdir,/directory,/executable,/read) then begin + kcwi_print_info,ppar,pre, $ + 'calib dir not readable',/error + return,3 + endif + endif else begin + if not file_test(cdir,/directory,/executable,/read,/write) then begin + kcwi_print_info,ppar,pre, $ + 'calib dir not read/writable',/error + return,3 + endif + endelse + ; + ; check if ddir exists and is readable + if not file_test(ddir,/directory,/executable,/read) then begin + kcwi_print_info,ppar,pre,'cannot access data dir',ddir,/error + return,4 + endif + ; + ; if we get here all is well + return,0 +end diff --git a/pro/kcwi_verify_geom.pro b/pro/kcwi_verify_geom.pro new file mode 100644 index 0000000..57db5c9 --- /dev/null +++ b/pro/kcwi_verify_geom.pro @@ -0,0 +1,56 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_VERIFY_GEOM +; +; PURPOSE: +; This function verifies the input KCWI_GEOM struct. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; Result = KCWI_VERIFY_GEOM(Kgeom) +; +; INPUTS: +; Kgeom - KCWI_GEOM struct +; +; RETURNS: +; The status of the input KCWI_GEOM struct: +; 0 - verified without problems +; 1 - a malformed or uninitialized KCWI_GEOM struct was passed +; +; KEYWORDS: +; INITIALIZED - set to check if KCWI_GEOM struct is initialized +; SILENT - set to silence output +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-SEP-15 Initial version +;- +function kcwi_verify_geom,kgeom,initialized=initialized,silent=silent + ; + ; setup + pre = 'KCWI_VERIFY_GEOM' + ; + ; check input + stat = 0 + sz = size(kgeom) + if sz[0] ne 1 or sz[1] lt 1 or sz[2] ne 8 then begin + if not keyword_set(silent) then $ + print,pre+': Error - malformed KCWI_[D]GEOM struct' + stat = 1 + endif else begin + if keyword_set(initialized) then begin + if kgeom.initialized ne 1 then begin + if not keyword_set(silent) then $ + print,pre+': Error - KCWI_[D]GEOM struct not initialized, run geometry fitting program first' + stat = 1 + endif + endif + endelse + ; + return,stat +end diff --git a/pro/kcwi_verify_ppar.pro b/pro/kcwi_verify_ppar.pro new file mode 100644 index 0000000..83f2a84 --- /dev/null +++ b/pro/kcwi_verify_ppar.pro @@ -0,0 +1,56 @@ +; +; Copyright (c) 2013, California Institute of Technology. All rights +; reserved. +;+ +; NAME: +; KCWI_VERIFY_PPAR +; +; PURPOSE: +; This function verifies the input KCWI_PPAR struct. +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; Result = KCWI_VERIFY_PPAR(Ppar) +; +; INPUTS: +; Ppar - array of struct KCWI_PPAR +; +; RETURNS: +; The status of the input KCWI_PPAR struct: +; 0 - verified without problems +; 1 - a malformed or uninitialized KCWI_PPAR struct was passed +; +; KEYWORDS: +; INITIALIZED - set to check if KCWI_PPAR struct is initialized +; SILENT - set to silence output +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2013-SEP-13 Initial version +; 2013-SEP-14 Added initialized keyword +;- +function kcwi_verify_ppar,ppar,initialized=initialized,silent=silent + ; + ; setup + pre = 'KCWI_VERIFY_PPAR' + ; + ; check input + stat = 0 + sz = size(ppar) + if sz[0] ne 1 or sz[1] lt 1 or sz[2] ne 8 then begin + if not keyword_set(silent) then $ + print,pre+': Error - malformed KCWI_PPAR struct array' + stat = 1 + endif else begin + if keyword_set(initialized) then begin + if ppar.initialized ne 1 then begin + print,pre+': Error - KCWI_PPAR struct not initialized' + stat = 1 + endif + endif + endelse + ; + return,stat +end diff --git a/pro/kcwi_write_dgeom.pro b/pro/kcwi_write_dgeom.pro new file mode 100644 index 0000000..ba612fd --- /dev/null +++ b/pro/kcwi_write_dgeom.pro @@ -0,0 +1,82 @@ +; +; Copyright (c) 2016, California Institute of Technology. All rights reserved. +;+ +; NAME: +; KCWI_WRITE_DGEOM +; +; PURPOSE: +; Writes out the KCWI_DGEOM struct as an IDL save file +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_WRITE_DGEOM,Ppar,Kdgeom +; +; INPUTS: +; Kdgeom - KCWI_DGEOM struct +; Ppar - KCWI_PPAR pipeline parameter struct +; +; INPUT KEYWORDS: +; TEST - only write out if kdgeom.status=0 (good fit) +; +; OUTPUTS: +; None. +; +; PROCEDURE: +; Uses the tag kdgeom.dgeomfile to write out the struct as an IDL +; save file. Checks if ppar.clobber is set and takes appropriate +; action. +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2014-SEP-11 Initial Revision +;- +; +pro kcwi_write_dgeom,ppar,kdgeom, test=test +; +; startup +pre = 'KCWI_WRITE_DGEOM' +q = '' +; +; check inputs +if n_params(0) lt 1 or keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', Ppar, Kdgeom' + return +endif +; +; Check structs +if kcwi_verify_geom(kdgeom,/init) ne 0 then return +if kcwi_verify_ppar(ppar,/init) ne 0 then return +; +; check fit status +if keyword_set(test) and kdgeom.status ne 0 then begin + kcwi_print_info,ppar,pre,'Kdgeom fit no good, nothing written.',/error + return +endif +; +; write it out +; check if it exists already +if file_test(kdgeom.dgeomfile) then begin + ; + ; clobber it, if requested + if ppar.clobber eq 1 then begin + file_delete,kdgeom.dgeomfile,verbose=ppar.verbose + kcwi_print_info,ppar,pre,'deleted existing dgeom file', $ + kdgeom.dgeomfile,format='(a,a)' + endif else begin + kcwi_print_info,ppar,pre, $ + 'existing dgeom file undisturbed', $ + kdgeom.dgeomfile,format='(a,a)' + return + endelse +endif +; +; write it out if we get here +save,kdgeom,filename=kdgeom.dgeomfile +kcwi_print_info,ppar,pre,'wrote dgeom file',kdgeom.dgeomfile,format='(a,a)' +; +return +end diff --git a/pro/kcwi_write_geom.pro b/pro/kcwi_write_geom.pro new file mode 100644 index 0000000..cc3a7b4 --- /dev/null +++ b/pro/kcwi_write_geom.pro @@ -0,0 +1,100 @@ +; +; Copyright (c) 2014, California Institute of Technology. All rights reserved. +;+ +; NAME: +; KCWI_WRITE_GEOM +; +; PURPOSE: +; Writes out the KCWI_GEOM struct as a FITS file +; +; CATEGORY: +; Data reduction for the Keck Cosmic Web Imager (KCWI). +; +; CALLING SEQUENCE: +; KCWI_WRITE_GEOM,Ppar,Kgeom +; +; INPUTS: +; Kgeom - KCWI_GEOM struct +; Ppar - KCWI_PPAR pipeline parameter struct +; +; INPUT KEYWORDS: +; TEST - only write out if kgeom.status=0 (good fit) +; +; OUTPUTS: +; None. +; +; PROCEDURE: +; Uses the tag kgeom.geomfile to write out the struct as a FITS +; table file. Checks if ppar.clobber is set and takes appropriate +; action. +; +; EXAMPLE: +; +; MODIFICATION HISTORY: +; Written by: Don Neill (neill@caltech.edu) +; 2014-SEP-11 Initial Revision +; 2017-JUN-29 Converted to writing FITS instead of save file +;- +; +pro kcwi_write_geom,ppar,kgeom, test=test +; +; startup +pre = 'KCWI_WRITE_GEOM' +q = '' +; +; check inputs +if n_params(0) lt 1 or keyword_set(help) then begin + print,pre+': Info - Usage: '+pre+', Ppar, Kgeom' + return +endif +; +; Check structs +if kcwi_verify_geom(kgeom,/init) ne 0 then return +if kcwi_verify_ppar(ppar,/init) ne 0 then return +; +; check fit status +if keyword_set(test) and kgeom.status ne 0 then begin + kcwi_print_info,ppar,pre,'Kgeom fit no good, nothing written.',/error + return +endif +; +; write it out +; check if it exists already +if file_test(kgeom.geomfile) then begin + ; + ; clobber it, if requested + if ppar.clobber eq 1 then begin + file_delete,kgeom.geomfile,verbose=ppar.verbose + kcwi_print_info,ppar,pre,'deleted existing geom file', $ + kgeom.geomfile,format='(a,a)' + endif else begin + kcwi_print_info,ppar,pre, $ + 'existing geom file undisturbed', $ + kgeom.geomfile,format='(a,a)' + return + endelse +endif +; +; get header +if n_elements(tag_names(kgeom)) lt 90 then $ + hdr = headfits(kgeom.arcfname) $ ; direct +else hdr = headfits(kgeom.cbarsfname) ; dispersed +sxaddpar,hdr,'SIMPLE','F' +sxdelpar,hdr,'NAXIS' +sxdelpar,hdr,'NAXIS1' +sxdelpar,hdr,'NAXIS2' +sxdelpar,hdr,'BITPIX' +sxdelpar,hdr,'BZERO' +sxdelpar,hdr,'BSCALE' +; +; write it out if we get here +mwrfits,kgeom,kgeom.geomfile +; +; update header +modfits,kgeom.geomfile,0,hdr +; +; log +kcwi_print_info,ppar,pre,'wrote geom file',kgeom.geomfile,format='(a,a)' +; +return +end From 6579270ecfd8360625f38842dda9ef77b76ae415 Mon Sep 17 00:00:00 2001 From: Zheng Cai Date: Thu, 2 Nov 2017 16:54:06 -0700 Subject: [PATCH 02/11] find_center.py calculate the center of the QSO position for each target. --- pyp_kcwi/find_center.py | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 pyp_kcwi/find_center.py diff --git a/pyp_kcwi/find_center.py b/pyp_kcwi/find_center.py new file mode 100644 index 0000000..bf5b3a7 --- /dev/null +++ b/pyp_kcwi/find_center.py @@ -0,0 +1,41 @@ +def find_center(image): + ''' + Find the center of QSOs for each target, + we will return the center of the QSO; + and later, we will use the center to + calculate the offset. (Z. Cai) + + input: datacube + ''' + # construct a 2D broadband image using median stacking along the wavelength + sci_med=np.median(image, axis=0) + + # interpolate to a finer grid: + x_ori= np.linspace(0, len(sci_med[:,0])-1, len(sci_med[:,0])) # 0--23 (24 pixels) + y_ori= np.linspace(0, len(sci_med[0,:])-1, len(sci_med[0,:])) # 0--69 (70 pixels) + x_new= np.linspace(0, len(sci_med[:,0])-1, 3.*len(sci_med[:,0])-2) # new grid, 3x pixel number + y_new= np.linspace(0, len(sci_med[0,:])-1, 3.*len(sci_med[0,:])-2) # new grid, 3x pixel number + + f = interpolate.interp2d(y_ori, x_ori, sci_med,kind='cubic') + sci_interp = f(y_new,x_new) + + # find center using pixels > 0.93* the peak value: + indices = np.where(sci_med > 0.93*np.amax(sci_med)) + indices_interp = np.where(sci_interp > 0.93*np.amax(sci_interp)) + + # center using the original pixel grid + x_mean0= np.mean(indices[0]) + y_mean0= np.mean(indices[1]) + + # center using the new finer pixel grid (1/3 x original pixel scale) + x_mean_interp= np.mean(indices_interp[0])*1/3. + y_mean_interp= np.mean(indices_interp[1])*1/3. + + # return center determined using both the original and new pixel grid. + center= [x_mean0, y_mean0, x_mean_interp,y_mean_interp] + return center + + +img1, ih1 = fits.getdata("kb171021_00077_icuber.fits", header=True) +center1=find_center(img1) + From 2204f2d0972db4ece68372b235fc9e7c3eea4ad5 Mon Sep 17 00:00:00 2001 From: Zheng Cai Date: Thu, 2 Nov 2017 16:54:51 -0700 Subject: [PATCH 03/11] Update find_center.py --- pyp_kcwi/find_center.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/pyp_kcwi/find_center.py b/pyp_kcwi/find_center.py index bf5b3a7..1a767f8 100644 --- a/pyp_kcwi/find_center.py +++ b/pyp_kcwi/find_center.py @@ -1,3 +1,18 @@ +# This function is used to subtract sky median +# +# INPUT: - sky_file +# - img_file +# - mask_limit mask limit +# - cut_ch cut channel number +# +# INPUT EXAMPLE: +# +# img= fits.getdata("kb171021_00092_icuber.fits", header=True) +# centers= find_center(img) +# +# MODIFICATION HISTORY: +# 2017-11-01 Initial version (ZC) + def find_center(image): ''' Find the center of QSOs for each target, From 8b25cf3ea91128ae57a0f675a027be310d8c16dc Mon Sep 17 00:00:00 2001 From: Zheng Cai Date: Thu, 2 Nov 2017 17:05:05 -0700 Subject: [PATCH 04/11] Create cube_shift.py --- pyp_kcwi/cube_shift.py | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 pyp_kcwi/cube_shift.py diff --git a/pyp_kcwi/cube_shift.py b/pyp_kcwi/cube_shift.py new file mode 100644 index 0000000..3340f3d --- /dev/null +++ b/pyp_kcwi/cube_shift.py @@ -0,0 +1,42 @@ +# This function is used to offset a datacube +# +# INPUT: - img_file +# - offset x, offset y (in pixels) +# - offset can be determined using find_center.py +# +# OUTPUT: - shifted datacube for combining +# +# INPUT EXAMPLE: +# +# img, ih = fits.getdata("kb171021_00077_icuber.fits", header=True) +# img_shift= cube_shift(img,shift_x,shift_y) +# +# EXAMPLE 2: +# +# img1, ih1 = fits.getdata("kb171021_00077_icuber.fits", header=True) +# img2, ih2 = fits.getdata("kb171021_00082_icuber.fits", header=True) +# +# center1=find_center(img1) +# center2=find_center(img2) +# shift_x= center2[0]-center1[0] +# shift_y= center2[1]-center1[1] +# +# img2_shift= cube_shift(img2,shift_x,shift_y) +# center_new= find_center(img2_shift) +#'''One can see new center for img2 is the same as that of img1. Ready for combining''' +# +# +# MODIFICATION HISTORY: +# 2017-11-02 Initial version (ZC) +# + +def cube_shift(image,shift_x,shift_y): + + image_sh= image # construct a new datacube called image_sh (image_shift) + for i in np.arange(0, len(image[:,1,1]), 1): # loop layers + # offset each layer using shift_x, shift_y + image_sh[i,:,:]= np.roll(image[i,:,:], -int(round(shift_x)), axis=0) + image_sh[i,:,:]= np.roll(image[i,:,:], -int(round(shift_y)), axis=1) + + return image_sh + From 64dcf33a9182f8e1443d25241134250fb14ae2e9 Mon Sep 17 00:00:00 2001 From: Zheng Cai Date: Thu, 2 Nov 2017 17:05:56 -0700 Subject: [PATCH 05/11] Update find_center.py --- pyp_kcwi/find_center.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pyp_kcwi/find_center.py b/pyp_kcwi/find_center.py index 1a767f8..8a811e2 100644 --- a/pyp_kcwi/find_center.py +++ b/pyp_kcwi/find_center.py @@ -1,9 +1,8 @@ -# This function is used to subtract sky median +# This function is used to determine the center of the QSO for each datacube # -# INPUT: - sky_file +# INPUT: # - img_file -# - mask_limit mask limit -# - cut_ch cut channel number +# # # INPUT EXAMPLE: # From 109043a296b724d6a00cbae96254b07b884c1ccb Mon Sep 17 00:00:00 2001 From: Zheng Cai Date: Thu, 2 Nov 2017 17:07:56 -0700 Subject: [PATCH 06/11] Update find_center.py --- pyp_kcwi/find_center.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyp_kcwi/find_center.py b/pyp_kcwi/find_center.py index 8a811e2..123bb56 100644 --- a/pyp_kcwi/find_center.py +++ b/pyp_kcwi/find_center.py @@ -50,6 +50,4 @@ def find_center(image): return center -img1, ih1 = fits.getdata("kb171021_00077_icuber.fits", header=True) -center1=find_center(img1) From ca896f4d65e0f36d796c95852d0e25a80b698db5 Mon Sep 17 00:00:00 2001 From: Zheng Cai Date: Fri, 3 Nov 2017 11:46:33 -0700 Subject: [PATCH 07/11] Update cube_shift.py --- pyp_kcwi/cube_shift.py | 65 +++++++++++++++++++++--------------------- 1 file changed, 33 insertions(+), 32 deletions(-) diff --git a/pyp_kcwi/cube_shift.py b/pyp_kcwi/cube_shift.py index 3340f3d..aab6049 100644 --- a/pyp_kcwi/cube_shift.py +++ b/pyp_kcwi/cube_shift.py @@ -1,37 +1,38 @@ -# This function is used to offset a datacube -# -# INPUT: - img_file -# - offset x, offset y (in pixels) -# - offset can be determined using find_center.py -# -# OUTPUT: - shifted datacube for combining -# -# INPUT EXAMPLE: -# -# img, ih = fits.getdata("kb171021_00077_icuber.fits", header=True) -# img_shift= cube_shift(img,shift_x,shift_y) -# -# EXAMPLE 2: -# -# img1, ih1 = fits.getdata("kb171021_00077_icuber.fits", header=True) -# img2, ih2 = fits.getdata("kb171021_00082_icuber.fits", header=True) -# -# center1=find_center(img1) -# center2=find_center(img2) -# shift_x= center2[0]-center1[0] -# shift_y= center2[1]-center1[1] -# -# img2_shift= cube_shift(img2,shift_x,shift_y) -# center_new= find_center(img2_shift) -#'''One can see new center for img2 is the same as that of img1. Ready for combining''' -# -# -# MODIFICATION HISTORY: -# 2017-11-02 Initial version (ZC) -# + def cube_shift(image,shift_x,shift_y): - +''' + This function is used to offset a datacube + + INPUT: - img_file + - offset x, offset y (in pixels) + - offset can be determined using find_center.py + + OUTPUT: - shifted datacube for combining + + INPUT EXAMPLE: + + img, ih = fits.getdata("kb171021_00077_icuber.fits", header=True) + img_shift= cube_shift(img,shift_x,shift_y) + + EXAMPLE 2: + + img1, ih1 = fits.getdata("kb171021_00077_icuber.fits", header=True) + img2, ih2 = fits.getdata("kb171021_00082_icuber.fits", header=True) + + center1=find_center(img1) + center2=find_center(img2) + shift_x= center2[0]-center1[0] + shift_y= center2[1]-center1[1] + + img2_shift= cube_shift(img2,shift_x,shift_y) + center_new= find_center(img2_shift) +Then One can see new center for img2 is the same as that of img1. Ready for combining + + + MODIFICATION HISTORY: + 2017-11-02 Initial version (ZC) +''' image_sh= image # construct a new datacube called image_sh (image_shift) for i in np.arange(0, len(image[:,1,1]), 1): # loop layers # offset each layer using shift_x, shift_y From 9e697ceb020397519c0b06cd32a506bb02647568 Mon Sep 17 00:00:00 2001 From: Zheng Cai Date: Fri, 3 Nov 2017 11:49:55 -0700 Subject: [PATCH 08/11] Update find_center.py --- pyp_kcwi/find_center.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/pyp_kcwi/find_center.py b/pyp_kcwi/find_center.py index 123bb56..18bd3a7 100644 --- a/pyp_kcwi/find_center.py +++ b/pyp_kcwi/find_center.py @@ -19,7 +19,15 @@ def find_center(image): and later, we will use the center to calculate the offset. (Z. Cai) - input: datacube + This function is used to determine the center of the QSO for each datacube + + INPUT EXAMPLE: + + img= fits.getdata("kb171021_00092_icuber.fits", header=True) + centers= find_center(img) + + MODIFICATION HISTORY: + 2017-11-01 Initial version (ZC) ''' # construct a 2D broadband image using median stacking along the wavelength sci_med=np.median(image, axis=0) From 04096a185e646d3063f96790596560340bc352d9 Mon Sep 17 00:00:00 2001 From: Zheng Cai Date: Fri, 3 Nov 2017 11:50:09 -0700 Subject: [PATCH 09/11] Update find_center.py --- pyp_kcwi/find_center.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/pyp_kcwi/find_center.py b/pyp_kcwi/find_center.py index 18bd3a7..5b9b79b 100644 --- a/pyp_kcwi/find_center.py +++ b/pyp_kcwi/find_center.py @@ -1,16 +1,3 @@ -# This function is used to determine the center of the QSO for each datacube -# -# INPUT: -# - img_file -# -# -# INPUT EXAMPLE: -# -# img= fits.getdata("kb171021_00092_icuber.fits", header=True) -# centers= find_center(img) -# -# MODIFICATION HISTORY: -# 2017-11-01 Initial version (ZC) def find_center(image): ''' From 7582bbef8de718513c65ee2c7b04c2aaffe9f975 Mon Sep 17 00:00:00 2001 From: Qiong Li Date: Wed, 8 Nov 2017 18:35:36 -0800 Subject: [PATCH 10/11] Add files via upload simple test for resampling --- resample.pro | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 resample.pro diff --git a/resample.pro b/resample.pro new file mode 100644 index 0000000..ad887b7 --- /dev/null +++ b/resample.pro @@ -0,0 +1,33 @@ +; the simplest example to understand resampling image + +pro resample + +; Create the display window +WINDOW, /FREE, XSIZE=20000, YSIZE=10000 + +; Set up the arrays of original points +yo = INDGEN(100) +xo = fix(sqrt(yo)) + +; Create and display the original image: +img = FLTARR(100,100) +xt1 = fix(sqrt(yo))+30 +xt2 = fix(sqrt(yo))+60 +img[yo,xo]=1 +img[yo,xt1]=1 +img[yo,xt2]=1 + +; X and Y coordinates to be fit +xi = fix(FLTARR(100) + median(xo)) +yi = yo + +; Run POLYWARP to obtain a Kx and Ky: +polywarp, xi, yi, xo, yo, 3, kx, ky + +; Create a warped image based on Kx and Ky with POLY_2D: +wimg = poly_2d(img, kx, ky, 2,cubic=-0.5) + +TVSCL, img, 0 +TVSCL, wimg, 1 + +end \ No newline at end of file From c3200e69c5b73a9dc1b38447be062c415b8713fa Mon Sep 17 00:00:00 2001 From: Qiong Li Date: Sat, 11 Nov 2017 22:57:44 -0800 Subject: [PATCH 11/11] Update cube_shift.py line 36: if a=b, a and b will change together. It's better to avoid using it. line 40: 'image_sh' not 'image' --- pyp_kcwi/cube_shift.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyp_kcwi/cube_shift.py b/pyp_kcwi/cube_shift.py index aab6049..c36482c 100644 --- a/pyp_kcwi/cube_shift.py +++ b/pyp_kcwi/cube_shift.py @@ -33,11 +33,13 @@ def cube_shift(image,shift_x,shift_y): MODIFICATION HISTORY: 2017-11-02 Initial version (ZC) ''' - image_sh= image # construct a new datacube called image_sh (image_shift) + ###image_sh= image # construct a new datacube called image_sh (image_shift) + image_sh= np.copy(image) # construct a new datacube called image_sh (image_shift) for i in np.arange(0, len(image[:,1,1]), 1): # loop layers # offset each layer using shift_x, shift_y image_sh[i,:,:]= np.roll(image[i,:,:], -int(round(shift_x)), axis=0) - image_sh[i,:,:]= np.roll(image[i,:,:], -int(round(shift_y)), axis=1) + ###image_sh[i,:,:]= np.roll(image[i,:,:], -int(round(shift_y)), axis=1) + image_sh[i,:,:]= np.roll(image_sh[i,:,:], -int(round(shift_y)), axis=1) return image_sh