;+ ; NAME: spectobs__define ; ; ; ; PURPOSE: Definiton of the spectobs ("spectrograph observation") superclass; ; This class contains spectrograph specific code common for all ; spectrographs, and is inherited by spectrograph classes like ; irisobs and eisobs. Code that is common for all missions and ; instruments can be found in the superclass anyobs. Note that ; spectobs doesn't inherit anyobs (or vice verca), but irisobs/eisobs ; inherit both anyobs and spectobs. ; ; (sdc_util) ; | ; anyobs SPECTOBS ; / \ / \ ; hinobs /\ \ ; / | \ / \ \ ; (...) eisobs irisobs \ ; | \ \ ; irissjiobs irisspectobs ; ; ; METHODS: ; ; ; SPECIAL CALLS: none ; ; ; WRITTEN: Terje Fredvik, UiO, May 2013 ; ; MODIFICATION HISTORY: ; 28. May 2013 Terje Fredvik: First draft ; 05. Nov 2013 Terje Fredvik: lots of restructuring... ; 19. Nov 2013 Terje Fredvik: ::_image_final_touch accepts keyword ; photosphere, if set do not invert ; colour table. ; ::_extract_image_slit_exp: ; skip the outer 5% of the window to ; avoid the wing of very strong lines ; leaking into the window, or spurious ; edge effects. ; 22. Nov 2013 Terje Fredvik: Rewrote code in several methods ; handling sit-and-stare observations ; and raster observations with few ; exposures (i.e. observations that ; should be displayed as a spectrum ; (lam vs y) rather than a raster image ; (X vs y)). ; 20. Jan 2014 Terje Fredivk: New methods ::_get_meanprofile and ; ::_get_instrument_specific_missing_value. ; In ::_auto_adef pick middle pixel as ; the only peak if no peaks are ; found. Minor tweaks. ; 31. Jan 2014 Terje Fredvik: New method ::_dynamictitle_addition. ; ::_extract_cfit_gausslam: extract a ; string of size 7 instead of 6. ; ::_auto_adef: pick middle pixel of ; int array when no peaks are found, ; not middle pixel of lam array. Call ; method getcorrectedwavelengthcube, a ; method that FOR the time being only ; exists in irisspectobs__define ; (ouch!). ; 10. Mar 2014 Terje Fredvik ::extract_image now returns -1 if ; keywords VELOCITY or WIDTH is set, or ; if ii is greater than 0. Else return an ; integrated intensity image. ; 12. Mar 2014 Terje Fredvik Minor changes in plot titles ; 19. Mar 2014 Terje Fredvik New method ::_get_lamdiffmax - returns ; the maximum allowed difference ; between a fitted wavelength and the ; rest wavelength of a line in a line ; list. ; 24. Mar 2014 Terje Fredvik Use strtrim instead of trim when ; trimming fitted wavelength in ; ::_auto_adef to prevent "00" after ; period to be removed. ; 09. Apr 2014 Terje Fredvik ::_extract_image_slit_exp: skip outer ; 15 pixels of window when calculating ; mean line profile instead of 5% of ; window width. ; 25. Apr 2014 Terje Fredvik ::extract_image now uses only non-missing ; values when calculating the median ; value of the image. Call ; _get_meanprofile in ; _extract_image_slit_integrated_intensity. ; 27. Jul 2022 Terje Fredvik Call cfit_block with keyword QUIET set ; ;- ; FUNCTION spectobs::init, file;, _ref_extra=extra o = self not_set = -1 o.slit = not_set o.raster = not_set o.nwin = not_set o.ii = not_set o.missing = -1e3 ; To be used in gaussian fitting o.firstdim = 'MISSING' o.dynamictitle = '' return, 1 END PRO spectobs::cleanup print,'Nothing to clean up' END ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;; <_ensure_spect_properties & help methods, called from a spectobs child> ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ PRO spectobs::_ensure_spect_properties o = self o->private o->_ensure_slit o->_ensure_raster o->_ensure_firstdim END PRO spectobs::_ensure_raster o = self o->private IF ~o->is_set(/raster) THEN BEGIN nraster = o->_gt_nraster() o.raster = (nraster GT 1) ? 1 : 0 END END PRO spectobs::_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 = o->_gt_firstdimraster();'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 ;; ------------------------------------------------------------------------ ;; <_ensure_spect_properties & help methods, called from a spectobs child> ;; ------------------------------------------------------------------------ ;; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;; ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ PRO spectobs::_ensure_nicons ;; nicons is an anyobs tag self->private self.nicons = gt_fkw(*self.header,'NWIN') END PRO spectobs::_ensure_xtitle ;; xtitle is an anyobs tag o = self o->private IF ~o->is_set(/xtitle) THEN o.xtitle = o->get(/firstdim) END ;; <_ensure_image_scale and help methods PRO spectobs::_ensure_image_scale ;; xtitle is an anyobs tag o = self o->private IF ~o->is_set(/image_scale) THEN BEGIN xtitle = strmid(o->get(/xtitle),0,4) CASE xtitle OF 'Inst': image_scale = [ptr_new(o->spectobs::_image_scale_x()), $ ptr_new(o->spectobs::_image_scale_y())] 'Sola': image_scale = [ptr_new(o->spectobs::_image_scale_x()), $ ptr_new(o->spectobs::_image_scale_y())] ; 'Time': image_scale = ptr_new(o->_image_scale_time_y()) 'Time': image_scale = [ptr_new(o->spectobs::_image_scale_time()), $ ptr_new(o->spectobs::_image_scale_y())] ; 'Wave': image_scale = ptr_new(o->_image_scale_wavelength_y()) 'Wave': image_scale = [ptr_new(o->spectobs::_image_scale_wavelength()),$ ptr_new(o->spectobs::_image_scale_y())] ELSE: message, 'Unknown xtitle' ENDCASE o.image_scale = ptr_new(image_scale) END END FUNCTION spectobs::_image_scale_x ;; spectobs::_ensure_image_iscale help method self->private dataobj = *(self->get(/data)) return, dataobj->getfovx()/dataobj->getnraster() END FUNCTION spectobs::_image_scale_y ;; spectobs::_ensure_image_scale help method self->private dataobj = *(self->get(/data)) return, dataobj->getfovy()/dataobj->getnaxis2() END FUNCTION spectobs::_image_scale_time ;; spectobs::_ensure_image_scale help method o = self o->private header = *(o->get(/header)) dataobj =*(o->get(/data)) DATE_END = gt_fkw(header,'DATE_END') DATE_OBS = gt_fkw(header,'DATE_OBS') nexp = dataobj->getnexp() nexp_prp = dataobj->getnexp_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 END FUNCTION spectobs::_image_scale_wavelength ;; spectobs::_ensure_image_scale help method o = self o->private header = *(o->get(/header)) nicons = o->get(/nicons) scale_wavelength = fltarr(nicons) scale_wavelength[*] = 1 dataobj = *(o->get(/data)) FOR i=0,nicons-1 DO BEGIN n = trim(i+1) wavemin = gt_fkw(header,'TWMIN'+n) wavemax = gt_fkw(header,'TWMAX'+n) nwavepix = dataobj->getxw(i) scale_wavelength[i] = $ [(wavemax-wavemin)/(nwavepix-1)] ENDFOR return, scale_wavelength END ;; PRO spectobs::_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 'Inst': image_origin = [ptr_new(o->spectobs::_image_origin_x()),$ ptr_new(o->spectobs::_image_origin_y())] 'Sola': image_origin = [ptr_new(o->spectobs::_image_origin_x()),$ ptr_new(o->spectobs::_image_origin_y())] 'Time': image_origin = [ptr_new(o->spectobs::_image_origin_time()),$ ptr_new(o->spectobs::_image_origin_y())] 'Wave': image_origin = [ptr_new(o->spectobs::_image_origin_wavelength()),$ ptr_new(o->spectobs::_image_origin_y())] ELSE: message, 'Unknown xtitle' ENDCASE o.image_origin = ptr_new(image_origin) END END FUNCTION spectobs::_image_origin_x self->private dataobj = *(self->get(/data)) xcen = dataobj->getxcen() fovx = dataobj->getfovx() return, xcen - fovx/2. END FUNCTION spectobs::_image_origin_y self->private dataobj = *(self->get(/data)) ycen = dataobj->getxcen() fovy = dataobj->getfovy() return, ycen - fovy/2. END FUNCTION spectobs::_image_origin_time o = self o->private nicons = o->get(/nicons) time = findgen(nicons) ;; Get the time and y coordinates from the fits header time[*] = *(self->get(/image_scale))[0]/2 ;t = *(self->get(/image_scale))[0]/2 ;time[*] = t[(where(finite(t)))] box_message,'CHECK IF IRIS DATA NEEDS DETECTOR OFFSET CORRECTION' return, time END FUNCTION spectobs::_image_origin_wavelength o = self o->private nicons = o->get(/nicons) ; origin_y = o->_image_origin_y() header = *(o->get(/header)) wavelength = fltarr(nicons) FOR i=0,nicons-1 DO $ wavelength[i] = [gt_fkw(header,'TWMIN'+trim(i+1))] ;; Correct for detector offset ;; wavelength_y = o->_correct_detector_offset(wavelength_y,/y) box_message,'CHECK IF IRIS DATA NEEDS DETECTOR OFFSET CORRECTION' return, wavelength END FUNCTION spectobs::_image_origin_wavelength_yFETTFJES 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]] ;; Correct for detector offset ;; wavelength_y = o->_correct_detector_offset(wavelength_y,/y) box_message,'CHECK IF IRIS DATA NEEDS DETECTOR OFFSET CORRECTION' return, wavelength_y END ;; ;; --------------------------------------------------------------------------- ;; ;; --------------------------------------------------------------------------- ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;; is_set(/data) THEN o->read_data returndata = (n_params() EQ 2) ? 1 : 0 o->_image_number_exist,i firstdimwave = o->get(/firstdim) EQ 'Wavelength (Å)' slit = o->get(/slit) ;; Decide what kind of observation we're looking at, call the appropriate ;; extract_image_xxxx method. ;; First check if it's a slit or slot observation IF ~slit THEN $ ;; It's a slot observation im = o->_extract_image_slot_raster(i) $ ELSE BEGIN ;; It's a slit observation IF ~o->emission_only(i) AND NOT firstdimwave THEN $ ;; It's a line with absorption features, e.g. Mg II h/k im = o->_extract_image_slit_absorption(i, ii=ii, _extra = extra) $ ELSE BEGIN ;; It's a clean emission line extraset = (size(extra,/tname) EQ 'UNDEFINED' ) ? 0 : 1 IF firstdimwave 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 BEGIN IF extraset AND NOT o->photosphere(i) THEN BEGIN ;; Fit gaussians, the output image is a I, v or w map depending ;; on the _extra keywords im = o->_extract_image_slit_cfit(i, ii = ii, gausslam = $ gausslam, _extra=extra, $ errorinfo=errorinfo) ENDIF ELSE BEGIN ;; Get a map of the integrated intensity IF extraset THEN BEGIN ;; This is a photospheric line and extra is set, i.e. int, ;; vel or width. The ideal thing to do would be to short ;; circuit and return -1, but when the SDC archive scripts ;; are running the code an intensity image has to be ;; produced. We therefore plot an integrated ;; intensity image when keyword INT is set and ii=0, short ;; circuit for VEL and WID and for ii gt 0. message,'No gaussian fitting to photospheric lines. Sorry.',/info intix = (where(strmid(tag_names(extra),0,2) EQ 'IN'))[0] intensity = (intix NE -1) ? extra.(intix) : 0 IF intensity AND ii EQ 0 THEN $ message,'Plotting integrated intensity image.',/info ELSE $ return, -1 ENDIF im = o->_extract_image_slit_integrated_intensity(i) ENDELSE ENDELSE ENDELSE ENDELSE IF arg_present(medianvalue) AND im[0] NE -1 THEN BEGIN missing = o->get(/missing) goodix = where(im GT missing, ngood) IF ngood GT 0 THEN medianvalue = median(im[goodix]) ELSE medianvalue = missing ENDIF IF returndata THEN data = im ELSE $ im = o->spectobs::_image_final_touch(im, errorinfo=errorinfo, maxvel=maxvel, $ errvelmax=errvelmax, $ gausslam = gausslam, $ photosphere = o->photosphere(i), $ _extra=extra) return,im END FUNCTION spectobs::_extract_image_slit_lambda, i o = self o->private dataobj = (*(o->get(/data))) dataobj->getwin,i,data,pos,/load sz = size(data) nexp_prp = dataobj->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 ;; Replace the instrument specific missing value with the anyobs global ;; missing value missing = dataobj->missing() missingix = where(im EQ missing, nmissing) IF nmissing NE -1 THEN im[missingix] = o->get(/missing) return, im END FUNCTION spectobs::_extract_image_slit_integrated_intensity, i o = self o->private dataobj = (*(o->get(/data))) dataobj->getwin,i,data,pos,/load missing = dataobj->missing() sz = size(data) meanprofile = o->_get_meanprofile(data, sz, edgeskip = 15) ;; Find pixelnumber in wavelength dimension where mean emission line peaks ; meanline = rebin(data,sz[1],1,1) ; ; ;; Skip the outer 15 pixels of the window to avoid edge effects (mostly due ; ;; to the missing pixels because of slit tilt correction, also significant ; ;; but unwanted contribution from wings of strong lines. ; ;; edgeskip = round(sz[1]*(5./100)) ; edgeskip = 15 ; dummy = max(meanline[edgeskip:sz[1]-edgeskip], mx,/nan) ; mx += edgeskip dummy = max(meanprofile, mx) integratedintensity = o->_integratedintensity(data, mx, missing) ;; Get the number of exposure times per raster position nexp_prp = dataobj->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 ;FOR EIS: ; rot = (o->get(/raster)) ? 1 : 4 rot = o->get_datacube_rotation() 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, im END FUNCTION spectobs::_integratedintensity, data, mx, missing self->private ;; We assume that we get a significant fraction of the line profile by ;; stepping 5 pixels before and 5 pixels after the peak of the mean ;; profile. This is a very very rough estimate, but for quicklook purposes ;; this seems to work. delta = 5 dmn = mx-delta > 0 ;; Don't go beyond 0! dmx = mx+delta < (size(data))[1]-1 ; Don't go beyond the array boundry integratedintensity = total(data[dmn:dmx, *,*], 1) ;; Sum from dmn to dmx ;; If we have missing data at specific pixels along the slit (i.e. not ;; taking into account missing exposures), the integratedintensity at such a ;; pixel will be missing*(dmx-dmn+1), i.e. missing*10 in most cases. We ;; identify such regions of missing data, and give them the spectobs ;; internal object missing value. Note - we have two different missing ;; values - the value of missing pixels in the fits files (-200 for IRIS) ;; and one that is used when line fitting or line feature recognition ;; fails (-1e3). In image_final_touch we need to exclude values that are ;; missing, and if the image we are looking at is velocity map, then a value ;; of -200 could be real, we therefore set the missing pixels to this ;; extreme value also in the case for integrated intensities. missix = where(integratedintensity EQ missing*(dmx-dmn+1),/ NULL) IF missix NE !NULL THEN integratedintensity[missix] = self->get(/missing) return, integratedintensity END FUNCTION spectobs::_extract_image_slot_raster, i ;; slit jaw images return,-1 END FUNCTION spectobs::emission_only, i ;; overridden in irisspectobs__define ;; As default all lines are emission lines. For lower transition region and ;; chromosphere some lines may have absorption features, e.g. the Mg II h/k ;; lines observed by IRIS. return,1 END FUNCTION spectobs::_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, medianvelcorr=medianvelcorr,$ quiet=quiet, ref__extra=extra ;;+ ;; 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 ;; kamio now in _ref_extra! ;; default, kamio, 1 NEEDED FOR EIS!!!!!!!!!!!!!!!!!!!!! default, newadef, 0 ; Get data and wavelength scale for window i dataobj = *(o->get(/data)) dataobj->getwin, i, data, pos,/load sz = size(data) ;; 4D if nexp_prp GT 1; replace the sz[4] dimensions with mean of all nexp_prp ;; exposures for that slit position nexp_prpGT1 = (sz[0] EQ 4) ? 1 : 0 IF nexp_prpGT1 THEN data = rebin(data, sz[1], sz[2], sz[3], 1) sz = size(data) lam = dataobj->getlam(i) lam = o->getcorrectedwavelengthcube(lam, sz, region=dataobj->getregion(i)) ;;;; T H I S W I L L G O I N T O E I S _ O B S ! ! ! ; ;; 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. ; extraset = (size(extra,/tname) EQ 'UNDEFINED' ) ? 0 : 1 ; IF keyword_set(kamio) THEN BEGIN ; lamcorr = eis_data->gethkwavecorr(median(lam), wvl_cube=lam) ; lam = lamcorr.cube ; lamcorrtime = lamcorr.time ; ENDIF ELSE BEGIN ; lamcorr = o->_gt_warren_lamcorr(lamcorrtime = lamcorrtime, $ ; lam_i = median(lam)) ; ; IF nexp_prpGT1 THEN 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 ; lam = lam3D ; ENDELSE ; data = transpose(data,[0,2,1]) ana = o->_gt_cfit_ana(data, lam, i, newadef=newadef, fittwice=fittwice, $ _extra = extra) 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(gaussind, result, lam, missing, $ intensity=intensity, velocity=velocity, width=width, $ lambda = lambda, lamcorrtime = lamcorrtime, $ medianvelcorr=medianvelcorr) missingix = where(im EQ missing,/NULL) IF missingix NE !NULL THEN im[missingix] = o->get(/missing) return,im END FUNCTION spectobs::_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) instrumentmissing = o->_get_instrument_specific_missing_value() anafilename = o->find_prepfile('ana',i = i, _extra=extra) ;; extra eq kamio=1 for eis IF keyword_set(newadef) OR ~file_test(anafilename) THEN BEGIN adeffile = o->find_prepfile('adef',i=i, _extra=extra) 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 = o->_get_meanprofile(data, sz, edgeskip = 15) ;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,instrumentmissing,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 spectobs::_get_instrument_specific_missing_value self->private return, self->get(/missing) END FUNCTION spectobs::_get_meanprofile, data, sz ;; overridden in irisspectobs ;; Returns the mean profile (summed over all exposures and along the slit) ;; of the input 3D data cube self->private return, total(total(data,2),2)/(sz[2]*sz[3]) END FUNCTION spectobs::_auto_adef_get_smoothing ;; Variable used by spectobs::_auto_adef. Default is no smoothing. self->private return,1 END FUNCTION spectobs::_auto_adef_get_meanmedian ;; Variable used by spectobs::_auto_adef. Default is 4.5 which is the value ;; found to be fitting for EIS observations. self->private return, 4.5 END FUNCTION spectobs::_auto_adef, lam, int ;; Automatically generate cfit analysis definitions based on input intensity and ;; wavelength arrays self->private sm = self->_auto_adef_get_smoothing() minmedian = self->_auto_adef_get_meanmedian() peakinds = gt_peaks(smooth(int,sm), fwhm=fwhm, minmedian=minmedian,/sort,/plot) IF peakinds[0] EQ -1 THEN BEGIN box_message,'Picking the middle pixel as the only "peak".' peakinds = n_elements(int)/2 ENDIF 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]) ;; Find out the forrect number of digits we want, for thislam lt 1000 we ;; need 3 + 2 (e.g. 195.12), for gt 1000 lt 10000 we need 4 + 2 ;; (e.g. 1335.80). thislam = lam0[i] pointpos = strpos(strtrim(thislam,2),'.') fortranformat = '(F'+trim(pointpos+3)+'.2)' lam0txt = trim(thislam, fortranformat) 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 spectobs::_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,7) 3: gausslam = strmid(fit.igauss3.name,9,7) 6: gausslam = strmid(fit.igauss4.name,9,7) 9: gausslam = strmid(fit.igauss5.name,9,7) 12: gausslam = strmid(fit.igauss6.name,9,7) 15: gausslam = strmid(fit.igauss7.name,9,7) 18: gausslam = strmid(fit.igauss8.name,9,7) 21: gausslam = strmid(fit.igauss9.name,9,7) 24: gausslam = strmid(fit.igauss10.name,9,7) 27: gausslam = strmid(fit.igauss11.name,9,7) 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 box_message,trim(gausslam) return, gausslam END FUNCTION spectobs::_extract_cfit_im, gaussind, result, lam, 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 spectobs::_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 spectobs::_image_final_touch, im, errorinfo=errorinfo, maxvel=maxvel, $ errvelmax=errvelmax, snratio=snratio, $ intensity=intensity, $ velocity=velocity, width=width, $ gausslam=gausslam, photosphere=photosphere,$ _extra=extra 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) extraset = (size(extra,/tname) EQ 'UNDEFINED' ) ? 0 : 1 IF extraset THEN self->set, dynamictitle=(tag_names(extra))[0] + ' intensity' $ ELSE self->set, dynamictitle='' self->set, rightytitle='' IF keyword_set(photosphere) THEN im = bytscl(im,/nan) ELSE im = 255-bytscl(im,/nan) 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, _extra = extra 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, _extra = extra 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, _extra = extra 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 IF keyword_set(photosphere) THEN im = bytscl(im,/nan) $ ELSE im = 255-bytscl(im,/nan) 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 spectobs::_erase_temp_variables self->private self.ii = -1 END FUNCTION spectobs::_extract_image_scale_slot, i self->private print,'fix spectobs::get_image_scale for slot observations!' return,[1,1] END FUNCTION spectobs::_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) image_scale_x = *image_scale[0] image_scale_y = *image_scale[1] ix = (n_elements(image_scale_x) EQ 1) ? 0 : i ;; will be arrays for iris, not for eis return,[image_scale_x[ix],image_scale_y[ix]] ENDELSE END FUNCTION spectobs::_extract_image_origin,i o = self o->private IF ~o->get(/slit) THEN return, o->_extract_image_scale_slot(i) $ ELSE BEGIN image_origin = o->get(/image_origin) image_origin_x = *image_origin[0] image_origin_y = *image_origin[1] ix = (n_elements(image_origin_x) EQ 1) ? 0 : i ;; will be arrays for iris, not for eis return,[image_origin_x[ix],image_origin_y[ix]] ENDELSE END PRO spectobs::_cenAndfov, xcen, ycen, fovx, fovy o = self o->private dataobj = *(o->get(/data)) xcen = dataobj->getxcen() ycen = dataobj->getycen() fovx = dataobj->getfovx() ;; SHOULD SUPPLY WINDOW NUMBER! fovy = dataobj->getfovy() ;; SHOULD SUPPLY WINDOW NUMBER! END FUNCTION spectobs::_extract_title,i o = self o->private dynamictitle = o->get(/dynamictitle);+o->_dynamictitle_addition() IF o->is_set(/ii) THEN title = dynamictitle ELSE BEGIN dynamictitle = (dynamictitle EQ '') ? '' : ' ('+dynamictitle+')' title = (*(o->get(/title)))[i] + dynamictitle ENDELSE return, title END ;; --------------------------------------------------------------------------- ;; private lamdiffmax = self->_get_lamdiffmax() IF keyword_set(gausslam) THEN BEGIN llist = self->_gt_commonblock_linelist() lam = llist.lam ion = llist.ion ind = gt_closest_index(gausslam, lam, mindiff=mindiff,/noprint) IF abs(mindiff) LT lamdiffmax 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.dynmictitle = trim(gausslam)+' Å ('+ion[ind]+')' ENDIF ELSE BEGIN llist_backup = self->_gt_commonblock_linelist(/backup) llist = llist_backup lam = llist.lam ion = llist.ion ind = gt_closest_index(gausslam, lam, mindiff=mindiff,/noprint) IF abs(mindiff) LT lamdiffmax THEN $ self.dynamictitle = ion[ind] + ' ' + trim(lam[ind]) + ' Å' ELSE $ self.dynamictitle = 'Fit at ' +gausslam+' Å' ENDELSE ENDIF END FUNCTION spectobs::_get_lamdiffmax ;; The maximum difference between the fitted wavelength of a line and the ;; rest wavelength in a linelist self->private return,0.05 END ;; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;; get_XXX functions for all spectobs structure tags. Each function checks if ;; header is a valid pointer, if not, call read_header which again ensures the ;; tag is set. ;; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FUNCTION spectobs::get_slit IF ~ptr_valid(self.header) THEN self->read_header return, self.slit END FUNCTION spectobs::get_raster IF ~ptr_valid(self.header) THEN self->read_header return, self.raster END ;; --------------------------------------------------------------------------- ;; get_XXX functions for all spectobs structure tags. ;; --------------------------------------------------------------------------- PRO spectobs::_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 ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; <_color_table and help methods> ;; ++++++++++++++++++++++++++++++++++++++++++++++++++++++ FUNCTION spectobs::_color_table_txt, i, velocity=velocity, width=width,$ altintensity=altintensity self->private ; wanted_ct_txt = (keyword_set(altintensity)) ? 'ALTINTENSITY' : 'INTENSITY' ;; We choose the colour table for intensity images... wanted_ct_txt = 'INTENSITY' ;; ...unless any of the keywords are set IF keyword_set(altintensity) THEN wanted_ct_txt = 'ALTINTENSITY' IF keyword_set(velocity) THEN wanted_ct_txt = 'VELOCITY' IF keyword_set(width) THEN wanted_ct_txt = 'WIDTH' return, wanted_ct_txt END PRO spectobs__define INT = 0 STR = '' d = {spectobs, $ slit:INT, $ raster:INT, $ nwin:INT, $ ii:INT, $ missing:INT, $ firstdim:STR, $ dynamictitle:STR $ } END