;+ ; NAME: eisobs__define ; ; ; ; PURPOSE: Definition of the eisobs (EIS observation) class, a subclass of ; hinobs (Hinode observation) which is a subclass of anyobs (any ; observation). The class definition contains methods that differs ; from those of hinobs and ; anyobs. Since the data tag of eis is a pointer to an object ; (eis_data) instead of a pointer to a data cube (this is the case ; for e.g. xrtobs and sotobs), there are significant differences ; between eisobs methods and the methods of hinobs/anyobs and its ; subclasses. ; ; ; ; (sdc_util) ; | ; anyobs ; / \ ; hinobs irisobs ; / | \ ; xrtobs sotobs EISOBS ; / \ ; sp1dobs spl2obs ; ; ; ; METHODS: See hinobs__define and anyobs__define ; ; ; ; SPECIAL CALLS: hinobs__define.pro, gt_fkw.pro, prune_filename.pro, ; eis_slotwin_oshadow.pro, gt_peaks.pro, gt_closest_index.pro, ; gt_fkw.pro, gt_snratio.pro, eis_linelist_PYoung.txt, ; eis_linelist_MSSL.txt ; ; ; ; ; WRITTEN: Terje Fredvik, UiO, December 2007 ; ; ; ; ; MODIFICATION HISTORY: ; ; Tue Jun 7 14:24:53 2011, Terje Fredvik ; 24. Aug 2010 Terje Fredvik: Removed old modification history. ; 24. Aug 2010 Terje Fredvik: # ::_extract_image_slit_cfit; ; made Kamio-san's house keeping ; temperature wavelength correction ; method default. # Removed several ; obsolete methods (the code has been moved ; to eis_data__define). # ; ::_gt_warren_lamcorr; fixed bug in slit ; tilt for the long wavelength detector when ; using eis_wave_corr_old (i.e. when keyword ; kamio is set to 0). # ::get_expmissing ; fixed bug occuring when there are no ; missing data. ; 07. Sep 2010 Terje Fredvik: # Added method :::_auto_adef() (code used ; to be in an external function). # Call ; SSW function trim() instead of home made ; flt2str(). # Added keyword medianlam to ; ::_extract_cfit_gausslam and added code to ; deal with user defined gaussian ; components. # Renamed variables and ; routines: autogauss_lam -> gausslam ; 15. Oct 2010 Terje Fredvik: # Get linelists from easy-to-edit txt ; files instead of IDL save files. The line ; lists are stored in commonblocks. ; 19. Oct 2010 Terje Fredvik: ::_extract_image_slit_cfit: fixed indexing ; bug ; 21. Jan 2010 Terje Fredvik: # call eis_slotwin instead of ; eis_slotwin_sdc, eis_wave_corr_old instead ; of eis_wave_corr_sdc (due to renaming) ; 31. Jan 2011 Terje Fredvik: # Do not set version (SVN revision is set ; in hinobs__define) ; 24. May 2011 Terje Fredvik: Yet another renaming: eis_slotwin -> ; eis_slotwin_oshadow ; 07. Jun 2011 Terje Fredvik: Fixed path to line list files ; 17. Jun 2011 Terje Fredvik: Search for line list files if not ; found in default directory ; 25. Jul 2012 Terje Fredvik: Deleted _plot_helio, code moved to hinobs::plot_helio ; 11. Apr 2013 Terje Fredvik: added info about new superclass anyobs ; 17. Oct 2013 Terje Fredvik: ::_extract_image_slit_cfit; got rid of the ; nexp_prp test, rebin works anyhow. Minor ; changes in variable names and syntax. ; ::_extract_cfit_im; removed unused ; input parameter lam. ; 11. Nov 2013 Terje Fredvik: ::_extract_image_slit_cfit: if the eis_data ; object fails in finding a wavelength ; correction, do not crash, but keep on ; working with non-corrected data. ; 22. Oct 2014 Terje Fredvik: Call eis_data::getxcen and getycen to ; allow for variable x/ycen for different ; detector windows. Removed the old and ; inaccurate method ; ::_correct_detector_offset. ; 27. Jul 2022 Terje Fredvik: call cfit_block with keyword quiet set ; 04. May 2023 Terje Fredvik: In ::xcen and ::ycen, call new function ::_eis_aia_offsets, ; which calls stand-alone function ; eis_aia_offsets to improve ; pointing accuracy. ; ;- ; FUNCTION eisobs::init, file, destdir=destdir, hk_cmn_dir=hk_cmn_dir, $ hk_status_dir=hk_status_dir, keepprep=keepprep, $ keeperror=keeperror, refill=refill, _ref_extra=extra o = self IF self->hinobs::init(file, destdir=destdir, _extra=extra) EQ 0 THEN return, 0 not_set = -1 o.slit = not_set o.raster = not_set o.widewindows = not_set o.subtitle = 'MISSING' o.firstdim = 'MISSING' o.singleexp = not_set o.missing = -1e3 ; To be used in gaussian fitting default, refill, 0 o.refill = refill o.hk_cmn_dir = (keyword_set(hk_cmn_dir)) ? hk_cmn_dir : '' o.hk_status_dir = (keyword_set(hk_status_dir)) ? hk_status_dir : '' default, keepprep,1 o.rmprepfiles = ~keepprep default, keeperror, 0 o.rmerrorfile = ~keeperror o.ii = not_set return,1 END PRO eisobs::cleanup o = self IF ptr_valid(o.data) THEN BEGIN obj_destroy,*o.data ptr_free, o.expgood ptr_free, o.expmissing ;; If data is prepped a level 1 file and an error file has been created ;; and must be deleted if keyword keep_prep is not set. IF o->get(/prep) AND o->get(/rmprepfiles) THEN BEGIN file_delete, file_search(o->find_prepfile('level1')+'*'), $ o->find_prepfile('lamcorr'), o->find_prepfile('fit'), $ /allow_nonexistent FOR i=0,o->get(/nicons)-1 DO BEGIN file_delete, o->find_prepfile('ana',i=i), $ o->find_prepfile('adef',i=i), /allow_nonexistent ENDFOR IF o->get(/rmerrorfile) THEN file_delete, o->get(/errorfile),$ /allow_nonexistent ENDIF ENDIF o->hinobs::cleanup print,"Cleaning up (eisobs)" END ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;; <_ensure_instrument_specific_properties & help methods> ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ PRO eisobs::_ensure_instrument_specific_properties o = self o->private o->_ensure_slit o->_ensure_raster o->_ensure_widewindows o->_ensure_singleexp o->_ensure_firstdim o->_ensure_expinfo END PRO eisobs::_ensure_slit o = self o->private IF ~o->is_set(/slit) THEN BEGIN slit_id = strtrim(gt_fkw(*(o->get(/header)),'SLIT_ID')) o.slit = (slit_id EQ '1"' OR slit_id EQ '2"') ? 1:0 END END PRO eisobs::_ensure_raster o = self o->private IF ~o->is_set(/raster) THEN BEGIN nraster = gt_fkw(*(o->get(/header)),'NRASTER') o.raster = (nraster GT 1) ? 1 : 0 END END PRO eisobs::_ensure_widewindows;_when_data o = self self->private IF ~o->is_set(/widewindows) THEN BEGIN tdetxw1 = gt_fkw(*(o->get(/header)),'TDETXW1') o.widewindows = (tdetxw1 GT 360) ? 1 : 0 ENDIF END PRO eisobs::_ensure_firstdim o = self o->private IF ~o->is_set(/firstdim) THEN BEGIN slit = o->get(/slit) ; eq 1 if 1" or 2" slit raster = o->get(/raster) ; eq 1 if number of unique raster positions gt 1 ;; Calculate the aspect ratio (FOVY/FOVX) for the observation. If the ;; aspect ratio is high we will use the wavelength as the first ;; dimension. For sit-and-stare observations we use NEXP instead of ;; FOVX. Also, for SAS we may use non-quadratic pixels, so we can stretch ;; pixels in the x direction, i.e. we can tolerate fewer x pixels than ;; for solar X VS solar Y images. IF slit THEN BEGIN header = *(o->get(/header)) fovx = gt_fkw(header, 'FOVX') fovy = gt_fkw(header, 'FOVY') IF raster THEN BEGIN firstdim = 'Solar X (arc sec)' IF fovx LT 1 THEN BEGIN ;; Something's wrong! (this rarely happens!) Let's estimate FOVX: slitsz = fix(strmid(gt_fkw(header, 'SLIT_ID'),0,1)) nraster = gt_fkw(header, 'NRASTER') fovx = slitsz*nraster ;; Assuming stepsize of one slit width... ENDIF aspectratio = fovy/fovx ENDIF ELSE BEGIN firstdim = 'Time (min)' ;; Sit-and-stare; we may have non-quadratic pixels, so we can stretch ;; pixels in the x direction, i.e. we can tolerate a much greater ;; aspect ratio than for x VS y images. "Much" is here defined to be 3. nexp = gt_fkw(header, 'NEXP') aspectratio = fovy/nexp/3 ENDELSE ENDIF ELSE BEGIN ;; Slot observations (40" or 266") firstdim = (o->get(/widewindows)) ? 'Wavelength (Å)': 'Solar X (arc sec)' aspectratio = -1 ENDELSE ;; aspectratio gt 35 corresponds to e.g. a FOVX gt 15" if the 1"/2" X ;; 512" slit is used IF aspectratio GT 35 THEN firstdim = 'Wavelength (Å)' o.firstdim = firstdim END END PRO eisobs::_ensure_expinfo o = self o->private eis_data = *(o->get(/data)) expmissing = where(eis_data->check_ti() EQ 0, nexpmissing, $ complement=expgood) self.expgood = ptr_new(expgood) self.expmissing = ptr_new(expmissing) END PRO eisobs::_ensure_singleexp o = self o->private IF ~o->is_set(/singleexp) THEN BEGIN nexp = gt_fkw(*(o->get(/header)), 'NEXP') o.singleexp = (nexp EQ 1) ? 1 : 0 ENDIF END ;; ----------------------------------------------------- ;; ;; ----------------------------------------------------- ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;; private self.nicons = gt_fkw(*self.header,'NWIN') END PRO eisobs::_ensure_title o = self o->private nicons = o->get(/nicons) o.title = ptr_new(strarr(nicons)) header = *(o->get(/header)) FOR i=0,nicons-1 DO (*o.title)[i] = $ gt_fkw(header,'ttype'+trim(i+1)) END PRO eisobs::_ensure_xtitle o = self o->private IF ~o->is_set(/xtitle) THEN o.xtitle = o->get(/firstdim) END PRO eisobs::_ensure_subtitle o = self o->private IF ~o->is_set(/subtitle) THEN BEGIN header = *(o->get(/header)) nexp_prp = gt_fkw(header,'NEXP_PRP') exp_txt = (nexp_prp GT 1) ? trim(nexp_prp)+' exp. summed' : '' nraster = gt_fkw(header,'NRASTER') pos_add = (nraster GT 1 AND o->get(/widewindows)) ? ' ('+trim(nraster)+ $ ' unique slit positions)' : '/'+trim(nraster) pos_txt = (nraster GT 1 AND strmid(o->get(/xtitle),0,7) NE 'Solar X') $ ? 'slit pos 0'+ pos_add : '' comma = (exp_txt NE '' AND pos_txt NE '') ? ', ' : '' o.subtitle = exp_txt + comma + pos_txt ENDIF END ;; <_ensure_image_scale and help methods> PRO eisobs::_ensure_image_scale o = self o->private IF ~o->is_set(/image_scale) THEN BEGIN xtitle = strmid(o->get(/xtitle),0,4) CASE xtitle OF 'Sola': image_scale = ptr_new(o->_image_scale_x_y()) ; in hinobs__define 'Time': image_scale = ptr_new(o->_image_scale_time_y()) 'Wave': image_scale = ptr_new(o->_image_scale_wavelength_y()) ELSE: message, 'Unknown xtitle' ENDCASE o.image_scale = image_scale END END FUNCTION eisobs::_image_scale_x self->private header = *(self->get(/header)) return, gt_fkw(header,'fovx')/gt_fkw(header,'nraster') END FUNCTION eisobs::_image_scale_time_y o = self o->private header = *(o->get(/header)) DATE_END = gt_fkw(header,'DATE_END') DATE_OBS = gt_fkw(header,'DATE_OBS') nexp = gt_fkw(header, 'NEXP') nexp_prp = gt_fkw(header, 'NEXP_PRP') IF strmid(DATE_END,10,1) EQ 'T' AND strmid(DATE_OBS,10,1) EQ 'T' THEN $ scale_time = ((anytim2tai(DATE_END) - $ anytim2tai(DATE_OBS))/NEXP*NEXP_PRP)/60. $ ELSE BEGIN scale_time = 10/60. print,'DATE_END and/or DATE_OBS not defined. X-image_scale of sit-and-stare observation will not be correct. Time between exposures assumed to be 10 sec.' ENDELSE return, [scale_time, o->_image_scale_y()] END FUNCTION eisobs::_image_scale_wavelength_y o = self o->private header = *(o->get(/header)) nicons = o->get(/nicons) scale_y = o->_image_scale_y() scale_wavelength_y = fltarr(nicons,2) scale_wavelength_y[*, *] = 1 FOR i=0,nicons-1 DO BEGIN n = trim(i+1) wavemin = gt_fkw(header,'TWMIN'+n) wavemax = gt_fkw(header,'TWMAX'+n) nwavepix = gt_fkw(header,'TDETXW'+n) scale_wavelength_y[i,*] = [(wavemax-wavemin)/(nwavepix-1), scale_y] ENDFOR return, scale_wavelength_y END ;; ;; <_ensure_image_origin and help mehtods> PRO eisobs::_ensure_image_origin o = self o->private IF ~o->is_set(/image_origin) THEN BEGIN xtitle = strmid(o->get(/xtitle),0,4) CASE xtitle OF 'Sola': image_origin = ptr_new(o->_image_origin_x_y()) 'Time': image_origin = ptr_new(o->_image_origin_time_y()) 'Wave': image_origin = ptr_new(o->_image_origin_wavelength_y()) ELSE: message, 'Unknown xtitle' ENDCASE o.image_origin = image_origin END END FUNCTION eisobs::_image_origin_x_y o = self o->private nicons = o->get(/nicons) xy = findgen(nicons,2) ; Array with the x and y coordinates of the image origins ;; Get the x and y coordinates from the fits header xy[*,0] = o->_image_origin_x() xy[*,1] = o->_image_origin_y() return, xy END FUNCTION eisobs::_image_origin_x self->private xcen = self->_xcen() fovx = gt_fkw(*(self->get(/header)),'FOVX') return, xcen-fovx/2. END FUNCTION eisobs::_eis_aia_offsets, x=x, y=y ;; Correct for long-term pointing variations by alligning with AIA images ;; http://solarb.mssl.ucl.ac.uk/eiswiki/Wiki.jsp?page=EISPointing date_obs = gt_fkw(*(self->get(/header)),'DATE_OBS') offset = eis_aia_offsets(date_obs) return, (keyword_set(x)) ? offset[0] : offset[1] END FUNCTION eisobs::_xcen self->private data = *(self->get(/data)) xcen = data->getxcen(/raster) x_offset = self->_eis_aia_offsets(/x) xcen += x_offset return,xcen END FUNCTION eisobs::_image_origin_y self->private ycen = self->_ycen() fovy = gt_fkw(*(self->get(/header)),'FOVY') return, ycen-fovy/2. END FUNCTION eisobs::_ycen self->private data = *(self->get(/data)) n = self->get(/nicons) ycen = fltarr(n) FOR i=0,n-1 DO ycen[i] = data->getycen(i,/raster) y_offset = self->_eis_aia_offsets(/y) ycen += y_offset return,ycen END FUNCTION eisobs::_image_origin_time_y o = self o->private nicons = o->get(/nicons) time_y = findgen(nicons,2) ;; Get the time and y coordinates from the fits header time_y[*,0] = (*(self->get(/image_scale)))[0]/2 time_y[*,1] = o->_image_origin_y() return, time_y END FUNCTION eisobs::_image_origin_wavelength_y o = self o->private nicons = o->get(/nicons) origin_y = o->_image_origin_y() header = *(o->get(/header)) wavelength_y = fltarr(nicons,2) FOR i=0,nicons-1 DO $ wavelength_y[i,*] = [gt_fkw(header,'TWMIN'+trim(i+1)), origin_y[i]] return, wavelength_y END ;; ;; ----------------------------------------------------- ;; get(/prep) THEN BEGIN file = file_search(o->find_prepfile('level1')+'*') IF file EQ '' THEN $ eis_prep, o.file, /default, outdir=o->get(/destdir), $ refill=o->get(/refill), datfile=file, errfile=errorfile, $ /noabs, /quiet, /save $ ELSE BEGIN ; file = o->find_prepfile('level1') errorfile = o->find_prepfile('error') ENDELSE o.errorfile = errorfile o.level1file = file ENDIF ELSE file = o.file ; data will be a pointer to an object reference ; header will be a pointer to a string array o->set, data=ptr_new(obj_new('eis_data',file)) o->set, header=ptr_new((*(o.data)->gethdr())->getbte_hdr()) END FUNCTION eisobs::_gt_eis_linelist_file, backup=backup default, backup, 0 ;; The main linelist is a manually edited main linelist based on ;; Young et al 2007, PASJ 59, 857-864 ;; The backup linelist is based the huge MSSL linelist with data from ;; QS, AR and flare. IF ~backup THEN linelist = 'eis_linelist_PYoung.txt' ELSE $ linelist = 'eis_linelist_MSSL.txt' file = '$SSW/hinode/eis/idl/atest/osdc/'+linelist IF file_test(file) THEN return, file ELSE BEGIN message,'Line list file has been moved. Searching...',/info return, file_search('$SSW/hinode/eis',linelist) ENDELSE END FUNCTION eisobs::_gt_commonblock_linelist, backup=backup self->private default, backup, 0 COMMON eis_linelists_sdc, list_main, list_backup IF (~backup && n_elements(list_main) EQ 0) OR $ (backup && n_elements(list_backup) EQ 0) THEN BEGIN IF ~backup THEN BEGIN template = { version:1.0,$ datastart:4L,$ delimiter:44b,$ missingvalue:!values.f_nan,$ commentsymbol:'',$ fieldcount:2L,$ fieldtypes:[4L,7],$ fieldnames:['lam','ion'],$ fieldlocations:[0L,7],$ fieldgroups:lindgen(2)} list_main = $ read_ascii(self->_gt_eis_linelist_file(),template=template) ENDIF ELSE BEGIN template = { version:1.0,$ datastart:3L,$ delimiter:44b,$ missingvalue:!values.f_nan,$ commentsymbol:'',$ fieldcount:3L,$ fieldtypes:[4L,7,4],$ fieldnames:['lam','ion','arcount'],$ fieldlocations:[0L,7,12],$ fieldgroups:lindgen(3)} list_backup = $ read_ascii(self->_gt_eis_linelist_file(/backup),template=template) ENDELSE ENDIF return, (~backup) ? list_main : list_backup END FUNCTION eisobs::find_prepfile, preptype, i=i, splinelambda0=splinelambda0, $ kamio=kamio file = prune_filename(self->get(/file)) kamiotxt = (keyword_set(kamio)) ? 'kamio_' : '' CASE preptype OF 'level1': file = 'eis_l1_' + strmid(file,7,15)+'.fits' 'error': file = 'eis_er_' + strmid(file,7,15)+'.fits' 'fit' : file = 'eis_fit' + trim(i) + '_'+kamiotxt+strmid(file,7,15)+'.save' 'adef' : file = 'eis_l1_' + strmid(file,7,15) +'_'+kamiotxt+trim(i)+'.adef' 'ana' : file = 'eis_l1_' + strmid(file,7,15) +'_'+kamiotxt+trim(i)+'.ana' 'lamcorr':file = 'eis_l1_' + strmid(file,7,15)+ '_'+ kamiotxt+'wavecorr.save' 'spline': file = self->_find_bestsplinefile(splinelambda0=splinelambda0) ELSE: file = -1 ENDCASE return, self->get(/destdir) + file END FUNCTION eisobs::_find_bestsplinefile, splinelambda0=splinelambda0 ;; Returns the name of the file containing the best fitted spline. For now ;; we hardcode this line to be the Fe XII 195.12 A line. self->private splinelambda0 = 195.12 file = self->get(/file) i_FeXII_195 = eis_get_iwin(file, splinelambda0) IF i_FeXII_195 NE -1 THEN $ file = 'eis_fit' + trim(i_FeXII_195) + '_' + $ strmid(prune_filename(file),7,15)+'.save' $ ELSE message,'Hm. No Fe XII 195 in this file...' return,file END ;; ----------------------------------------------------- ;; is_set(/data) THEN o->read_data returndata = (n_params() EQ 2) ? 1 : 0 o->_image_number_exist,i slit = o->get(/slit) IF ~slit THEN im = o->_extract_image_slot_raster(i) ELSE BEGIN extraset = (size(extra,/tname) EQ 'UNDEFINED' ) ? 0 : 1 IF o->get(/firstdim) EQ 'Wavelength (Å)' THEN BEGIN IF n_elements(ii) EQ 0 THEN im = o->_extract_image_slit_lambda(i) ELSE BEGIN im = -1 message,'Gaussians will not be fitted to this data set',/info ENDELSE IF extraset THEN undefine,extra ENDIF ELSE $ IF extraset THEN $ im = o->_extract_image_slit_gaussian(i, ii = ii, gausslam = $ gausslam, _extra=extra, $ errorinfo=errorinfo) $ ELSE $ im = o->_extract_image_slit_exp(i) ENDELSE IF arg_present(medianvalue) THEN IF (size(im))[0] NE 0 THEN BEGIN IF size(ii,/tname) EQ 'UNDEFINED' THEN medianvalue = median(im) ELSE BEGIN expgood = o->get(/expgood) medianvalue = (n_elements(expgood) GT 1) ? median(im[expgood,*]) : 0 ENDELSE ENDIF IF returndata THEN data = im ELSE $ im = o->_image_final_touch(im, errorinfo=errorinfo, maxvel=maxvel, $ errvelmax=errvelmax, $ gausslam = gausslam, $ _extra=extra) return, im END FUNCTION eisobs::_extract_image_slot_raster, i o = self o->private eisdata = (*(o->get(/data))) eisdata->getwin,i,data,pos nraster = eisdata->getnraster() lambda = eisdata->getlambda() lam = [lambda.scale_b,lambda.scale_a] dx_size = eisdata->getdx_size() rot = (o->get(/widewindows)) ? 0 : 5 ; Do not rotate if widewindow im = rotate(eis_slotwin_oshadow(data,pos,nraster,dx_size,lam),rot) return, sigrange(im) END FUNCTION eisobs::_extract_image_slit_exp, i o = self o->private eisdata = (*(o->get(/data))) eisdata->getwin,i,data,pos sz = size(data) ;; Find pixelnumber in wavelength dimension where mean emission line peaks meanline = fltarr(sz[1]) FOR y=0,sz[2]-1 DO FOR exp=0,sz[3]-1 DO meanline += data[*,y,exp] dummy = max(meanline, mx) integratedintensity = total(data[(mx-5) > 0 : (mx+5) < sz[1]-1, *, *],1) ;; Get the number of exposure times per raster position nexp_prp = (*(o->get(/data)))->getnexp_prp() ;; If nexp_prp gt 1 im is at first a 3D cube im = fltarr(sz[3]/nexp_prp, sz[2], nexp_prp) ;; EIS scans from west to east, i.e. for a raster image the first exposure ;; is the west most. When plotting the image the y*x data array must be ;; rotated with rotate=1. For sit-and-stare-observations the first exposure ;; should be the first column in the 2D time vs y image, and the y*x data ;; array must be rotated with rotate=4. Fill image 2D/3D cube (use ;; [0,0,nexp] instead of [*,*,nexp] to speed up ;fmir_ss = gt_fkw(*(o->get(/header)),'FMIR_SS') ;rot = (fmir_ss NE 0) ? 1 : 4 rot = (o->get(/raster)) ? 1 : 4 FOR nexpct=0,nexp_prp-1 DO im[0, 0, nexpct] = $ rotate(reform(integratedintensity[*, nexpct:*:nexp_prp]),rot) ;; Get rid of 3rd dimension (due to nexp_prp) by using rebin. im will now be ;; the mean image of all nexp_prp exsposures im = rebin(im, (size(im))[1], (size(im))[2], 1) return, sigrange(im) END FUNCTION eisobs::_extract_image_slit_gaussian, i, ii=ii, slow=slow, $ gausslam=gausslam, _ref_extra=extra self->private IF keyword_set(slow) THEN $ im = self->_extract_image_slit_eis_auto_fit(i, _extra = extra) ELSE $ im = self->_extract_image_slit_cfit(i, ii=ii, gausslam=gausslam,$ _extra = extra) return,im END FUNCTION eisobs::_extract_image_slit_cfit, i, ii=ii, errorinfo=errorinfo, $ intensity=intensity, velocity=velocity, width=width, lambda=lambda, $ newadef=newadef, inspect=inspect, gausslam=gausslam, $ fittwice=fittwice, snratio=snratio, kamio=kamio, medianvelcorr=medianvelcorr,$ quiet=quiet ;;+ ;; Returns 2D image of a specified gauss parameter of gauss component ;; number ii in data window number i. Check if there exists a wavelength ;; correction file, an .adef definition file for data window i and a .ana ;; file for window i. Restore if exsisting, generate if not. ;; ;; INPUT: ;; i: extract image from data window number i ;; KEYWORDS: ;; ii: extract image of a gauss parameter from gauss profile number ii ;; errorinfo: output from gt_snratio ;; intensity/velocity/width: extract this gauss parameter ;; lambda: IF set to 1: return wavelength array, IF set to 2 return ;; wavelength array where only the slit tilt is corrected for. ;; newadef: if set to 1, automatically generate new analysis ;; definition, if set to 2 interactively generate one ;; inspect: call xcfit_block to inspect the results of the cfit ;; gausslam: median wavelength of gauss component ii ;; fittwice: call cfit twice, the second time with the result from ;; the first fit as input ;; snratio: the signal to noise ratio to get from gt_snratio ;; medianvelcorr; after the fit is done, subtract the west/east ;; velocity trend from the image. CAUTION! USE WITH CARE! THIS MIGHT ;; HIDE BAD FITTING! ;; quiet: reduce the number of messages to the terminal ;; OUTPUT: 2D array ;; ;; We don't know how many lines there are in data window i, therefore we ;; will fit an arbitrary number of gaussians (or a maximum number set in ;; gt_peaks (called by gt_adef/auto_adef) using S.V.H. Haugan's cfit system. o = self o->private default, kamio, 1 default, newadef, 0 ; Get data and wavelength scale for window i eis_data = *(o->get(/data)) windata = eis_data->getwindata(i) data = windata.int sz = size(data) ;; Data may be 4D, let's get rid of the fourth dimension (that is caused by ;; nexp_prp gt 1). The best thing would of course be to be able to plot an ;; image of each of the nexp_prp exposures, and for any serious scientific ;; use this is the only way to do it. But this is a simple quick look ;; software. data = rebin(data, sz[1], sz[2], sz[3]) sz = size(data) lam = windata.wvl ;; Cannot use 1D wavelength scale for EIS data due to the orbital variation ;; of the line position and tilt of slit; must create 3D wavelength cube ;; where each pixel has its own 1D wavelength scale. IF keyword_set(kamio) THEN BEGIN lamcorr = eis_data->gethkwavecorr(median(lam), /wvl_cube) IF typename(lamcorr) EQ 'INT' THEN BEGIN ;; gethkwavecorr returned -1, i.e. eis_data.sethkpixcorr failed, ;; i.e. the wavelength correction failed. We still want to be able to ;; display a velocity image, so we skip wavelength correction lam3D = fltarr(sz[1],sz[2],sz[3]) FOR j = 0,sz[1]-1 DO lam3D[j,*,*] = lam[j] box_message,'DISPLAYING VELOCITY IMAGE WITHOUT WAVELENGTH CORRECTION' ENDIF ELSE BEGIN ;; When nexp_prp gt 1 the lamcorr.cube is 4D and lamcorr.time is 2D. We ;; take the easy way out and get rid of the extra dimension, not taking ;; into account the multiple exposures per raster position. lam3D = rebin(lamcorr.cube, sz[1],sz[2],sz[3]) lamcorrtime = rebin(lamcorr.time, sz[2]) ENDELSE ENDIF ELSE BEGIN lamcorr = o->_gt_warren_lamcorr(lamcorrtime = lamcorrtime, $ lam_i = median(lam)) lamcorr = rebin(lamcorr, sz[2], sz[3]) lam3D = fltarr(sz[1],sz[2],sz[3]) FOR j = 0,sz[1]-1 DO lam3D[j,*,*] = lam[j] - lamcorr ENDELSE ana = o->_gt_cfit_ana(data, lam3D, i, newadef=newadef, fittwice=fittwice, $ kamio = kamio) IF keyword_set(inspect) THEN BEGIN ;; Visually and interactively inspect the result loadct,3 xcfit_block,ana=ana ENDIF missing = ana.missing IF n_elements(snratio) NE 0 THEN errorinfo = gt_snratio(ana,missing=missing) handle_value, ana.result_h, result handle_value, ana.fit_h, fit ;; We want gaussian profile number ii. Check if gaussian profile number ;; exist, if so translate ii to correct ana.result_h index, if not return to ;; caller ii_set = (n_elements(ii) EQ 1) ? 1 : 0 IF ii_set THEN o.ii = ii IF ii_set && ii GT 0 THEN BEGIN IF NOT o->_check_cfit_sub_image_number_exists(ii, result, ngaussians, $ quiet=quiet) THEN $ BEGIN delete_analysis, ana return,-1 ENDIF gaussind = ii*3 ENDIF ELSE gaussind = 0 delete_analysis,ana ;; extract the rest wavelength of the fitted profile (may be returned to caller) medianlam = fltarr(((size(result))[1]-2)/3) FOR j=0,n_elements(medianlam)-1 DO medianlam[j] = median(result[j*3+1,*, *]) gausslam = o->_extract_cfit_gausslam(fit, gaussind, medianlam) ;; extract the I/v/w image im = o->_extract_cfit_im(eis_data, gaussind, result, missing, $ intensity=intensity, velocity=velocity, width=width, $ lambda = lambda, lamcorrtime = lamcorrtime, $ medianvelcorr=medianvelcorr) return,im END FUNCTION eisobs::_gt_warren_lamcorr, lamcorrtime=lamcorrtime, lam_i=lam_i o = self o->private lamcorrfile = o->find_prepfile('lamcorr') IF file_test(lamcorrfile) THEN restore, lamcorrfile ELSE BEGIN ;; Correct for the slit tilt and orbit variation If the duration of the ;; raster is less than 300 s, eis_wave_corr_old will crash when trying to ;; fit a spline, i.e. set keyword nospline header = *(o->get(/header)) time = anytim2tai(gt_fkw(header,'DATE_END')) - $ anytim2tai(gt_fkw(header,'DATE_OBS')) nospline = (time LT 600) ? 1 : 0 eis_wave_corr_old, o->get(/level1file), lamcorrtiltshort, $ lamcorrtime, nospline=nospline ;,/do_fit ;; eis_wave_corr_old returns the slit tilt for lines at the short ;; wavelength detector. We need to calculate the slit for the other ;; detector as well, and save both tilts eis_data = *o.data yws = eis_data->getinfo('yws') ; y window start ny = median((eis_data->getyw())) ; number of elements along y date = eis_data->getinfo('date_obs') slit = eis_data->getinfo('slit_ind') lamcorrtiltlong = eis_slit_tilt(yws, ny, date=date, $ slit=eis_data->getinfo('slit_ind'), /long) IF (size(where(finite(lamcorrtiltshort))))[0] EQ 0 THEN $ lamcorrtiltshort[*] = 0 IF (size(where(finite(lamcorrtiltlong))))[0] EQ 0 THEN $ lamcorrtiltlong[*] = 0 IF (size(where(finite(lamcorrtime))))[0] EQ 0 THEN lamcorrtime[*] = 0 nx = n_elements(lamcorrtime) ny = n_elements(lamcorrtiltshort) lamcorrshort = fltarr(nx, ny) lamcorrlong = fltarr(nx, ny) FOR i=0,nx-1 DO BEGIN lamcorrshort[i,*] = lamcorrtiltshort lamcorrlong[i,*] = lamcorrtiltlong ENDFOR FOR j=0,ny-1 DO BEGIN lamcorrshort[*,j] += lamcorrtime lamcorrlong[*,j] += lamcorrtime ENDFOR tai = (*(o->get(/data)))->ti2tai() save,file=lamcorrfile, lamcorrshort, lamcorrlong, lamcorrtiltshort, $ lamcorrtiltlong, lamcorrtime, tai END lamcorr = (lam_i GT 230) ? lamcorrlong : lamcorrshort return, lamcorr END FUNCTION eisobs::_gt_cfit_ana, data, lam, i, $ newadef=newadef, fittwice=fittwice, kamio=kamio ;;+ ;; Create or restore the ana file for the data set. ;;- o = self o->private sz = size(data) missing = o->get(/missing) anafilename = o->find_prepfile('ana',i = i, kamio = kamio) IF keyword_set(newadef) OR ~file_test(anafilename) THEN BEGIN adeffile = o->find_prepfile('adef',i=i, kamio = kamio) IF keyword_set(newadef) OR ~file_test(adeffile) THEN BEGIN ;; Interactively create an analysis definition (newadef eq 2) or let ;; gt_adef find/generate one (newadef eq 1 or adeffile does not exist) meanprofile = total(total(data,2),2)/(sz[2]*sz[3]) IF newadef EQ 2 THEN BEGIN l = median(median(lam,dim=2),dim=2) ;; median 1D lam array xcfit,l,meanprofile,adef ENDIF ELSE $ adef = o->_auto_adef(lam,meanprofile) save,file=adeffile,adef ENDIF ELSE restore,file=adeffile wts = 1./(data > 1) ; error is proportional to sqrt of counts ;; Perform the cfit! ana = mk_analysis(lam,data,wts,adef,missing,result,residual) cfit_block, analysis=ana, /quiet IF keyword_set(fittwice) THEN cfit_block, analysis=ana, /quiet ana.filename = anafilename save_analysis,ana ENDIF ELSE ana = restore_analysis(anafilename) return, ana END FUNCTION eisobs::_auto_adef, lam, int ;; Automatically generate cfit analysis definitions based on input intensity and ;; wavelength arrays self->private peakinds = gt_peaks(int, fwhm=fwhm, minmedian=4.5,/sort,/plot) npeaks = n_elements(peakinds) gaussians = replicate(mk_comp_gauss([0,0,0]), npeaks) int0 = int[peakinds] lampeak = lam[peakinds,*,*] lamfwhm = lam[peakinds-fwhm,*,*] lam0 = fltarr(npeaks) wid0 = fltarr(npeaks) FOR i=0,npeaks-1 DO lam0[i] = median(lampeak[i,*,*]) FOR i=0,npeaks-1 DO wid0[i] = lam0[i] - median(lamfwhm[i,*,*]) v = 75. ; Max shift in km/s dlam = v*lam0/3.e5 ; Max shift in Aangstrom intmin = fltarr(npeaks) ; minimum intensity is 0 lammin = (lam0 - dlam) > min(lam); v0 - v widmin = (wid0 - 0.04) > 0.02 ; random guess... intmax = int0*100;30000 ; More random guessing lammax = (lam0 + dlam) < max(lam) ; v0 + v widmax = wid0 + 0.04 ; A final shot in the dark FOR i=0,n_elements(peakinds)-1 DO BEGIN gauss = mk_comp_gauss([int0[i],lam0[i],wid0[i]], $ max_arr=[intmax[i],lammax[i],widmax[i]], $ min_arr=[intmin[i],lammin[i],widmin[i]], $ trans_a=[1,1,0.424661], trans_b=[0,0,0], $ const=[0b,0b,0b]) lam0txt = trim(lam0[i],'(F6.2)') gauss.name = 'AutoGauss'+lam0txt gaussians[i] = gauss ENDFOR bg = mk_comp_poly([0.5*median(int)], max_arr=[30000],min_arr=[0],trans_a=[1],$ trans_b=[0],const=[0b]) bg.name = 'Background' IF npeaks EQ 1 THEN adef = {igauss2:gaussians[0], bg:bg} IF npeaks EQ 2 THEN adef = {igauss2:gaussians[0], igauss3:gaussians[1], bg:bg} IF npeaks EQ 3 THEN adef = {igauss2:gaussians[0], igauss3:gaussians[1], $ igauss4:gaussians[2], bg:bg} IF npeaks EQ 4 THEN adef = {igauss2:gaussians[0], igauss3:gaussians[1], $ igauss4:gaussians[2], igauss5:gaussians[3], bg:bg} IF npeaks EQ 5 THEN adef = {igauss2:gaussians[0], igauss3:gaussians[1], $ igauss4:gaussians[2], igauss5:gaussians[3], $ igauss6:gaussians[4], bg:bg} IF npeaks EQ 6 THEN adef = {igauss2:gaussians[0], igauss3:gaussians[1], $ igauss4:gaussians[2], igauss5:gaussians[3], $ igauss6:gaussians[4], igauss7:gaussians[5], bg:bg} IF npeaks EQ 7 THEN adef = {igauss2:gaussians[0], igauss3:gaussians[1], $ igauss4:gaussians[2], igauss5:gaussians[3], $ igauss6:gaussians[4], igauss7:gaussians[5], $ igauss8:gaussians[6], bg:bg} IF npeaks EQ 8 THEN adef = {igauss2:gaussians[0], igauss3:gaussians[1], $ igauss4:gaussians[2], igauss5:gaussians[3], $ igauss6:gaussians[4], igauss7:gaussians[5], $ igauss8:gaussians[6], igauss9:gaussians[7], $ bg:bg} IF npeaks GE 9 THEN adef = {igauss2:gaussians[0], igauss3:gaussians[1], $ igauss4:gaussians[2], igauss5:gaussians[3], $ igauss6:gaussians[4], igauss7:gaussians[5], $ igauss8:gaussians[6], igauss9:gaussians[7], $ igauss10:gaussians[8], bg:bg} IF npeaks GE 9 THEN box_message,'Found '+trim(npeaks)+', fitting the 9 highest...' return, adef END FUNCTION eisobs::_extract_cfit_gausslam, fit, gaussind, medianlam ;;+ ;; extract the rest wavelength of the fitted profile ;;- self->private CASE gaussind OF ;; If the name of the gaussian component doesn't start with AutoGauss, ;; the component is not automatically genererated by the eisobs object. 0: autogauss = (strmid(fit.igauss2.name,0,9) EQ 'AutoGauss') ? 1 : 0 3: autogauss = (strmid(fit.igauss3.name,0,9) EQ 'AutoGauss') ? 1 : 0 6: autogauss = (strmid(fit.igauss4.name,0,9) EQ 'AutoGauss') ? 1 : 0 9: autogauss = (strmid(fit.igauss5.name,0,9) EQ 'AutoGauss') ? 1 : 0 12: autogauss = (strmid(fit.igauss6.name,0,9) EQ 'AutoGauss') ? 1 : 0 15: autogauss = (strmid(fit.igauss7.name,0,9) EQ 'AutoGauss') ? 1 : 0 18: autogauss = (strmid(fit.igauss8.name,0,9) EQ 'AutoGauss') ? 1 : 0 21: autogauss = (strmid(fit.igauss9.name,0,9) EQ 'AutoGauss') ? 1 : 0 24: autogauss = (strmid(fit.igauss10.name,0,9) EQ 'AutoGauss') ? 1 : 0 27: autogauss = (strmid(fit.igauss11.name,0,9) EQ 'AutoGauss') ? 1 : 0 ELSE: autogauss = '10 gaussians ought to do it!' ENDCASE IF autogauss THEN BEGIN ;; The component is automatically generated, and we know that the ;; wavelength of the line is a part of the fit.igaussX.name: CASE gaussind OF 0: gausslam = strmid(fit.igauss2.name,9,6) 3: gausslam = strmid(fit.igauss3.name,9,6) 6: gausslam = strmid(fit.igauss4.name,9,6) 9: gausslam = strmid(fit.igauss5.name,9,6) 12: gausslam = strmid(fit.igauss6.name,9,6) 15: gausslam = strmid(fit.igauss7.name,9,6) 18: gausslam = strmid(fit.igauss8.name,9,6) 21: gausslam = strmid(fit.igauss9.name,9,6) 24: gausslam = strmid(fit.igauss10.name,9,6) 27: gausslam = strmid(fit.igauss11.name,9,6) ELSE: gausslam = '10 gaussians ought to do it!' ENDCASE ENDIF ELSE BEGIN CASE gaussind OF ;; The component is NOT automatically generated, and we need to find ;; the wavelength from the param structure: 0: gausslam = trim(medianlam[0]) 3: gausslam = trim(medianlam[1]) 6: gausslam = trim(medianlam[2]) 9: gausslam = trim(medianlam[3]) 12: gausslam = trim(medianlam[4]) 15: gausslam = trim(medianlam[5]) 18: gausslam = trim(medianlam[6]) 21: gausslam = trim(medianlam[7]) 24: gausslam = trim(medianlam[8]) 27: gausslam = trim(medianlam[9]) ELSE: gausslam = '10 gaussians ought to do it!' ENDCASE ENDELSE return, gausslam END FUNCTION eisobs::_extract_cfit_im, eis_data, gaussind, result, missing, $ intensity=intensity, velocity=velocity, $ width=width, lambda=lambda, $ lamcorrtime = lamcorrtime, $ medianvelcorr=medianvelcorr ;;+ ;; Return the intensity/velocity/linewidth or wavelength array. If lamba is ;; set to 1 return the wavelength array where slit tilt and the variation of ;; the line centre due to the orbital motion of the satellite is corected ;; for. If lambda is set to 2 only the slit tilt is subtracted from the ;; data. If keyword velocity is set then any remaining velocity trend along ;; the x axis is subtracted if also medianvelcorr is set. ;;- self->private ;; Find the number of exposures in the ;; observations and the indices of missing exposures (if any). Set ;; missing exposures to badvalue. imagecolumnnumber = (self->get(/firstdim) NE 'Time (min)') ? 1 : 0 expmissing = self->get_expmissing(imagecolumnnumber=imagecolumnnumber) dummy = where(expmissing GT -1, nexpmissing) IF keyword_set(intensity) THEN BEGIN im = reform(result[gaussind+0,*,*]) IF nexpmissing GT 0 THEN im[expmissing,*] = missing ENDIF IF keyword_set(width) THEN BEGIN im = reform(result[gaussind+2,*,*]) IF nexpmissing GT 0 THEN im[expmissing,*] = missing ENDIF IF keyword_set(velocity) THEN BEGIN lam = reform(result[gaussind+1,*,*]) IF nexpmissing GT 0 THEN lam[expmissing,*] = missing notfinite = where(~finite(lam),count) IF count GT 0 THEN lam[where(~finite(lam))] = missing ;; goodpix are all pixels that are finite, not marked as missing in the ;; ana file or not in a missing column according to eis_data->check_ti(): goodpix = where(lam NE missing, ngoodpix, $ complement=badpix, ncomplement=nbadpix) IF ngoodpix GT 0 THEN lam_med = median(lam[goodpix]) ELSE lam_med = 0 ;; Wavelength or velocity? IF NOT keyword_set(lambda) THEN im = 3.0e5*(lam-lam_med)/lam_med ELSE BEGIN IF lambda EQ 1 THEN im = lam IF lambda EQ 2 THEN BEGIN sz = size(result) im = self->_extract_cfit_uncorrected_lambda(lam, lamcorrtime, sz) ENDIF ENDELSE IF keyword_set(medianvelcorr) THEN BEGIN IF keyword_set(lambda) THEN $ message,'Keyword lambda ignored when keyword medianvelcorr is set',$ /info ELSE BEGIN IF ngoodpix EQ 0 THEN BEGIN im[*,*] = 0 ;; Set to 0 if all is Inf or NaN or missing badpix = where(im,nbadpix) ENDIF IF nbadpix GT 0 THEN im[badpix] = 0 sz = size(im) v = median(im,dimension=2) vim = fltarr(sz[1],sz[2]) FOR j=0,sz[2]-1 DO vim[*,j] = v im -= vim ENDELSE ENDIF IF nbadpix GT 0 THEN im[badpix] = missing ENDIF notfinite = where(~finite(im),count) IF count GT 0 THEN im[where(~finite(im))] = missing return, im END FUNCTION eisobs::_extract_cfit_uncorrected_lambda, lam, lamcorrtime, sz ;; Reintroduce the time dependent wavlength correction that has previously ;; been subtracted from the data. The 2D lam array is now corrected for slit ;; tilt only. o = self o->private lamcorrtime2D = dblarr(sz[2],sz[3]) FOR i =0,sz[3]-1 DO lamcorrtime2D[*,i] = lamcorrtime lam += lamcorrtime2D return, lam END FUNCTION eisobs::_check_cfit_sub_image_number_exists, ii, result, ngaussians, $ quiet=quiet self->private ;; a cfit result block consists of I, v and w for each component pluss ;; background and residuals. The number of gaussians is then the number of 2D ;; images in result - 2 (background and residuals) divided by 3 (I, v and w). ngaussians = ((size(result))[1]-2)/3 exist = (ii GE ngaussians) ? 0 : 1 IF NOT keyword_set(quiet) THEN IF NOT exist THEN $ message, 'Gaussian profile number must be in the range [0,'+$ trim(ngaussians-1)+']',/info return,exist END FUNCTION eisobs::_extract_image_slit_eis_auto_fit, i, intensity=intensity, $ velocity=velocity, width=width, errorinfo=errorinfo, _ref_extra=extra o = self o->private ;; Can't trust the header: TWAVEi is often eq TWMIN; must guess where the ;; line actually is... First guess is in the middle of the data window. header = *(o->get(/header)) twmin = gt_fkw(header,'TWMIN'+trim(i+1)) lambda0 = twmin + (gt_fkw(header,'TWMAX'+trim(i+1)) - twmin)/2. ;; Check if gaussians and orbit spline already have been fitted to the ;; data. If so restore the file, if not do the fitting. gaussfitfile = o->find_prepfile('fit',i=i, _extra=extra) IF ~file_test(gaussfitfile) THEN BEGIN ;; Extract data, perform eis_auto_fit and slit tilt correction. Get the ;; spline structure of the strongest line (either by fitting or restoring ;; file). Use the best spline as input, or fit spline if i is the ;; strongest line. Save fit and spline structures. fit = o->_fit_and_correct(lambda0) spline = o->_best_spline(lambda0) ;spline eq -1 if i is the strongest line IF size(spline,/tname) NE 'STRUCT' THEN undefine,spline fit = eis_orbit_spline(fit,/setzero, spline=spline, /plotvel) save,file=gaussfitfile, fit, spline ENDIF ELSE restore, gaussfitfile, /relaxed_structure_assignment IF keyword_set(intensity) THEN im = fit.int IF keyword_set(velocity) THEN im = fit.vel IF keyword_set(width) THEN im = fit.wid errorinfo = fit undefine,fit return,im END FUNCTION eisobs::_fit_and_correct, lambda0 ;; Extract data for window i, perform eis_auto_fit and slit tilt correction o = self o->private data = eis_getwindata(o->get(/level1file), lambda0) ; data is a structure wvlpix = o->_gt_wvlpix(lambda0) IF n_elements(wvlpix) NE 2 THEN undefine,wvlpix eis_auto_fit, data, fit, refwvl=lambda0, wvlpix=wvlpix fit = eis_tilt_correction(data, fit) return, fit END FUNCTION eisobs::_gt_wvlpix, lambda0 self->private l0 = round(lambda0*10)/10. IF l0 GT 256.00 AND l0 LT 256.5 THEN w = [256.10,256.55] ELSE $ ; He II 256.32 IF l0 GT 192.5 AND l0 LT 193 THEN w = [192.70,192.90] ELSE $ ; Ca XVII 192.82 w = -1 return,w END FUNCTION eisobs::_best_spline, lambda0 ;; Return the spline of the strongest line in the data set. o = self o->private bestsplinefile = o->find_prepfile('spline',i=i, splinelambda0=splinelambda0) file = o->get(/file) ilam0 = eis_get_iwin(file, lambda0) ispllam0 = eis_get_iwin(file, splinelambda0) ;; Check if current line is the strongest line. If so, set spline to -1, if ;; not restore spline file if it already exist or fit if if doesn't exist. IF ilam0 EQ ispllam0 THEN spline = -1 ELSE BEGIN IF file_test(bestsplinefile) THEN restore,bestsplinefile ELSE BEGIN fit = o->_fit_and_correct(splinelambda0) fit = eis_orbit_spline(fit,/setzero, spline=spline, /plotvel) save,file=o->find_prepfile('fit',i=ispllam0), fit, spline ENDELSE ENDELSE return, spline END FUNCTION eisobs::_extract_image_slit_lambda, i o = self o->private eisdata = (*(o->get(/data))) eisdata->getwin,i,data,pos sz = size(data) nexp_prp = eisdata->getnexp_prp() IF nexp_prp GT 1 THEN $ im = reform((rebin(data, sz[1], sz[2], sz[3]/nexp_prp))[*,*,0]) ELSE $ im = data im = reform(im[*,*,0]) ; We will display the first of n images return, sigrange(im, fraction=0.99) ; default fraction is 0.99 END FUNCTION eisobs::_image_final_touch, im, errorinfo=errorinfo, maxvel=maxvel, $ errvelmax=errvelmax, snratio=snratio, $ intensity=intensity, $ velocity=velocity, width=width, $ gausslam=gausslam, kamio=kamio self->private IF (size(im))[0] EQ 0 THEN return,-1 ;; return if im isn't a 2D image missing = self->get(/missing) goodpix = where(im NE missing, ngoodpix, $ complement=badpix, ncomplement=nbadpix) IF n_elements(snratio) EQ 0 THEN BEGIN ;; ::plot_image keyword snratio is not set; display all pixels of the ;; image IF ngoodpix EQ 0 THEN BEGIN im[*,*] = 0 ;; Set to 0 if all is Inf or NaN or missing badpix = where(im,nbadpix) ENDIF IF nbadpix GT 0 THEN BEGIN IF ngoodpix GT 0 THEN med = median(im[goodpix]) ELSE med = 0 im[badpix] = med ENDIF errtxt = '0' ENDIF IF ~(keyword_set(intensity) OR keyword_set(velocity) OR keyword_set(width)) $ THEN BEGIN im = sigrange(im) self->set, dynamictitle='' self->set, rightytitle='' im = 255-bytscl(im) IF nbadpix GT 0 THEN im[badpix] = 255 ENDIF ELSE BEGIN IF n_elements(snratio) NE 0 THEN BEGIN IF NOT keyword_set(errorinfo) THEN return,-1 ;; The image is based on gaussian fitting of the line profiles, either ;; using P. Young's eis_auto_fit or S.V.H. Haugan's cfit. cfit = (size(errorinfo,/tname) EQ 'STRUCT') ? 0 : 1 IF ~cfit THEN BEGIN default, errvelmax, 8 errtxt = errvelmax velerr = errorinfo.velerr badpix = where(velerr GT errvelmax OR velerr EQ errorinfo.missing, $ nbadpix, complement=goodpix, ncomplement=ngoodpix) ENDIF ELSE BEGIN default, snratio, 3 errtxt = snratio missing = self->get(/missing) badpix = where(errorinfo LT snratio OR errorinfo EQ missing, nbadpix,$ complement=goodpix, ncomplement=ngoodpix) ENDELSE ENDIF IF keyword_set(velocity) THEN BEGIN default, maxvel, 15 ; Velocities GT 15 or LT -15 is red/blue IF ngoodpix GT 0 THEN BEGIN im -= median(im[goodpix]) ;; Adjust max/min vel to be used when displaying image if necessary. mx = maxvel < (abs(min(im[goodpix])) < max(im[goodpix])) ;; Set velocities exceeding max/min limit equal to the limit IF mx NE 0 THEN im = im > (-mx) < mx ELSE $ im = im > min(im[goodpix]) < max(im[goodpix]) ;mx = 40;!!!!!!!!!!!!!!!!!!!!!!!!!! ENDIF ELSE mx = 0 ;; Set bad pixels to 0 in order not to mess up the bytescaling. Do not ;; use color indices 0 or 255 (they are black and white, not red/blue) IF nbadpix GT 0 THEN im[badpix] = 0 im = bytscl(im) im = im > 1 < 254 self->_dynamictitle_set, gausslam = gausslam self->_rightytitle_set, 'V', 'km s!u-1!n', $ [trim(-round(mx)),trim(round(mx))], errtxt, $ gausslam = gausslam ENDIF IF keyword_set(width) THEN BEGIN IF nbadpix GT 0 THEN BEGIN IF ngoodpix GT 0 THEN med = median(im[goodpix]) ELSE med = 0 im[badpix] = med ENDIF im = im-median(im) im = sigrange(im) absmin = abs(min(im)) < max(im) ix = where(im LT -absmin, count) IF count GT 0 THEN im[ix] = -absmin ix = where(im GT absmin, count) IF count GT 0 THEN im[ix] = absmin self->_dynamictitle_set, gausslam = gausslam self->_rightytitle_set, 'W', 'Å', $ [strmid(trim(min(im[goodpix])),0,5), $ strmid(trim(max(im[goodpix])),0,5)],$ errtxt, gausslam = gausslam im = bytscl(im) im = im > 1 < 254 ENDIF IF keyword_set(intensity) THEN BEGIN IF nbadpix GT 0 THEN BEGIN IF ngoodpix GT 0 THEN med = median(im[goodpix]) ELSE med = 0 im[badpix] = med ENDIF im = sigrange(im) self->_dynamictitle_set, gausslam = gausslam self->_rightytitle_set, '!17I!3', 'erg cm!u-2!ns!u-1!nsr!u-1!n', $ [strmid(trim(min(im[goodpix])),0,5), $ strmid(trim(max(im[goodpix])),0,5)],$ errtxt, gausslam = gausslam im = 255-bytscl(im) END ;; Bad pixels are 255 (white) (vel/width images are bytscl ranging from 1 ;; to 254, in all three colortables are 0=black, 255=white) IF nbadpix GT 0 THEN im[badpix] = 255 ENDELSE ;; If all pixels have the same value, we define them all as bad. Bad pixels! IF max(im) EQ min(im) THEN im[*,*] = 255 return, im END ;; PRO eisobs::_erase_temp_variables self->private self.ii = -1 END ;; PRO eisobs::_cenAndfov, xcen, ycen, fovx, fovy self->private xcen = self->_xcen() ycen = self->_ycen() header = *(self->get(/header)) fovx = gt_fkw(header,'FOVX') fovy = gt_fkw(header,'FOVY') END FUNCTION eisobs::_extract_image_position, sz o = self o->private IF strmid((o->get(/xtitle)),0,5) NE strmid((o->get(/ytitle)),0,5) OR $ (~o->get(/slit) AND o->get(/raster)) THEN $ return, o->_image_position(sz) ELSE return,-1 END FUNCTION eisobs::_extract_image_scale,i o = self o->private IF ~o->get(/slit) THEN return, o->_extract_image_scale_slot(i) $ ELSE BEGIN image_scale = *(o->get(/image_scale)) IF (size(image_scale))[0] EQ 2 THEN image_scale = reform(image_scale[i,*]) return, image_scale ENDELSE END FUNCTION eisobs::_extract_image_scale_slot, i self->private print,'fix eisobs::get_image_scale for slot observations!' return,[1,1] END FUNCTION eisobs::_extract_title,i o = self o->private dynamictitle = o->get(/dynamictitle) IF o->is_set(/ii) THEN title = dynamictitle ELSE BEGIN dynamictitle = (dynamictitle EQ '') ? '' : ' ('+dynamictitle+')' title = (*(o->get(/title)))[i] + dynamictitle ENDELSE return, title END PRO eisobs::_dynamictitle_set, gausslam=gausslam self->private IF keyword_set(gausslam) THEN BEGIN ;; Manually edited linelist, based on Young et al 2007, PASJ 59, 857-864 list = self->_gt_commonblock_linelist() lam = list.lam ion = list.ion ind = gt_closest_index(gausslam, lam, mindiff=mindiff,/noprint) IF abs(mindiff) LT 0.05 THEN BEGIN IF strpos(ion[ind],'/') EQ -1 THEN $ ;; i.e. if line is not a blend self.dynamictitle = ion[ind] + ' ' + trim(lam[ind]) + ' Å' ELSE $ self.dynamictitle = trim(gausslam)+' Å ('+ion[ind]+')' ENDIF ELSE BEGIN ;; Huge MSSL linelist with data from QS, AR and flare. list_backup = self->_gt_commonblock_linelist(/backup) list = list_backup ix = where(list.arcount GT 0) lam = float(list.lam[ix]) ion = list.ion[ix] ind = gt_closest_index(gausslam, lam, mindiff=mindiff,/noprint) IF abs(mindiff) LT 0.05 THEN $ self.dynamictitle = ion[ind] + ' ' + trim(lam[ind]) + ' Å' ELSE $ self.dynamictitle = 'Fit at ' +gausslam+' Å' ENDELSE ENDIF END PRO eisobs::_rightytitle_set, gaussiantxt, units, mx, mxerr, $ gausslam=gausslam IF n_params() EQ 4 THEN BEGIN IF keyword_set(gausslam) THEN $ self.rightytitle = gaussiantxt + $ '!9 e!3 [' + mx[0] + ', ' + mx[1] + '] ' + units + $ ', min(!S!U S!R!S-!R!D N!N)='+trim(mxerr) $ ELSE $ self.rightytitle = ', ' + gaussiantxt + $ '!9 e!3 [' + mx[0] + ', ' + mx[1] + $ '] '+units+', max(!4r!3!D vel' + $ '!n)='+trim(mxerr) + 'km s!u-1!n' ENDIF ELSE self.rightytitle = '' END ;; FUNCTION eisobs::_thumb_final_touch, thumb, intensity=intensity, $ width=width self->private IF keyword_set(intensity) OR keyword_set(width) THEN $ thumb = bytscl(thumb,/nan) return, thumb END ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; <_color_table and help methods> ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++ FUNCTION eisobs::_color_table_txt, i, velocity=velocity, width=width self->private header = self->get(/header) wanted_ct_txt = 'INTENSITY' IF keyword_set(velocity) THEN wanted_ct_txt = 'VELOCITY' IF keyword_set(width) THEN wanted_ct_txt = 'WIDTH' return, wanted_ct_txt END ;; ------------------------------------------------------ ; ;; ------------------------------------------------------ ;; --------------------------------------------------------------------- ;; ;; --------------------------------------------------------------------- FUNCTION eisobs::_read_header_fxmethod, lun self->private fxbopen, lun, self.file, 1, hdr,/reopen return, hdr END FUNCTION eisobs::nmissing ;; Return the number of missing exposures ;dummy = where(*(self->get(/data))->check_ti() EQ 0, count) expmissing = self->get_expmissing() return, n_elements(expmissing) END ;; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;; get_XXX functions for all eisobs structure tags. Each funciton checks if ;; header is a valid pointer, if not, call read_header which again ensures the ;; tag is set. ;; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FUNCTION eisobs::get_slit IF ~ptr_valid(self.header) THEN self->read_header return, self.slit END FUNCTION eisobs::get_raster IF ~ptr_valid(self.header) THEN self->read_header return, self.raster END FUNCTION eisobs::get_widewindows IF ~ptr_valid(self.header) THEN self->read_header return,self.widewindows END FUNCTION eisobs::get_expgood IF ~ptr_valid(self.header) THEN self->read_header return, *(self.expgood) END FUNCTION eisobs::get_expmissing, imagecolumnnumber=imagecolumnnumber IF ~ptr_valid(self.header) THEN self->read_header expmissing = *(self.expmissing) IF expmissing[0] NE -1 THEN IF keyword_set(imagecolumnnumber) THEN $ expmissing = reverse(gt_fkw(*(self->get(/header)), 'NEXP')-1-expmissing) return, expmissing END FUNCTION eisobs::get_singleexp IF ~ptr_valid(self.header) THEN self->read_header return, self.singleexp END FUNCTION eisobs::get_image_scale IF ~ptr_valid(self.header) THEN self->read_header return, self.image_scale END ;; ------------------------------------------------------------------------ ;; ;; ------------------------------------------------------------------------ FUNCTION eisobs::emission_only, i ;; All children of the spectobs class should be able to tell if a ;; window contains pure emission lines only. In the case of EIS this is the ;; case for the entire wavelength range. In a future update of eisobs, when ;; it actually inherits spectobs, this method should be removed and instead ;; the spectobs version of this method will be used. return,1 END PRO eisobs::mk_jpg2000, outdir=outdir message,'JPG2000 generation of EIS files not implemented yet.',/info END PRO eisobs__define INT = 0 STR = '' PTR = ptr_new() d = {eisobs, $ inherits hinobs, $ slit:INT,$ raster:INT,$ widewindows:INT ,$ firstdim:STR, $ singleexp:INT, $ level1file:STR,$ errorfile:STR, $ gaussfitfile:STR, $ missing:INT, $ expgood:PTR,$ expmissing:PTR,$ dynamictitle:STR, $ rmprepfiles:INT, $ rmerrorfile:INT, $ ii:INT, $ refill:INT, $ hk_cmn_dir:STR,$ hk_status_dir:STR} END