;+ ; NAME: ; eis_prep ; ; PURPOSE: ; ; CATEGORY: ; Hansteen/Wikstol Data analysis SW ; ; CALLING SEQUENCE: ; ; INPUTS: ; ; input:Name of FITS file with EIS data to be calibrated. If the file ; is already open a data object can be sent instead of a filename. ; ; KEYWORD PARAMETERS: ; ; -- Dark current data -- ; default: Use the default file containing dark current data. ; ; dc_file: Optional name of the dark current file. ; ; -- Saving -- ; save: Save the data and the error to files. The data is saved ; to "eis_l1_YYYYMMDD_HHMMSS.fits" The error is saved to ; "eis_er_YYYYMMDD_HHMMSS.fits". The default is to save to ; the current directory. ; ; outdir: Optional directory for storing the fits files. ; ; datfile: Optional name of the data fits file. ; ; errfile: Optional name of the error fits file. ; ; -- Calibration options -- ; retain: Retain negative pixel values after dark current subtraction. ; ; nocr: Don't remove cosmic rays (despiking). ; ; nohp: Don't remove hot pixels. ; ; nowp: Don't remove warm pixels. ; ; norwp: Don't remove residual warm pixels. ; ; nodp: Don't remove dusty pixels. ; ; noabs: Don't perform absolute calibration. ; ; photons: Return absolutely calibrated data in units of photons rather ; than ergs. ; ; refill: Fill missing mixels and errors with average values of the ; neighbouring valid pixels according to the method of P.Young. ; Only up to three consecutive pixels in the same column are ; filled. The method is intended to recover scientific ; information mostly lost because of the warm pixels. ; ; -- Despiking -- ; neighbours: Mark the neighbours of bad pixels as bad in EIS_DESPIKE. ; ; n_iter: ; ; _extra: ; ; -- General -- ; quiet: Don't issue xwindow messages. ; ; cleanup: Cleanup data if set. ; ; original: If set, then the original version of eis_prep will ; be used (i.e., prior to Dec-2019). ; ; OUTPUTS: ; ; If save is set, a data file and an error file are written: their names ; are "eis_l1_YYYYMMDD_HHMMSS.fits" and "eis_er_YYYYMMDD_HHMMSS.fits" by ; default. These names can be superseded by setting the datfile and ; errfile keywords.` ; ; CALLS: ; ; ; COMMON BLOCKS: ; ; ; PROCEDURE: ; ; ; RESTRICTIONS: ; ; ; MODIFICATION HISTORY: ; 22-Jun-2006 Oivind Wikstol - Version 1.0 ; 12-Feb-2007 Viggo Hansteen - Slight change in structure, added ; despiking. Still to do are error ; estimates, bad pixel marking ; and quality control (ie debugging)! ; 28-Mar-2007 Harry P Warren - New despiking keywords, new output ; keywords, fix to DC subtraction ; 02-Apr-2007 Viggo Hansteen - Force conversion to float even if only ; doing DC subtraction ; 10-Oct-2007 A. Gardini - Modified the naxis1 computation; ; set nohp=nohp and noff=noff ; 10-Oct-2007 Harry P Warren - fixed no_neighbours despiking keyword ; 19-Oct-2007 Harry P Warren - added delvarx to help with memory ; problems ; 23-Oct-2007 A. Gardini - Marking of the dusty pixels ; 20-Nov-2007 A. Gardini - CR removal performed after hot/dusty ; pixels correction ; 28-Nov-2007 A. Gardini - Added a check on the number of ; missing pixels during CR removal ; 30-Nov-2007 Viggo Hansteen - Inserted getexp method in place of ; difference between t1 and t2 ; 7-Jan-2008 A. Gardini - Check if exposure time = 0 in abs cal. ; 8-Jan-2008 A. Gardini - Modified the CAL_ keywords in l1 files. ; 10-Jan-2008 A. Gardini - Corrections in the HP calibration ; 1-Feb-2008 Mike Marsh - Added photons keyword, errors ; corrected based on photon stats. ; 2-Feb-2008 Mike Marsh - Set photon calibration status. ; 7-Feb-2008 Peter Young - Changed parameter calling eis_despike ; 12-Feb-2008 A. Gardini - Checked that the data file level is 0, ; and that ABS did not precede CR. ; 19-Feb-2008 A. Gardini - Added warm pixels calibration. ; Uploaded on 1-Apr-2004. ; 4-Apr-2008 A. Gardini - Moved CR removal before HP. Set xbox=1. ; 6-May-2008 A. Gardini - Removed the texp=0 check in abs cal. ; 16-May-2008 Peter Young - Previously eis_prep set any pixels with ; photon numbers < 1 after dc subtraction ; to be missing as the photon errors on ; such pixels did not make sense. ; For these pixels the error is now set i ; to the read noise and they are no ; longer flagged as missing. ; 19-May-2008 A. Gardini - Implemented the option above according ; to the new keyword retain. ; 9-Jul-2008 A. Gardini - Set the new values of the read out ; noise according to the quadrant. ; 20-Feb-2009 A. Gardini - Moved the error computation in DN from ; the CR removal to the ABS calibration. ; 2-Apr-2009 A. Gardini - Fixed datfile and errfile. ; 6-Apr-2009 A. Gardini - Added the refill option from P. Young. ; 28-Apr-2009 A. Gardini - Fixed an 1 pixel error in the pixel- ; wvl conversion for the CCD-A (LW) in ; the absolute calibration. ; 30-Apr-2009 A. Gardini - Set keywords to eis_cal. ; 5-May-2009 A. Gardini - Modified the refill option to deal ; with retained negative pixels. ; 6-May-2009 A. Gardini - Fixed the infinite loops in hp/wp. ; 24-Jun-2009 A. Gardini - Fixed the photon unit tag. ; 26-Jun-2009 A. Gardini - Roll back of the previous change. ; 16-Apr-2010 T. Fredvik - Added keyword hkwavecorr: if set correct ; for orbital variations of the line ; centre due to temperature variations, ; using Kamio-san's code. ; 30-Jun-2010 V. Hansteen - Added keyword ; "/correct_sensitivity" which ; should correct for CCD ; degradation as a funciton of time. ; 11-Aug-2010 T. Fredvik - Added check for missing data when ; calculating error arrays ; 19-Apr-2011 V. Hansteen - Mark ZERO DN as "missing" data ; 14-Jan-2014 V. Hansteen/H. Warren - New (correct) error ; estimate for non-absolute ; calibrated data. ; 22-Jan-2015 H. Warren - Removed stop ; 29-Jan-2015 H. Warren - corrected "correct_sensistivity" bug on line 615 ; 17-Oct-2019 H. Warren - corrected i_not_missing bug ; 17-Oct-2019 P. Young - flag 2048 DN points as missing; ; modified input behavior so that ; if input is a string then it is ; not returned as an object. ; 21-Oct-2019 P. Young - modified 2048DN check to require ; entire column have this value; ; also pixel values are set to 0 ; 04-Dec-2019 P. Young - introduced residual warm pixel ; removal (RWP); changed how CRs ; are flagged; modified order in ; which DPs, HPs, WPs and CRs and ; flagged. ; 06-Dec-2019 P. Young - fixed error when handling ; nexp_prp>1 files. This only ; affected the two most recent ; versions of this routine. ; 17-Dec-2019 P. Young - now works for slot data; ; I don't check for RWPs, ; and I use eis_despike for CRs. ; 13-Sep-2021, P. Young - added /verbose keyword. ; 13-May-2025, P. Young - corrected err for the /noabs case. ; 06-Aug-2025, P. Young - fixed array indexing problems when computing the ; error arrays (for /noabs, /phot and regular cases. ; 11-Aug-2025, P. Young - removed comment from file; fixed syntax problem ; for eis_replicate_vector. ; 14-Aug-2025, P. Young - interpolation of effective area is now done on ; the logarithm of the curve. ; 09-Mar-2026, P. Young - fixed bug in the case /nowp is set. ; ; $Id: eis_prep.pro 748 2011-04-19 10:03:15Z viggoh $ ; ;- pro eis_prep,input, $ ; default=default, dc_file=dc_file, ff_file=ff_file,$ default=default, dc_file=dc_file,$ save=save,outdir=outdir,datfile=datfile,errfile=errfile,$ retain=retain,$ refill=refill,$ nocr=nocr,nohp=nohp,nowp=nowp,nodp=nodp,noabs=noabs,$ norwp=norwp, $ photons=photons,$ correct_sensitivity=correct_sensitivity,$ ; noff=noff,$ quiet=quiet,cleanup=cleanup,$ neighbours=neighbours,n_iter=n_iter, hkwavecorr=hkwavecorr,$ _extra=ex, original=original, verbose=verbose ; ver_no='$Revision: 748 $' ; ; PRY, 5-Dec-2019 ; IF keyword_set(original) THEN BEGIN print,'% EIS_PREP: Using older version of EIS_PREP.' eis_prep_original,input, $ default=default, dc_file=dc_file,$ save=save,outdir=outdir,datfile=datfile,errfile=errfile,$ retain=retain,$ refill=refill,$ nocr=nocr,nohp=nohp,nowp=nowp,nodp=nodp,noabs=noabs,$ photons=photons,$ correct_sensitivity=correct_sensitivity,$ quiet=quiet,cleanup=cleanup,$ neighbours=neighbours,n_iter=n_iter, hkwavecorr=hkwavecorr,$ _extra=ex return ENDIF if n_params() ne 1 then begin message,'eis_prep,data [data,file], /default,/save, $',/info message,' /nohp,/nodp,/nocr,/noabs,/photons,/correct_sensitivity, $',/info message,' /nowp, $',/info ; message,' /noff, $',/info message,' dc_file=dc_file, $',/info ; message,' ff_file=ff_file, $',/info message,' outdir=outdir,datfile=datfile,errfile=errfile, $',/info message,' _extra=ex,/quiet',/info return endif if n_elements(default) eq 0 then default=0 if n_elements(save) eq 0 then save=0 if n_elements(retain) eq 0 then retain=0 if n_elements(nocr) eq 0 then nocr=0 if n_elements(nohp) eq 0 then nohp=0 if n_elements(nowp) eq 0 then nowp=0 if n_elements(norwp) eq 0 then norwp=0 if n_elements(nodp) eq 0 then nodp=0 ; if n_elements(noff) eq 0 then noff=0 if n_elements(noabs) eq 0 then noabs=0 if n_elements(photons) eq 0 then photons=0 if n_elements(correct_sensitivity) eq 0 then correct_sensitivity=0 if n_elements(quiet) eq 0 then quiet=0 IF n_elements(hkwavecorr) EQ 0 THEN hkwavecorr = 0 if default and (n_elements(dc_file) ne 0) then begin message,'WARNING: /default overrides dc_file',/info dc_file='' endif ; if default and (n_elements(dc_file) ne 0 or n_elements(ff_file) ne 0) then begin ; message,'WARNING: /default overrides dc_file and/or ff_file',/info ; dc_file='' ; ff_file='' ; endif if n_elements(dc_file) eq 0 then dc_file='' ; if n_elements(ff_file) eq 0 then ff_file='' ; open data file if necessary ; ; PRY 17-Oct-2019: introduced "file_type" below and renamed "data" to ; "input". ; - allows routine to remember if input was given as a string or ; object (see end of routine). ; case datatype(input) of 'STR': begin file_type=0b data=obj_new('eis_data',input,datasource='fits') fitslev = data-> getfitslev() if fitslev ne 0 then begin message,'eis_prep runs on level 0 files only',/info return end end 'OBJ': BEGIN file_type=1b data=input END else : begin message,'data parameter must be either eis_data object or filename',/info return end endcase fitslev=1 data-> setfits_reformat, 'ver_rf1', 'EIS_prep '+ver_no get_utc, date, /ccsds home_inst = data-> gethome_inst() data-> setfitslev,fitslev if fitslev eq 1 then begin data-> setfits_reformat, 'ver_rf1', 'EIS_prep Ver '+ver_no data-> setfits_reformat, 'date_rf1', date data-> setfits_reformat, 'orig_rf1', home_inst endif ; read in DC data: either manually pick, given in dc_file, or default set ; in eis_cal object if not default and dc_file eq '' then begin dc_dir=getenv('EIS_CAL_DATA')+path_sep()+'dc'+path_sep() title = 'Pick Dark Current file' dc_file = dialog_pickfile(path = dc_dir, filter = '*.idlsave', title = title) endif ;; read in Flat Field (FF) data (FF file must be manually picked) ; if not default and ff_file eq '' then begin ; ff_dir=getenv('EIS_CAL_DATA')+path_sep()+'ff'+path_sep() ; title = 'Pick Flat Field' ; ff_file = dialog_pickfile(path = ff_dir, filter = '*.idlsave', title = title) ; endif ;; define calibration object with dc_file and ff_file ; cal = obj_new('eis_cal', dc_file = dc_file, ff_file = ff_file, /gen_aeff) ; define calibration object with dc_file time_obs=(data->gethdr())->getdate_obs() slit_ind=data->getinfo('slit_ind') cal = obj_new('eis_cal', dc_file = dc_file, /gen_aeff,$ time_obs = time_obs, slit_ind = slit_ind) calstat='' ; ; CHECK for SATURATED DATA and for data set to ZERO ; ; PRY 18-Oct-2019: added check for 2048 DN points ; missing=-100. if not (data->getcalstat()).dc then begin for iwin=0,data->getnwin()-1 do begin wd=float(data->getvar(iwin)) (data->getcal())->setupvar,wd,iwin ; set up error array err=wd*0.0 ; i=where(wd eq 16383,ni) if ni ne 0 then err[i]=missing ; i=where(wd eq 0,ni) if ni ne 0 then err[i]=missing ; ; PRY, 21-Oct-2019 ; The following handles 2048 DN pixels. It assumes that an entire ; column has 2048 values, so I check the min and max of each column, ; and if they're both 2048 then I set the column to ; missing. Note that I replace the WD ; value with 0, which gets rid of the bright column usually seen ; for these pixels. ; nexp=data->getnexp() FOR i=0,nexp-1 DO BEGIN wdmin=min(reform(wd[*,*,i]),dim=2) wdmax=max(reform(wd[*,*,i]),dim=2) k=where(wdmin EQ 2048 AND wdmax EQ 2048,nk) IF nk GT 0 THEN BEGIN err[k,*,i]=-50. ENDIF ENDFOR k=where(err EQ -50.,nk) IF nk GT 0 THEN BEGIN print,'% EIS_PREP: window '+trim(iwin)+' has '+trim(nk)+' DN2048 pixels. Set to missing.' err[k]=missing wd[k]=0 ENDIF data->setvar,wd,iwin data->seterr,err,iwin delvarx,wd,err endfor endif ; ; DARK CURRENT subtraction ; if not (data->getcalstat()).dc then begin if retain then $ message = ['Dark current subtraction... ','Retained values less than 0.'] $ else message = ['Dark current subtraction... '] if not quiet then xmessage,message,wbase=wbase,font='helvetica' for iwin=0,data->getnwin()-1 do begin wd=float(cal->dc(data,iwin,verbose=verbose)) if not retain then wd = wd > 0 data->setvar, wd, iwin err=data->geterr(iwin) if not retain then begin i=where(wd lt 1,ni) if ni ne 0 then err[i]=missing endif data->seterr,err,iwin delvarx,wd,err endfor data->setcalstat, 'dc', 1 ; set calibration status to 'true' for dc if retain then data->setcalstat, 'retain', 1 endif else begin message,'Dark current subtraction done previously',/info endelse if (data->getcalstat()).dc then calstat='DC' else begin message,'Dark current subtraction must be performed in order ot proceed',/info return end ; ; MARK DUSTY PIXELS ; if not (data->getcalstat()).dp and not nodp then begin message = ['Dusty pixel correction...'] if not quiet then xmessage,message,wbase=wbase,font='helvetica' for iwin=0,data->getnwin()-1 do begin err=data->geterr(iwin) ; returns error window iwin dp=cal->dusty_pixels(data,iwin) for iexp=0,data->getnexp()-1 do begin i=where(dp eq 1,ni) if ni ne 0 then begin edum=err[*,*,iexp] edum[i]=missing err[*,*,iexp]=edum endif endfor data->seterr,err,iwin endfor data->setcalstat, 'dp', 1 ; set calibration status to 'true' for dp endif else begin if (data->getcalstat()).dp then message,'Dusty pixels marked previously',/info endelse if (data->getcalstat()).dp then calstat=calstat+' DP' ;************* ; PRY, 6-Nov-2019 ; Previously, eis_prep had separate sections for flagging hot pixels, ; warm pixels and cosmic rays. In each of these the flagged intensity ; pixels were interpolated from neighboring pixels. In the new code, ; I now have a single for-loop for each data window, and I then flag ; HPs, WPs, RWPs and CRs consecutively. The locations of all the ; flagged pixels get stored in SWTCH. From the SWTCH array I then do ; the interpolation. This avoids, for example, RWPs or CRs being ; interpolated into the HP or WP locations. It results in the final ; WD intensity array having less artifacts. ;************* ; nwin=data->getnwin() nexp=data->getnexp() nprp=data->getnexp_prp() ; IF NOT (data->getcalstat()).hp AND NOT nohp THEN hp_go=1 ELSE hp_go=0 IF NOT (data->getcalstat()).wp AND NOT nowp THEN wp_go=1 ELSE wp_go=0 IF NOT (data->getcalstat()).rwp and (NOT norwp) then rwp_go=1 ELSE rwp_go=0 IF NOT (data->getcalstat()).cr and not (data->getcalstat()).abs and not nocr THEN cr_go=1 ELSE cr_go=0 ; ; Don't do RWP removal for wide slit data. ; IF (slit_ind/2)*2 NE slit_ind THEN rwp_go=0 ; FOR iwin=0,nwin-1 DO BEGIN wd=data->getvar(iwin) err=data->geterr(iwin) swtch=make_array(/byte,size(wd,/dim)) IF hp_go THEN hp=cal->hot_pixels(data,iwin) IF wp_go THEN wp=cal->warm_pixels(data,iwin) FOR iexp=0,nexp-1 DO BEGIN wdum=wd[*,*,iexp] ; ; Flag hot pixels (HPs) ; IF hp_go THEN BEGIN IF hp[0] NE -1 THEN BEGIN edum=err[*,*,iexp] sdum=swtch[*,*,iexp] i=where(hp eq 1,ni) sdum[i]=1b edum[i]=missing err[*,*,iexp]=edum swtch[*,*,iexp]=sdum ENDIF ENDIF ; ; Flag warm pixels (WPs) ; IF wp_go THEN BEGIN IF wp[0] NE -1 THEN BEGIN edum=err[*,*,iexp] sdum=swtch[*,*,iexp] i=where(wp eq 1,ni) ; ni == count sdum[i]=1b edum[i]=missing err[*,*,iexp]=edum swtch[*,*,iexp]=sdum ENDIF ENDIF ; delvarx,wdum,edum,sdum ENDFOR ; wdum=wd ; ; Flag residual warm pixels (RWPs). ; IF rwp_go THEN BEGIN edum=err sdum=swtch eis_prep_residual_wps, wdum, edum, swtch=swtch_out, missing=missing, /nofill i=where(swtch_out EQ 1,ni) IF ni NE 0 THEN BEGIN sdum[i]=1 err=temporary(edum) swtch=temporary(sdum) ENDIF ENDIF ; ; Flag cosmic rays (CRs). Note that a different method is used for ; wide slit data. ; IF cr_go THEN BEGIN IF (slit_ind/2)*2 NE slit_ind THEN BEGIN edum=err sdum=swtch i=where(edum EQ missing,ni) IF ni NE 0 THEN wdum[i]=missing IF keyword_set(neighbours) THEN no_neighbours=0 ELSE no_neighbours=1 eis_despike,wdum,wdds,nprp,swtch=swtch_out,no_neighbours=no_neighbours,n_iter=n_iter,_extra=ex, $ missing=missing,/clean_exp i=where(swtch_out EQ 1,ni) IF ni NE 0 THEN BEGIN sdum[i]=1 err=temporary(edum) swtch=temporary(sdum) ENDIF ENDIF ELSE BEGIN edum=err sdum=swtch eis_prep_clean_crs, wdum, edum, swtch=swtch_out, missing=missing, /nofill i=where(swtch_out EQ 1,ni) IF ni NE 0 THEN BEGIN sdum[i]=1 err=temporary(edum) swtch=temporary(sdum) ENDIF ENDELSE ENDIF ; delvarx,wdum ; ; Interpolation of missing data. ; The four types of missing data are flagged in SWTCH. I now ; replace the intensity value at these locations using the median ; of the immediate neighbors. ; FOR iexp=0,nexp-1 DO BEGIN wdum=reform(wd[*,*,iexp]) sdum=reform(swtch[*,*,iexp]) index=where(sdum EQ 1,ni) IF ni NE 0 THEN BEGIN edum=reform(err[*,*,iexp]) new_img=eis_prep_med_fill(wdum,index=index,/quiet,missing=missing) wd[*,*,iexp]=new_img ENDIF ENDFOR data->seterr,err,iwin data->setvar,wd,iwin delvarx,wd,err,swtch ENDFOR ; ; Now update calibration status for each of HP, WP, RWP and CR. ; IF hp_go THEN BEGIN data->setcalstat, 'hp', 1 calstat=calstat+' HP' ENDIF ; IF wp_go THEN BEGIN data->setcalstat, 'wp', 1 calstat=calstat+' WP' ENDIF ; IF rwp_go THEN BEGIN data->setcalstat, 'rwp', 1 calstat=calstat+' RWP' ENDIF ; IF cr_go THEN BEGIN data->setcalstat, 'cr', 1 calstat=calstat+' CR' ENDIF ; ; ABSOLUTE CALIBRATION and ERROR estimation ; if not (data->getcalstat()).abs and not noabs then begin message = ['Absolute calibration...'] if not quiet then xmessage,message,wbase=wbase,font='helvetica' exp=data->getexp() t0 = where(exp eq 0.,ct0) if (ct0 gt 0) then begin message,'Some exposure times are equal to 0. Set to 1 to complete calibration',/info exp[t0] = 1. endif case data->getinfo('slit_ind') of 0: begin slitsize=1.0*1.0 if not photons then begin data->setunit,['erg/cm!e2!N/s/sr/'+string("305B),'',''] endif else begin data->setunit,['Photon-Events'] endelse end 2: begin slitsize=1.0*2.0 if not photons then begin data->setunit,['erg/cm!e2!N/s/sr/'+string("305B),'',''] endif else begin data->setunit,['Photon-Events'] endelse end else: begin slitsize=-1.0 if not photons then begin data->setunit,['erg/cm!e2!N/s/sr','',''] endif else begin data->setunit,['Photon-Events'] endelse message,'This is slot data - using slitsize of 1.0 and dispersion of 1.0 to compute absolute calibration',/info end endcase for iwin=0,data->getnwin()-1 do begin data->getwin,iwin,dn,pos err=data->geterr(iwin) i_not_miss=where(err ne missing, count) if count eq 0 then begin intensity = err data->setvar, intensity, iwin print, ' ! no good data, skipping this window ' + trim(iwin) continue endif ;; --- read out noise definition quad=cal->quadrant(pos) case quad of 0: begin read_noise = 2.26 ;; in DN end 1: begin read_noise = 2.29 ;; in DN end 2: begin read_noise = 2.37 ;; in DN end 3: begin read_noise = 2.24 ;; in DN end endcase if pos[0] gt 2147 then begin detector='A' pos[0]=pos[0]-2148 lambda=((cal->getlambda()).scale_a)[pos[0]:pos[0]+pos[1]-1] endif else begin detector='B' lambda=((cal->getlambda()).scale_b)[pos[0]:pos[0]+pos[1]-1] endelse aeff = cal->getaeff(detector) ; ; PRY, 14-Aug-2025 ; Effective area interpolation now done on log of data. ; ; effarea=interpol(aeff.eff,aeff.lambda,lambda,/spline) effarea=10.^interpol(alog10(aeff.eff),aeff.lambda,lambda,/spline) sub=where(lambda lt aeff.lambda[0] or lambda gt aeff.lambda[n_elements(aeff.lambda)-1]) if sub[0] ne -1 then effarea[sub]=1.e-20 ; intensity in photons, for error calculation intensity_ph = cal->dn_to_photon(lambda,float(dn)) data->setcalstat,'sens',0 if not photons then begin days_since_launch=0.0 if correct_sensitivity then days_since_launch=data->days_since_launch() if days_since_launch ne 0.0 then begin if correct_sensitivity eq 1 then begin action=' Correcting ' (data->getcal())->setcorrected_sensitivity,1 data->setcalstat, 'sens', 1 endif ; else begin ; action=' Uncorrecting ' ; days_since_launch=-1*days_since_launch ; (data->getcal())->setcorrected_sensitivity,0 ; endelse message,string(days_since_launch,format='(f7.1)')+' days since lauch.' $ +action+'sensitivity with tau = ' $ +string((data->getcal())->gettau_sensitivity(),format='(f7.1)')+' days.',/info endif ; print,(data->getcal())->sensitivity_loss(days_since_launch) intensity = cal->photon_to_ergs(lambda*effarea,intensity_ph, $ exp=exp,slitsize=slitsize,days_since_launch=days_since_launch) endif else intensity=intensity_ph ;leave intensity in photons if photons keyword set ;; --- add in the read noise rnoise_dn = fltarr(size(dn,/dimensions))+read_noise ;sets read noise array in DN rnoise_ph = cal->dn_to_photon(lambda,float(rnoise_dn)) ;read noise array in photons ;; PRY 12-Feb-08 - modified line below by adding "> 0." ;; PRY 6-Aug-2025 - I inserted [i_not_miss] for rnoise_ph since it has the same ; size as err and intensity_ph err[i_not_miss] = sqrt(intensity_ph[i_not_miss] > 0. +rnoise_ph[i_NOT_miss]^2) ;photon variance = photon counts data->setvar, intensity,iwin ;; Changed method of calculating error using the photon_to_ergs method ;; Old method follows below (commented out). PRY if not photons then begin i_miss=where(err EQ missing) err = cal->photon_to_ergs(lambda*effarea,err, $ exp=exp,slitsize=slitsize,days_since_launch=days_since_launch) err[i_miss]=missing endif ; err[i_not_miss]=intensity[i_not_miss]*err[i_not_miss]/intensity_ph[i_not_miss] ;multiply by fractional error in photons data->seterr,err,iwin delvarx,err,intensity,i_not_miss ;print," MEMORY = "+trim(double(memory(/current))/1.0E+6) endfor data->setcalstat, 'abs', 1 ; set calibration status to 'true' for ff if photons then data->setcalstat, 'phot', 1 ;Mike Marsh 2/2/08 set photon calibration status to 'true' if keyword set endif else begin if (data->getcalstat()).abs then begin message,'Absolute calibration done previously',/info if correct_sensitivity eq -1 and (data->getcal())->getcorrected_sensitivity() then begin message,'Uncorrecting for sensitivity with tau = ' $ +string((data->getcal())->gettau_sensitivity(),format='(f7.1)')+' days.',/info for iwin=0,data->getnwin()-1 do begin data->getwin,iwin,intensity,pos err=data->geterr(iwin) ; days_since_launch=data->days_since_launch() sens_loss=(data->getcal())->sensitivity_loss(abs(days_since_launch)) intensity=intensity*sens_loss err=err*sens_loss data->setvar,intensity,iwin data->seterr,err,iwin ; (data->getcal())->setcorrected_sensitivity,0 data->setcalstat, 'sens', 0 ; delvarx,err,intensity endfor endif endif else begin ; error computation in DN for iwin=0,data->getnwin()-1 do begin data->getwin,iwin,dn,pos err=data->geterr(iwin) i_not_miss=where(err ne missing, nmissing) ;; --- read out noise definition quad=cal->quadrant(pos) case quad of 0: begin read_noise = 2.26 ;; in DN end 1: begin read_noise = 2.29 ;; in DN end 2: begin read_noise = 2.37 ;; in DN end 3: begin read_noise = 2.24 ;; in DN end endcase ;; [Start HPW 13-JAN-2016] ;; get wavelength if pos[0] gt 2147 then begin detector='A' pos[0]=pos[0]-2148 lambda=((cal->getlambda()).scale_a)[pos[0]:pos[0]+pos[1]-1] endif else begin detector='B' lambda=((cal->getlambda()).scale_b)[pos[0]:pos[0]+pos[1]-1] endelse ;; sets read noise array in DN rnoise_dn = fltarr(size(dn,/dimensions))+read_noise IF nmissing GT 0 THEN BEGIN photon_to_dn=cal->photon_to_dn(lambda,1) ; ; PRY, 6-Aug-2025: convert the 1D array photon_to_dn to an array ; that matches DN and RNOISE_DN. ; dn_siz=size(dn,/dim) photon_to_dn_arr=eis_replicate_vector(photon_to_dn,dn_siz[1:*]) ; ; PRY, 13-May-2025: The expression below was not quite right. ; PRY, 6-Aug-2025: rnoise_dn needed [i_not_miss] indices. ; err[i_not_miss] = sqrt(dn[i_not_miss] > 0. + rnoise_dn[i_not_miss]^2/photon_to_dn_arr[i_not_miss])*sqrt(photon_to_dn_arr[i_not_miss]) ; err[i_not_miss] = sqrt(dn[i_not_miss] > 0. + rnoise_dn^2/photon_to_dn[i_NOT_miss])*sqrt(photon_to_dn[i_not_miss]) ; err[i_not_miss] = sqrt(dn[i_not_miss] > 0. + rnoise_dn^2)*sqrt((cal->photon_to_dn(lambda,1))[i_not_miss]) endif ELSE err[*] = -100 data->seterr,err,iwin delvarx,err,dn,i_not_miss endfor ;; [Stop HPW 13-JAN-2016] endelse endelse if (data->getcalstat()).abs then calstat=calstat+' ABS' if (data->getcalstat()).sens then calstat=calstat+' SENS' if (data->getcalstat()).phot then calstat=calstat+' PHOT' if (data->getcalstat()).retain then calstat=calstat+' RETAIN' ; ; UPDATE data header ; hdr = data->gethdr() main_hdr = hdr->getmain_hdr() bte_hdr = hdr-> getbte_hdr() fits_reformat=data->getfits_reformat() fxaddpar, main_hdr, 'DATA_LEV', (data-> getfitslev()), $ 'FITS Level (0, 1, 2)' fxaddpar, main_hdr, 'DATE_RF1', fits_reformat.date_rf1, $ 'Date and time of Level 1 reformat' fxaddpar, main_hdr, 'ORIG_RF1', fits_reformat.orig_rf1, $ 'Institution where Level 1 reformat was done' fxaddpar, main_hdr, 'VER_RF1', fits_reformat.ver_rf1, $ 'Fits Level 1 reformatter version no.' fxaddpar, bte_hdr, 'DATA_LEV', (data-> getfitslev()), $ 'FITS Level (0, 1, 2)' fxaddpar, bte_hdr, 'DATE_RF1', fits_reformat.date_rf1, $ 'Date and time of Level 1 reformat' fxaddpar, bte_hdr, 'ORIG_RF1', fits_reformat.orig_rf1, $ 'Institution where Level 1 reformat was done' fxaddpar, bte_hdr, 'VER_RF1', fits_reformat.ver_rf1, $ 'Fits Level 1 reformatter version no.' calflag=data->getcalstat() fxaddpar,main_hdr,'CALSTAT',calstat,'Calibration status', after = 'CMIRR ' fxaddpar,bte_hdr,'CALSTAT', calstat,'Calibration status', after = 'CMIRR ' fxaddpar,main_hdr,'CAL_DC',calflag.dc,after='CALSTAT' fxaddpar,bte_hdr,'CAL_DC',calflag.dc,after='CALSTAT' fxaddpar,main_hdr,'CAL_HP',calflag.hp,after='CAL_DC' fxaddpar,bte_hdr,'CAL_HP',calflag.hp,after='CAL_DC' fxaddpar,main_hdr,'CAL_DP',calflag.dp,after='CAL_HP' fxaddpar,bte_hdr,'CAL_DP',calflag.dp,after='CAL_HP' fxaddpar,main_hdr,'CAL_WP',calflag.dp,after='CAL_DP' fxaddpar,bte_hdr,'CAL_WP',calflag.dp,after='CAL_DP' ; fxaddpar,main_hdr,'CAL_FF',calflag.ff ; fxaddpar,bte_hdr,'CAL_FF',calflag.ff fxaddpar,main_hdr,'CAL_RWP',calflag.rwp,after='CAL_WP' fxaddpar,bte_hdr,'CAL_RWP',calflag.rwp,after='CAL_WP' fxaddpar,main_hdr,'CAL_CR',calflag.cr,after='CAL_RWP' fxaddpar,bte_hdr,'CAL_CR',calflag.cr,after='CAL_RWP' ; fxaddpar,main_hdr,'CAL_WVL',calflag.wvl ; fxaddpar,bte_hdr,'CAL_WVL',calflag.wvl fxaddpar,main_hdr,'CAL_ABS',calflag.abs,after='CAL_CR' fxaddpar,bte_hdr,'CAL_ABS',calflag.abs,after='CAL_CR' fxaddpar,main_hdr,'CAL_PHOT',calflag.phot,after='CAL_ABS' fxaddpar,bte_hdr,'CAL_PHOT',calflag.phot,after='CAL_ABS' fxaddpar,main_hdr,'CAL_SENS',calflag.sens,after='CAL_PHOT' fxaddpar,bte_hdr,'CAL_SENS',calflag.sens,after='CAL_PHOT' fxaddpar,main_hdr,'CAL_RETA',calflag.retain,after='CAL_SENS' fxaddpar,bte_hdr,'CAL_RETA',calflag.retain,after='CAL_SENS' sxdelpar, main_hdr, ['CAL_FF','CAL_WVL'] sxdelpar, bte_hdr, ['CAL_FF','CAL_WVL'] ; if calflag.sens then begin fxaddpar,main_hdr,'TAU_SENS',(data->gettau_sensitivity()),$ 'Sensitivity decrease e-folding time (days)',after='CAL_RETA' fxaddpar,bte_hdr,'TAU_SENS',(data->gettau_sensitivity()),$ 'Sensitivity decrease e-folding time (days)',after='CAL_RETA' endif ; max and min values in line windows have changed. ; TDMIN and TDMAX must be set again nwin = data-> getnwin() naxis1=0l for iwin=0, nwin-1 do begin data-> getwin, iwin, wd, pos wdmin = min(wd) ; minimum value in line window wdmax = max(wd) ; maximum value in line window param1 = 'TDMIN'+strtrim(string(iwin+1), 2) param2 = 'TDMAX'+strtrim(string(iwin+1), 2) fxaddpar, bte_hdr, param1, wdmin fxaddpar, bte_hdr, param2, wdmax param1 = 'TFORM'+strtrim(string(iwin+1), 2) winsize=ulong((data->getxw())[iwin])*ulong((data->getyw())[iwin]) winstring=strtrim(string(winsize),2) case datatype(wd) of 'FLO': begin naxis1=winsize*4l+naxis1 fxaddpar,bte_hdr,param1,winstring+'E','Real*4 (floating point)' end 'INT': begin naxis1=winsize*2l+naxis1 fxaddpar,bte_hdr,param1,winstring+'I','Integer*2 (short integer)' end else : message,'Unexpected datatype for wd, all bets are off!',/info endcase endfor ; naxis1=naxis1+73l ; NB this is a HACK!!!! should count number of bytes added in auxilary data naxis1=naxis1+65l if ptr_valid((data->getaux_data()).mhc_hz_t10) then naxis1=naxis1+8l if ptr_valid((data->getaux_data()).v_sat) then naxis1=naxis1+12l ; by reading the bte_hdr as it exists. fxaddpar,main_hdr,'NAXIS1',naxis1,'Number of bytes per row' fxaddpar,bte_hdr,'NAXIS1',naxis1,'Number of bytes per row' hdr-> setmain_hdr, main_hdr hdr-> setbte_hdr, bte_hdr ; ; REFILLING missing pixels according to the method discussed by PRY in ; http://msslxr.mssl.ucl.ac.uk:8080/eiswiki/Wiki.jsp?page=WarmPixelStudy ; if keyword_set(refill) then begin message = ['Refilling missing pixel',$ 'and errors: up to three',$ 'pixels in a column are filled.'] if not quiet then xmessage,message,wbase=wbase,font='helvetica' for iwin=0,data->getnwin()-1 do begin data->getwin,iwin,wd,pos err=data->geterr(iwin) ; computation of wvl and ea if pos[0] gt 2147 then begin detector='A' pos[0]=pos[0]-2148 lambda=((cal->getlambda()).scale_a)[pos[0]:pos[0]+pos[1]-1] endif else begin detector='B' lambda=((cal->getlambda()).scale_b)[pos[0]:pos[0]+pos[1]-1] endelse aeff = cal->getaeff(detector) effarea=interpol(aeff.eff,aeff.lambda,lambda,/spline) ; the routine spl_interp returns slightly different results. y2 = spl_init(aeff.lambda,aeff.eff) effarea=spl_interp(aeff.lambda,aeff.eff,y2,lambda) sub=where(lambda lt aeff.lambda[0] or lambda gt aeff.lambda[n_elements(aeff.lambda)-1]) if sub[0] ne -1 then effarea[sub]=1.e-20 ;----- sz=size(wd) nl=sz[1] ny=sz[2] nt=sz[3] wvl=fltarr(nl,ny,nt) ea=fltarr(nl,ny,nt) for k=0,nl-1 do begin wvl[k,*,*]=lambda[k] ea[k,*,*]=effarea[k] endfor ;----- i=where(err eq missing,nmiss) if nmiss gt 0 then begin frac=[1.0,1.2,1.2,1.3,1.3] ; chck=bytarr(nl,ny,nt) miss_map=bytarr(nl,ny,nt) miss_map(i)=1b tofill_map=miss_map x=fltarr(nl,ny,nt) y=fltarr(nl,ny,nt) x=wd*ea y=err^2*ea^2*wvl j=where(err ne missing and wd gt 0.,nj) f=poly_fit(x[j],y[j],1) ; wdaux=fltarr(nl,ny,nt) wdout=fltarr(nl,ny,nt) wdout=wd errout=err ; chck[*]=1b wdaux[*,0:ny-2,*]=wd[*,1:ny-1,*] chck[*,0:ny-2,*]=miss_map[*,1:ny-1,*] k=where(chck[i] eq 0b,nk) frac_val=frac[4] if nk ne 0 then begin wdout[i[k]]=wdaux[i[k]] ; errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac_val errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $ (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val tofill_map[i[k]]=0b endif ; chck[*]=1b wdaux[*,1:ny-1,*]=wd[*,0:ny-2,*] chck[*,1:ny-1,*]=miss_map[*,0:ny-2,*] k=where(chck[i] eq 0b,nk) frac_val=frac[4] if nk ne 0 then begin wdout[i[k]]=wdaux[i[k]] ; errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[4] errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $ (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val tofill_map[i[k]]=0b endif ; chck[*]=1b wdaux[*,2:ny-3,*]=0.5*(wd[*,0:ny-5,*]+wd[*,4:ny-1,*]) chck[*,2:ny-3,*]=miss_map[*,0:ny-5,*]+miss_map[*,4:ny-1,*] k=where(chck[i] eq 0b,nk) frac_val=frac[3] if nk ne 0 then begin wdout[i[k]]=wdaux[i[k]] ; errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[3] errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $ (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val tofill_map[i[k]]=0b endif ; chck[*]=1b wdaux[*,3:ny-2,*]=(wd[*,0:ny-5,*]+3.5*wd[*,4:ny-1,*])/4.5 chck[*,3:ny-2,*]=miss_map[*,0:ny-5,*]+miss_map[*,4:ny-1,*] k=where(chck[i] eq 0b,nk) frac_val=frac[2] if nk ne 0 then begin wdout[i[k]]=wdaux[i[k]] ; errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[2] errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $ (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val tofill_map[i[k]]=0b endif ; chck[*]=1b wdaux[*,1:ny-4,*]=(3.5*wd[*,0:ny-5,*]+wd[*,4:ny-1,*])/4.5 chck[*,1:ny-4,*]=miss_map[*,0:ny-5,*]+miss_map[*,4:ny-1,*] k=where(chck[i] eq 0b,nk) frac_val=frac[2] if nk ne 0 then begin wdout[i[k]]=wdaux[i[k]] ; errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[2] errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $ (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val tofill_map[i[k]]=0b endif ; chck[*]=1b wdaux[*,2:ny-2,*]=(wd[*,0:ny-4,*]+2.*wd[*,3:ny-1,*])/3. chck[*,2:ny-2,*]=miss_map[*,0:ny-4,*]+miss_map[*,3:ny-1,*] k=where(chck[i] eq 0b,nk) frac_val=frac[1] if nk ne 0 then begin wdout[i[k]]=wdaux[i[k]] ; errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[1] errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $ (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val tofill_map[i[k]]=0b endif ; chck[*]=1b wdaux[*,1:ny-3,*]=(2.*wd[*,0:ny-4,*]+wd[*,3:ny-1,*])/3. chck[*,1:ny-3,*]=miss_map[*,0:ny-4,*]+miss_map[*,3:ny-1,*] k=where(chck[i] eq 0b,nk) frac_val=frac[1] if nk ne 0 then begin wdout[i[k]]=wdaux[i[k]] ; errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[1] errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $ (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val tofill_map[i[k]]=0b endif ; chck[*]=1b wdaux[*,1:ny-2,*]=0.5*(wd[*,0:ny-3,*]+wd[*,2:ny-1,*]) chck[*,1:ny-2,*]=miss_map[*,0:ny-3,*]+miss_map[*,2:ny-1,*] k=where(chck[i] eq 0b,nk) frac_val=frac[0] if nk ne 0 then begin wdout[i[k]]=wdaux[i[k]] ; errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0]) errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $ (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val tofill_map[i[k]]=0b endif ; For all remaining pixels, set intensity and error to be missing. k=where(tofill_map eq 1b,nk) if nk ne 0 then begin wdout[k]=missing errout[k]=missing endif ; data->setvar,wdout,iwin data->seterr,errout,iwin end endfor endif IF hkwavecorr THEN data->sethkpixcorr ; xkill,wbase obj_destroy,cal ; ; SAVE data to fits file ; if keyword_set(save) then begin file = data->getfilename() break_file,file,disk,dir,name,ext if not keyword_set(datfile) then $ datfile = str_replace(name,'l0','l1')+'.fits' if not keyword_set(errfile) then $ errfile = str_replace(name,'l0','er')+'.fits' if keyword_set(outdir) then begin datfile = concat_dir(outdir,datfile) errfile = concat_dir(outdir,errfile) endif message,' saving data to '+datfile,/info data->save,file=datfile message,' saving error to '+errfile,/info data->saveerr,file=errfile file = datfile endif ; ; PRY, 17-Oct-2019 ; Depending if the input FILE was a string (file_type=0) or an ; object (file_type=1), then I do a clean-up of DATA. ; IF file_type EQ 1 THEN file=temporary(data) ELSE obj_destroy,data if n_elements(cleanup) ne 0 then obj_destroy,data ; return end