;+ ; NAME: ; iris_sgflat_calc_fuv() ; ; PURPOSE: ; For making the FUV spectral flat field. ; Takes the temporal average of the data produced by ; iris_sgflat_prep and does the forward geometry tranformation. ; The average spatial and spectral patterns are found and ; transformed with the reverse transformation and applied to the ; original image to produce the flat field. ; ; NOTE: This requires two geometry files, one each for the FUV 1 and ; FUV 2 channels ; ; CALLING SEQUENCE: ; result=iris_sgflat_calc_fuv(filename=filename, [/pickfile, ; datadir=, /geofile]) ; ; INPUTS: ; filename - an IDL save file with the Level 1.5 data in the form ; of index, data ; ; OPTIONAL INPUTS: ; ; KEYWORD PARAMETERS: ; /pickfile - instead of supplying a filename variable, the user ; can pick the file starting from the datadir path ; datadir - data directory for temporary and permanent storage ; of the data. If nothing is specified then it uses ; the current working directory ; geofile - user can supply the string path and filename of an ; IDL save file produced by iris_spec_cal, or if set ; as a keyword the user can interactively select the file ; ; OUTPUTS: ; result - filename of the FITS file containing the results ; ; OPTIONAL OUTPUTS: ; None ; ; COMMON BLOCKS: ; None ; ; PROCEDURES USED: ; ; COMMENTS: ; ; EXAMPLES: ; IDL> result1=iris_bgcal_getdata(time_start, time_end, index, data, $ ; > datadir=datadir, /nosave, /fuv) ; > result2=iris_bgcal_prep(index, data, datadir=datadir) ; > result3=iris_bgcal_calc(filename=result2, datadir=datadir) ; ; MODIFICATION HISTORY: ; Started 2014-June-26 by Sarah A. Jaeggli, Montana State University ;- function iris_sgflat_calc_fuv, index, data, filename=filename, $ datadir=datadir, geofile=geofile, $ pickfile=pickfile if keyword_set(datadir) ne 1 then datadir='' if keyword_set(pickfile) then $ filename=dialog_pickfile(path=datadir, filter='*FUV_l1p5.sav', $ title='Please select a L1.5 file to process') if keyword_set(filename) then restore, filename if typename(geofile) eq 'INT' then $ geofile=dialog_pickfile(path=datadir, filter='*_spec_cal.sav', $ /multiple_files, $ title='Please select FUV geometry files') if index[0].img_path ne 'FUV' then begin print, 'Data provided is not FUV, exiting...' return, 0 endif rotidx=7 ;produce a mean flat image flat=median(data, dimension=3) nx=(size(data))[1] ny=(size(data))[2] ;set up buffer for master flat results master_flat=fltarr(nx,ny) master_spat_flat=fltarr(nx,ny) master_spat_flat_rm=fltarr(nx,ny) spec_profile=fltarr(nx) spat_profile=fltarr(ny,2) ;set up buffer for fiducial coordinates fx_coord=dblarr(2,2) fy_coord=dblarr(2,2) ;restore ccd flat ccd_flat_file = concat_dir(concat_dir('$SSW_IRIS','data'), $ 'FUV_flat.fits') fxread, ccd_flat_file, ccd_flat ccd_flat=rotate(ccd_flat, rotidx) for k=0,1 do begin ccd_side=strmatch(geofile[k], '*FUV2*') case ccd_side of 0 : begin x0=0 & x1=2071 dx0=x0 & dx1=x1 xmask0=34 & xmask1=2045 ymask0=20 & ymask1=1060 ;1075 end 1 : begin x0=2072 & x1=4143 dx0=x0 & dx1=x1 xmask0=700 & xmask1=2058 ymask0=20 & ymask1=1075 end endcase nx=x1-x0+1 ;restore geometry data ;geometry correction is measured from first pixel on each *CCD* restore, geofile[k] ;geo_solution ;calculate forward destortion map xo=rebin(dindgen(nx, 1), nx, ny) yo=rebin(dindgen(1, ny), nx, ny) xi_forward=replicate(0.D,nx,ny) yi_forward=replicate(0.D,nx,ny) for i=0,2 do begin for j=0,2 do begin xi_forward+=geo_solution.kx[i,j] * xo^j * yo^i yi_forward+=geo_solution.ky[i,j] * xo^j * yo^i endfor endfor ;calculate reverse distortion map xi_reverse=replicate(0.D,nx,ny) yi_reverse=replicate(0.D,nx,ny) for i=0,2 do begin for j=0,2 do begin xi_reverse+=geo_solution.kxi[i,j] * xo^j * yo^i yi_reverse+=geo_solution.kyi[i,j] * xo^j * yo^i endfor endfor ;mask overscan/bad pixels nanmask=replicate(!values.f_nan, nx, ny) nanmask[xmask0:xmask1,ymask0:ymask1]=1 flat1=flat[x0:x1,*]/ccd_flat[x0:x1,*]*nanmask fill_missing, flat1, !values.f_nan, 1, /extrapolate fill_missing, flat1, !values.f_nan, 2, /extrapolate ;upsample the image scl=2D flat2=iris_sgflat_fftrebin(flat1, scl) ;upsample the nanmask bigmask=rebin(nanmask, nx*scl, ny*scl) ;scale the geometry results xi_forward=rebin(xi_forward, nx*scl, ny*scl)*scl yi_forward=rebin(yi_forward, nx*scl, ny*scl)*scl xi_reverse=rebin(xi_reverse, nx*scl, ny*scl)*scl yi_reverse=rebin(yi_reverse, nx*scl, ny*scl)*scl ;transform image to straight coordinates flat2=interpolate(flat2, xi_forward, yi_forward, cubic=-0.5) bigmask=interpolate(bigmask, xi_forward, yi_forward, cubic=-0.5) ;Build spectral flat spec_flat0=rebin(mean(flat2*bigmask, dimension=2, /nan), nx*scl, ny*scl) spec_profile[x0:x1]=rebin(median(flat2*bigmask, dimension=2), nx) spec_flat = smooth(flat2*bigmask, [0,ny*scl/12.], /edge_truncate, /nan) ;Build spatial flat ;only consider positions with signal bad=where(spec_flat0*bigmask lt 10) ;pick some reasonable threshold in intensity linemask=bigmask linemask[bad]=!values.f_nan spat_flat = mean(flat2/spec_flat*linemask, dimension=1, /nan) spat_profile[*,k]=rebin(mean(flat2*linemask, dimension=1, /nan), ny) ;fit and remove fiducials spat_flat_rm=spat_flat ;estimate the fiducial position yarray=dindgen(ny*scl) df=60 ;allowed fit range for f=0,1 do begin a=geo_solution.fiducials.polyfit[f,0:1] fpos=scl*(a[0] + a[1] * nx/4.) fit=gaussfit(yarray[fpos-df:fpos+df], $ spat_flat[fpos-df:fpos+df], a, nterms=6) spat_flat_rm-=a[0]*exp(-0.5*((yarray-a[1])/a[2])^2) fx_coord[k,f]=nx/4.+x0 fy_coord[k,f]=a[1]/scl endfor ;sort coordinates for later fx_coord[k,*]=fx_coord[k,sort(fy_coord[k,*])] fy_coord[k,*]=fy_coord[k,sort(fy_coord[k,*])] spat_flat = rebin(reform(spat_flat, 1, ny*scl), nx*scl, ny*scl) spat_flat_rm = rebin(reform(spat_flat_rm, 1, ny*scl), nx*scl, ny*scl) ;transform images back to curved coordinates ;and bin back down to size spat_flat=interpolate(spat_flat, xi_reverse, yi_reverse, cubic=-0.5) spat_flat_rm=interpolate(spat_flat_rm, xi_reverse, yi_reverse, cubic=-0.5) spec_flat=interpolate(spec_flat, xi_reverse, yi_reverse, cubic=-0.5) linemask=interpolate(linemask, xi_reverse, yi_reverse, cubic=-0.5) spat_flat=congrid(spat_flat, nx, ny, cubic=-0.5) spat_flat_rm=congrid(spat_flat_rm, nx, ny, cubic=-0.5) spec_flat=congrid(spec_flat, nx, ny, cubic=-0.5) linemask=congrid(linemask, nx, ny, cubic=-0.5) ;divide and multiply to make master flat tmp_flat=flat[x0:x1,*] / spec_flat / spat_flat * nanmask * linemask tmp_flat[where(finite(tmp_flat) ne 1)]=1.0 master_flat[x0:x1,*]=tmp_flat tmp_spat_flat=spat_flat * nanmask tmp_spat_flat[where(finite(tmp_spat_flat) ne 1)]=1.0 master_spat_flat[x0:x1,*]=tmp_spat_flat tmp_spat_flat=spat_flat_rm * nanmask tmp_spat_flat[where(finite(tmp_spat_flat) ne 1)]=1.0 master_spat_flat_rm[x0:x1,*]=tmp_spat_flat endfor ;display results imexam, master_flat*master_spat_flat, 0, 1.2 ;write files dt=index[0].date_obs timestring=strmid(dt, 0,4) + $ ;year strmid(dt, 5,2) + $ ;month strmid(dt, 8,2) + $ ;day '_' + $ strmid(dt,11,2) + $ ;hour strmid(dt,14,2) + $ ;min strmid(dt,17,2) ;sec savename=datadir+timestring+'_FUV_flat.sav' fitsname=datadir+timestring+'_FUV_flat.fits' print, 'Saving FUV flat data in '+savename save, flat, ccd_flat, master_flat, master_spat_flat, master_spat_flat_rm, $ fx_coord, fy_coord, spec_profile, spat_profile, filename=savename ;build a fits file mkhdr, hdr, [[[ccd_flat]], [[master_spat_flat]], [[master_spat_flat_rm]]] spawn, 'echo $USER', creator spawn, 'echo $HOSTNAME', host sxaddpar, hdr, 'CREATOR', creator[0]+'@'+host[0] sxaddpar, hdr, 'FID1_X0', fx_coord[0,0] sxaddpar, hdr, 'FID1_X1', fx_coord[0,1] sxaddpar, hdr, 'FID1_Y0', fy_coord[0,0] sxaddpar, hdr, 'FID1_Y1', fy_coord[0,1] sxaddpar, hdr, 'FID2_X0', fx_coord[1,0] sxaddpar, hdr, 'FID2_X1', fx_coord[1,1] sxaddpar, hdr, 'FID2_Y0', fy_coord[1,0] sxaddpar, hdr, 'FID2_Y1', fy_coord[1,1] print, 'Saving FUV flats to fits file '+fitsname writefits, fitsname, $ [[[ccd_flat]], [[master_spat_flat]], [[master_spat_flat_rm]]], $ hdr return, fitsname end