pro mk_sfd, infil, outdir, prefix, interactive=interactive, run_time=run_time, $ ist=ist, filename=filename, outfil=outfil, $ qdebug=qdebug, qstop=qstop ; ;+ ;NAME: ; mk_sfd ;PURPOSE: ; Given a set of files, read the roadmaps and find where there ; is a long and short exposure. Make a composite image and remove ; the saturated pixels. Write images to an output file. ;INPUT: ; infil - an array of input file names ; outdir - the output directory where the file should be written ;OPTIONAL KEYWORD INPUT: ; interactive- If set, display the results to the screen as the ; processing occurs. ; outfil - The output file name (excluding the directory) ;KEYWORD OUTPUT: ; filename- the output file name as derived from the input file list ;HISTORY: ; draft version: 28 January 1992 (KTS) ; 20-Mar-92 (MDM) - Took Keith Strong's program "db_desat.pro" ; and made the procedure "mk_desat" ; 7-Apr-92 (KTS) - Put IFLAG in to avoid error if first image ; is rejected ; 14-Apr-92 (KTS) - 'PREFIX' parameter hardwired for 'sfd' ; - Added automatic week # feature ; - Better generated dark frame algorithm ; - Avoid night images algorithm ;Ver 2.0 16-Apr-92 (KTS) - Added dark frame library read ; 16-Apr-92 (MDM) - Changed file name definition to use the ; "progverno" variable to determine the extension ; - Put the prefix parameter back in (so it would work ; with old existing programs) ; 20-Apr-92 (MDM) - Added RUN_TIME parameter ; 21-Apr-92 (KTS) - Included LWA/MDM desat technique ; - Included MDM DC routine ;Ver 2.1 22-Apr-92 (MDM) - Added "roadmap" to the RD_SDA calls (speeds ; things up since it will not re-read the roadmap ; for each call) ; - Adjusted the format of the code so MDM can ; read it more easily ; - Fixed DC subtraction to use the whole image, ; not just one pixel ; - Put the call to GET_DC_IMAGE for the short ; exposure images outside the loop. ; - Moved the check for apropriate short exposure ; to before the RD_SDA (use the roadmap info) ; - Replaced RD_SDA with RD_XDA ; 23-Apr-92 (MDM) - Added "ist" input option to skip to image ; "ist" in the loop (for recovery from bombout) ;Ver 2.2 26-Apr-92 (MDM) - Corrected bug in background subtraction. For ; quarter res images, it was subtracting them ; wrong. data1-dark1(*,*,i) does not work if ; dark1 is 512x512xN and data1 is 256x256. But ; data1-dark1 will work even if dark1 is 512x512! ; - Renamed from "mk_desat" to "mk_sfd" ;Ver 2.21 26-Apr-92 (MDM) - Adjusted search to look at the image maximum ; and only choose images with max greater than 150 ; DN (high-8). There are images which are being ; taken after the sun has set, so they are dark ; images. ;Ver 2.22 26-Apr-92 (MDM) - Added the warm dark current image subtraction ; capability (calling "get_dc_warm") ;Ver 2.23 27-Apr-92 (MDM) - Fixed bug (induced by me) where the "m1" ; search was being performed on the exposure ; normalized short image. It was grabbing too ; many pixels. ;Ver 2.24 27-Apr-92 (MDM) - Added a check to see that the dark image ; subtraction is at least close when taking ; warm images. ;Ver 2.25 30-Jul-92 (MDM) - Check to see how many minutes before predicted ; night the image was taken. Exclude images taken ; less than 2 minute before night. Removed the ; img_max > 150 check (tried to do the same thing) ;Ver 3.00 17-Aug-92 (MDM) - Updated version number to 3.00 to make a new ; file extension because of 2.25 changes that ; result in a reduced number of images selected. ; - Added "filename" parameter ;Ver 3.01 24-Aug-92 (MDM) - Passed PROGVERNO and PROGNAME to SAV_SDA ;Ver 3.02 28-Aug-92 (MDM) - Return without doing anything if cannot find ; any proper images ; - Passed PROVERNO*1000 to SAV_SDA ;Ver 3.03 20-Oct-92 (MDM) - Added "OUTFIL" input parameter ;Ver 3.04 21-Oct-92 (MDM) - Modified selection of images to process ; changed exposure range and technique used ; OLD: Used ROADMAP.EXPLEVMODE directly (MBE) ; NEW: Use GT_DPE ; OLD Long Exposure Bracket: DPE 17-23 (MBE 10-16) ; OLD Short Exposure Bracket: DPE 9-14 (MBE 2- 7) ; New Long Exposure Bracket: DPE 17-25 ; New Short Exposure Bracket: DPE 9-15 ; Range changed because dagwood short/long exposures ; where changed from DPE 13/23 to 15/25 around 7-Oct. ;Ver 3.05 22-Oct-92 (MDM) - Futhur mod/fix from Ver 3.04 (missed one WHERE) ; 7-Jan-93 (MDM) - Added QDEBUG and QSTOP keywords ;Ver 4.00 7-Jan-93 (MDM) - Made change to GT_EXPDUR which corrected for ; returning a zero exposure duration (was causing ; saturated images) ; - Put in a call to LEAK_SUB to remove the thin Al ; pin hole scattered light leak ; - Put in code to try to eliminate images with bleed ; from insufficient flushes (from previous staturated ; image) ;Ver 4.01 8-Jan-93 (MDM) - Corrected code to remove bleed images to do so ; in the short exposure images also ;Ver 5.00 14-Jan-93 (MDM) - New patch to the LEAK_SUB routine to work properly ; with quarter resolution images ;Ver 6.00 3-Feb-93 (MDM) - LEAK_SUB had been updated to correct for AlMg leak ; also. Started use of /SAVE option for LEAK_SUB. ; ** BOGUS ** - see Ver 7.00 below ; 8-Feb-93 (MDM) - Trapped on case where files exist, but are empty ;Ver 7.00 26-Mar-93 (MDM) - Enabled using LEAK_SUB for AlMg filters ;Ver 8.00 21-Jul-93 (MDM) - Changed the insertion of the short exposure for ; pixels over 90 DN to only be done where the long ; exposure is saturated. ; ** This fixes the glitch in the SXL histograms ** ; 23-Jul-93 (SLF) - Use temporary function for memory management ;Ver 9.00 26-Aug-93 (MDM) - Enabled history records ; - Enabled interpolation for the dark images ; - Updated the SFC leak image to be used (was ; using Nov-92, now uses Aug-93 image) ;Ver 9.01 31-Aug-93 (MDM) - The "Ver 2.25 30-Jul-92" patch to check to see ; how many minutes before predicted night was not ; done in the final selection of the short exposure ; "where" command. It is fixed now. ;V 10.00 29-Nov-93 (MDM) - Select only images taken with Filter A open ; (don't select neutral density images) ;V 11.00 3-Dec-93 (MDM) - Added rough alignment of short and long exposures ;V 11.01 3-Jan-94 (MDM) - Corrected a bug which caused subscript out of range ; (was looking at m1 instead of m1s for -1) ; 24-Jan-94 (DMZ) - removed VMS version numbers from infile ;V 12.00 12-apr-95 (SLF) - Filter out fill data images ;V 12.01 9-oct-95 (SLF) - Filter out non-square images ;- ; if (n_elements(prefix) eq 0) then prefix='sfd' ;-- strip off VMS version numbers since they seem to cause ; problems further on down the line if !version.os eq 'vms' then begin break_file,infil,dsk,direc,file,ext infil=dsk+direc+file+ext endif qdebug = keyword_set(qdebug) qdebug2 = keyword_set(qstop) start_time = systime(1) run_time = 0. progname = 'MK_SFD' progverno=12.00 ;program Version number iflag=0 ;Flag to prevent error if first ;image is rejected at late stage if (infil(0) eq '') then return ;no SFR files if (n_elements(outdir) eq 0) then outdir = '' rd_roadmap, infil, roadmap, filidx=filidx if (get_nbytes(roadmap) lt 10) then return ;no data in the files - added 8-Feb-93 tim2orbit, roadmap, tim2night=tim2night if (keyword_set(interactive)) then window,0,xs=512,ys=512,retain=2 ; Find number of elements and set up any required arrays for longer ; exposures (omits 15 sec at present) ssla=where(roadmap.percentd eq 255 and gt_filtb(roadmap) ge 2 and $ gt_filtb(roadmap) le 3 and $ gt_filta(roadmap) eq 1 and $ ;mdm added 29-Nov-93 gt_res(roadmap) ge 1 and gt_expmode(roadmap) eq 0 and $ gt_dpe(roadmap) ge 17 and gt_dpe(roadmap) le 25 and $ (gt_shape_cmd(roadmap,/x) eq gt_shape_cmd(roadmap,/y)) and $ ; slf 9-oct tim2night gt 2.0) ; roadmap.img_max gt 150) ;;roadmap.explevmode ge 10 and roadmap.explevmode le 16 and $ ;removed 21-Oct-92 if (ssla(0) eq -1) then return ;no long exposure matches get_dc_image, roadmap(ssla), dci1, dark1, imap, /interp ; get dark frames for long exposure ;--- Beware that there is another WHERE command below sssa=where(roadmap.percentd eq 255 and gt_filtb(roadmap) ge 2 and $ ;MDM start adding 22-Apr-92 gt_filtb(roadmap) le 3 and $ gt_filta(roadmap) eq 1 and $ ;mdm added 29-Nov-93 gt_res(roadmap) ge 1 and gt_expmode(roadmap) eq 0 and $ gt_dpe(roadmap) ge 9 and gt_dpe(roadmap) le 15 and $ (gt_shape_cmd(roadmap,/x) eq gt_shape_cmd(roadmap,/y)) and $ ; slf 9-oct tim2night gt 2.0) ;;roadmap.explevmode ge 2 and roadmap.explevmode le 7 and $ ;removed 21-Oct-92 if (sssa(0) eq -1) then return ;no short exposure matches get_dc_image, roadmap(sssa), dci2, dark2, imap2, /interp ; get dark frames for short exposures imap2_arr = intarr(n_elements(roadmap)) imap2_arr(sssa) = imap2 ;MDM end adding 22-Apr-92 dla=mk_dset_str(filidx, ssla) tla=int2secarr(roadmap(ssla),roadmap(0)) nla=n_elements(tla) ; Read files sequentially and find nearest short exposure if (keyword_set(ist)) then iflag = 1 ;file already started if (n_elements(ist) eq 0) then ist = 0 for i = ist,nla-1 do begin time1_all = systime(1) rr = roadmap(ssla(i)) print, i+1, ' of ', nla,' images ', fmt_tim(rr), gt_res(rr,space=3) ; Print which frame we are on ss=where(roadmap.percentd eq 255 and $ gt_filtb(roadmap) eq gt_filtb(rr) and $ gt_filta(roadmap) eq 1 and $ ;mdm added 29-Nov-93 gt_res(roadmap) eq gt_res(rr) and $ gt_expmode(roadmap) eq 0 and $ (gt_shape_cmd(roadmap,/x) eq gt_shape_cmd(roadmap,/y)) and $ ; slf 9-oct gt_dpe(roadmap) ge 9 and gt_dpe(roadmap) le 15 and $ ;mdm added 22-Oct-92 tim2night gt 2.0) ;mdm added here 31-Aug-93 ;;roadmap.explevmode ge 2 and roadmap.explevmode le 7) ;mdm removed 22-oct-92 if (ss(0) lt 0) then goto, lab1 ;no short exp - loop ds=mk_dset_str(filidx,ss) ; Find the closest time of a short exposure tss = int2secarr(roadmap(ss),roadmap(0)) dt = abs(tss - tla(i)) ctm = min(dt) if (ctm gt 600.) then goto, lab1 ;too long - loop id=where(dt eq ctm) ;id is in reference to ss subset time1 = systime(1) rd_xda, infil, ds(id(0)), index2, data2, roadmap time2 = systime(1) if (qdebug) then print, 'RD_XDA for short exposure took ', (time2-time1), ' seconds' time1 = systime(1) rd_xda, infil, dla(i), index1, data1, roadmap ;read next long exp. time2 = systime(1) if (qdebug) then print, 'RD_XDA for long exposure took ', (time2-time1), ' seconds' siz = size(data1) ;added 26-Apr-92 (MDM) nx1 = siz(1) ny1 = siz(2) siz = size(data2) nx2 = siz(1) ny2 = siz(2) isum = 2^gt_res(index1) ;1,2,4 ;added Saturation check 7-Jan-93 x0 = 450/isum x1 = 550/isum y0 = 50/isum ;lower center of CCD readout y1 = 150/isum tmp = data1(x0:x1, y0:y1) if (total(tmp)/n_elements(tmp) gt 250) then begin print, 'MK_SFD: Evidence of bleed back from insufficient flush' print, 'MK_SFD: Skipping this image goto, lab1 end tmp = data2(x0:x1, y0:y1) if (total(tmp)/n_elements(tmp) gt 250) then begin print, 'MK_SFD: Evidence of bleed back from insufficient flush' print, 'MK_SFD: Skipping this image goto, lab1 end if ((histogram(data1))(170) ge gt_shape_cmd(index1,/x)) or $ ((histogram(data2))(170) ge gt_shape_cmd(index2,/x)) then begin print,'MK_SFD: Fill data limit exceeded' print,'MK_SFD: Skipping this image file_append,'filldata.dat',[get_info(index1,/non),get_info(index2,/non)] goto, lab1 endif time1 = systime(2) datal=sxt_decomp(data1) ; decompress short and long exposures datas=sxt_decomp(data2) time2 = systime(2) if (qdebug) then print, 'Decompression for two images took ', (time2-time1), ' seconds' ; Short exposure taken at night? - if so then loop ; if (max(datas) lt 5000.) then goto, lab1 ; if (max(datal) lt 5000.) then goto, lab1 ; Subtract background and normalize exposures time1 = systime(2) if (gt_temp_ccd(index1) lt -20) then begin dark10 = dark1(0:nx1-1, 0:ny1-1, imap(i)) ;needed because if 4x4 is nested in 2x2, the subtraction is done ;wrong with the command a=a-dark1(*,*,i) ;dark1 was already decompressed up above dci10 = dci1(imap(i)) end else begin ;get_dc_warm, index1, dci10, dark10 get_dc_image, index1, dci10, dark10, /interp end if (gt_temp_ccd(index2) lt -20) then begin dark20 = dark2(0:nx2-1,0:ny2-1,imap2_arr(ss(id(0)))) ;MDM added 22-Apr-92 ;dark2 was already decompressed up above dci20 = dci2(imap2_arr(ss(id(0)))) end else begin get_dc_image, index2, dci20, dark20, /interp end ;datal=datal-dark1(imap(i)) ;this was just subtracting a scalar value!! - MDM removed 22-Apr-92 ;datas=datas-back2d datal=datal-dark10 datas=datas-dark20 his_index, /enable his_index, index1, 0, 'time_dark', dci10.his.time_dark ;only save the long exposure dark subtraction information his_index, index1, 0, 'day_dark', dci10.his.day_dark time_compos = [gt_time(index1), gt_time(index2), 0] day_compos = [gt_day(index1), gt_day(index2), 0] his_index, index1, 0, 'q_composite', progverno*1000 his_index, index1, 0, 'time_compos', time_compos his_index, index1, 0, 'day_compos', day_compos if (gt_temp_ccd(index1) gt -20) then begin ;MDM added check that warm dark current subtraction was at least close tmp = datal(0:40,0:40) tmp = total(tmp)/n_elements(tmp) if (abs(tmp) gt 100) then goto, lab1 end ;if ( (gt_filtb(index1) eq 2) and (int2secarr(index1, '13-Nov-92 17:06') gt 0) ) then begin ;added 7-Jan-92 if (int2secarr(index1, '13-Nov-92 17:06') gt 0) then begin ;modified 26-Mar-93 datal = leak_sub(index1, datal, /save, /update_index) datas = leak_sub(index2, datas, /save, /update_index) end ;-------------------- Added rough registration of short and long images (MDM) 3-Dec-93 cen = sxt_cen([index1.gen,index2.gen]) dx = (cen(0,0) - cen(0,1)) / (2^gt_res(index1)) ;they are the same resolution dy = (cen(1,0) - cen(1,1)) / (2^gt_res(index1)) dx = fix(dx + 0.5*sgn(dx)) ;number of binned pixels to shift the short image over by dy = fix(dy + 0.5*sgn(dy)) ; ;positive dx ==> long exposure (index1) is to the right of index2 ;==> when extracting from the short exposure, we need to subtract from the subscripts print, 'Shifting short exposure by', dx, dy, ' original pixels' ;m1 = where(datas ge 90) ;find good pixels in short m1s = where( (datas ge 90) and (data1 ge 253)) ;find good pixels in short - Change made 21-Jul-93 (Ver 8.00) m1 = (m1s + dx + dy*nx1)>0<(nx1*ny1-1) ;do a simple shift since the saturated pixels are never ;on the edge of the ccd - would get wrap around in that case ;(MDM added 3-Dec-93) ;m1 is the pixel address in the long exposure, m1s is in the short exposure time2 = systime(2) if (qdebug) then print, 'Dark image subtraction took ', (time2-time1), ' seconds' if (qdebug) then begin print, 'Long: ' + fmt_tim(index1) + gt_res(index1,sp=2) + gt_dpe(index1,/conv,/str) + $ ' Dark: ' + fmt_tim(dci10) + gt_res(dci10,sp=2) + gt_dpe(dci10,/conv,/str) print, 'Short: ' + fmt_tim(index2) + gt_res(index2,sp=2) + gt_dpe(index2,/conv,/str) + $ ' Dark: ' + fmt_tim(dci20) + gt_res(dci20,sp=2) + gt_dpe(dci20,/conv,/str) end time1 = systime(2) datal=temporary(datal)*(1000./gt_expdur(index1)) ;slf, 23-Jul-93 datas=temporary(datas)*(1000./gt_expdur(index2)) ;use temporary function datal=temporary(datal)>1 ; datas=temporary(datas)>1 ; time2 = systime(2) if (qdebug) then print, 'Image normalization took ', (time2-time1), ' seconds' time1 = systime(2) s2=where(data1 ge 253) ;find sat. pixels s2s = (s2 - dx - dy*nx1)>0<(nx1*ny1-1) ;do a simple shift since the saturated pixels are never ;on the edge of the ccd (MDM added 3-Dec-93) ;--- MDM replaced case statement which check if it was HR or QR and divided by 512 or 256 - use "nx1" xx= s2 mod nx1 yy = s2/nx1 s2up = ((yy+1)<(nx1-1))*nx1 + xx s2dn = ((yy-1)>0 )*nx1 + xx s2sup = (s2up - dx - dy*nx1)>0<(nx1*ny1-1) s2sdn = (s2dn - dx - dy*nx1)>0<(nx1*ny1-1) ; Put the longest exposure in the working array ; and replace saturated pixels with smoothed med image ;;cimg=datal ;MDM removed and renamed all references to "cimg" to "datal" ;save memory/computing time if (s2(0) ne -1) then begin ;MDM added 22-Apr-92 (some images have no saturated pixels) datal(s2) = datas(s2s) datal(s2up) = datas(s2sup) datal(s2dn) = datas(s2sdn) timg=smooth(datal,5) datal(s2) = timg(s2) datal(s2up) = timg(s2up) datal(s2dn) = timg(s2dn) end if (m1s(0) ne -1) then datal(m1)=datas(m1s) ; Replace "m1" pixels in datal with unsmoothed med exposure values ; normalize for different resolutions if (gt_res(index1) eq 2) then datal=datal*0.25 ; compress the images to byte arrays print,'Maximum Count Rate: ',max(datal) datal=alog10(datal>1.) datal=bytscl(datal,max=6.,min=0.,top=255) print,'Maximum Count Rate: ',max(datal) time2 = systime(2) if (qdebug) then print, 'Desaturation/image combination/image scaling took ', (time2-time1), ' seconds' if (keyword_set(interactive)) then begin if (gt_res(index1) eq 2) then begin tv,rebin(datal,512,512,/sample) endif else begin tv,datal endelse text= string(i) + ' ' + fmt_tim(index1) + ' ' + $ gt_filtb(index1,/string) + ' ' + gt_res(index1,/string) xyouts,10,10,text,/device end time1 = systime(2) if ((iflag eq 0) or (i eq ist)) then begin anum=strlen(infil(0)) ;length of file id fileid = getobsid(startt=strmid(infil(0),anum-11,6)) ;create Week ID fileid=fileid + '.' + string(fix(progverno), format='(i2.2)') filename = concat_dir(outdir, prefix+fileid) end if (keyword_set(outfil)) then filename = concat_dir(outdir, outfil) sav_sda,filename,index1,datal,append=iflag, progname=progname, progverno=progverno*1000 iflag=1 ;reset flag time2 = systime(2) if (qdebug) then print, 'SAV_SDA took ', (time2-time1), ' seconds' if (qdebug2) then stop time2_all = systime(1) if (qdebug) then print, '*********************** One image process time = ', (time2_all-time1_all), ' seconds lab1: dumy=0 endfor end_time = systime(1) run_time = (end_time-start_time)/60. end