PRO idlines, specfile, linelist, xlo, xhi, xr, yr, atoms, thresh, rv, plotfile, spc_file, spc_broad, $ FIT = kft, EW = kew, NF = knf, ALL = kall, STIS = kst, PRINT = kpr, $ XTEND = kxt, XPAND = kxp, MULT = kmlt, DUPLEX = kdup, QUAD = kfour, QUIN = kfive, SEX = ksix, $ SPEC = kspc ;-------------------------------------------------------------------------- ; ; idlines ; ; IDL procedure to display ; a) a region of spectrum, together with a fit, if supplied, ; b) positions and relative strengths (gf or W_lambda) of lines in ; common atomic spectra may be marked in the lower panel. ; c) lines of sleceted atoms are identified in the top panel ; ; Command arguments: ; ; specfile - character string - name of file containing spectrum (and fit) ; assumed to be x,y file with 3-line header ; linelist - character string - name of file containing eq.widths ; assumed to be 'spectrum' output ; xlo - real variable - lower limit of plot ; xhi - real variable - upper limit of plot ; xr - real variable - range of individual plot ; yr - real variable - [ylo,yup] for each segment ; atoms - integer vector - atomic numbers of spectra to be identified ; [ thresh - real variable - threshold eq width or line strength to identify lines ] ; [ rv - real variable - radial velocity of spectrum ] ; [ plotfile - character string - name of file to which plot is to be written (default: 'idlines.eps') ] ; [ spc_file - character string - name of file containing comparison fit ] ; [ spc_broad - real variable - instrumental broadening to apply to comparison spectrum ] ; ;--------------------------------------------------------------------------- ann_idlines IF (N_PARAMS() LT 7 OR specfile EQ 'help') THEN BEGIN print, '****2018****' print, 'idlines: insufficient parameters' print, ' ' print, 'Usage:' print, ' ' print, 'idlines, specfile, linelist, xlo, xhi, xr, yr, ions, [thresh, [rv, ]], [keywords]' print, '' print, ' specfile - character string - name of file containing spectrum (and fit)' print, ' assumed to be x,y file with 3-line header' print, ' linelist - character string - name of file containing linelist' print, ' assumed to be in ''lte_lines'' format' print, ' xlo - real variable - lower limit of plot' print, ' xhi - real variable - upper limit of plot' print, ' xr - real variable - xrange for each segment ' print, ' yr - real array - [ylo,yup] for each segment ' print, ' atoms - integer vector - atomic numbers of spectra to be identified' print, '[thresh - real variable - threshold linestrength to identify lines]' print, '[rv - real variable - radial velocity of spectrum ]' print, '[plotfile - character string - name of file to which plot is to be written (default: idlines.eps ) ]' print, '[spc_file - character string - name of file containing comparison fit ]' print, '[spc_broad- real variable - instrumental broadening to apply to comparison spectrum ]' print, '' print, ' Keywords ' print, '/FIT - overplot the fit ' print, '/EW - read and select the lines using equivalent widths ' print, '/NF - read and select the lines using n*gf (central opacities) ' print, '/ALL - plot individual atomic spectra schematically ' print, '/MULT - add multiplet labels to line ids ' print, '/STIS - to allow for air/vacuum shift at 2000A ' print, '/DUPLEX - plot 2 ranges of spectrum on one panel ' print, '/QUAD - plot 4 ranges of spectrum on one panel ' print, '/XPAND - expand plot by factor 5 ' print, '/XTEND - extend line identification marks ' print, '/PRINT - plot hardcopy as a postscript file' print, '/SPECTRUM - overplot an additional spectrum' print, ' /YLO - real variable - lower y limit for each segment ' print, ' /YHI - real variable - upper y limit for each segment ' RETURN ENDIF ;; defaults for omitted parameters IF N_PARAMS() LT 8 THEN thresh = 0.0 IF N_PARAMS() LT 9 THEN rv = 0.0 IF N_PARAMS() LT 10 THEN plotfile= 'idlines.eps' IF N_PARAMS() LT 12 THEN spc_broad = 0.0 ; Check for keywords xft = 0 ; default IF KEYWORD_SET(kft) THEN xft = 1 xew = 0 ; default IF KEYWORD_SET(kew) THEN xew = 2 IF KEYWORD_SET(knf) THEN xew = 1 xall = 0 ; default IF KEYWORD_SET(kall) THEN xall = 1 xmlt = 0 ; default IF KEYWORD_SET(kmlt) THEN xmlt = 1 IF KEYWORD_SET(kmlt) THEN print,' Multiplet identification enabled: linelist must contain multiplet identifiers ' xstis = 0 IF KEYWORD_SET(kst) THEN xstis = 1 xxtd = 0 IF KEYWORD_SET(kxt) THEN xxtd = 1 xxpd = 0 IF KEYWORD_SET(kxp) THEN xxpd = 1 kone = 0 IF ( NOT KEYWORD_SET(kfour) AND NOT KEYWORD_SET (kfive) AND NOT KEYWORD_SET (ksix) ) THEN kone = 1 xspc = 0 IF KEYWORD_SET(kspc) THEN xspc = 1 ;; set y range defaults y1 = 0.45 y2 = 1.60 IF yr(0) GT 0. THEN y1 = yr(0) IF yr(1) GT 0. THEN y2 = yr(1) dev = 'X' xcp = 0 IF KEYWORD_SET(kpr) THEN BEGIN xcp = 1 dev='PS' set_plot,'ps' IF ( KEYWORD_SET(kone) ) THEN device,filename=plotfile,/color,/landscape,/encapsulated IF ( KEYWORD_SET(kfour) ) THEN device,filename=plotfile,/color,/portrait,xsize=20,ysize=20,yoffset=2,/encapsulated IF ( KEYWORD_SET(kfive) ) THEN device,filename=plotfile,/color,/portrait,xsize=20,ysize=25,yoffset=2,/encapsulated IF ( KEYWORD_SET(ksix) ) THEN device,filename=plotfile,/color,/portrait,xsize=20,ysize=25,yoffset=2,/encapsulated ENDIF idsize = 1.0 ; [default size for line identifiers] colors = GetColor(/Load) !p.multi=0 IF KEYWORD_SET(kdup) AND xall EQ 0 THEN !p.multi=[0,1,2,0,0] ;; get the spectrum to plot IF xft EQ 1 THEN BEGIN readspecfit, specfile, specnw, specwav, specflx, specfit ENDIF ELSE BEGIN readspec, specfile, specnw, specwav, specflx ENDELSE IF xxpd THEN specflx = specflx * 5 - 4 IF xxpd AND xft THEN specfit = specfit * 5 - 4 IF xstis THEN ww = specwav ;; shift spectrum to rest wavelength if velocity supplied IF ABS (rv) GT 1 THEN specwav = specwav * (1 - rv/2.997E5) ;; get the linelist with identifications IF (xew EQ 2) THEN $ readeqwid, linelist, listn, listwav, listlab, listion, listmlt, listgf, listew IF (xew EQ 1) THEN $ readlines, linelist, listn, listwav, listlab, listion, listmlt, listgf, listep, listfr IF (xew EQ 0) THEN $ readlte, linelist, listn, listwav, listlab, listion, listmlt, listgf, listep, listfr nx = fix((xhi - xlo) / xr) dx = xr cont = ' ' ;; IF xstis THEN spec0 = specwav ;; get the additional model if defined IF (xspc EQ 1) THEN BEGIN sp2rd, spc_file, spc_wav, spc_flx, spc_head IF (spc_broad GT 0) THEN $ spc_flx = qsmooth( spc_wav, spc_flx, spc_broad ) ENDIF ;; Trim the linelist to match the spectrum wok = where(specwav gt 0, nwok) wmin = min(specwav(wok)) wmax = max(specwav(wok)) lok = where( (listwav ge wmin) and (listwav le wmax), nlok ) listn = listn(lok) listwav=listwav(lok) listlab=listlab(lok) listion=listion(lok) listmlt=listmlt(lok) listgf =listgf(lok) IF (xew EQ 2) THEN listew =listew(lok) IF (xew LT 2) THEN listep =listep(lok) IF (xew LT 2) THEN listfr =listfr(lok) ;; ============================= ;; Loop over wavelength intervals FOR i = 0,nx-1 DO BEGIN IF (cont NE 'N' AND cont NE 'n' and cont NE 'q' AND cont NE 'Q') THEN BEGIN x1 = xlo + i * dx - dx/25. x2 = xlo + i * dx + dx + dx/25. IF xall THEN BEGIN y1 = -0.8 y2 = 1.8 ENDIF ;; plot the spectrum and fit ;; because linelist is in air for GE 2000 A ;; shift the observed spectrum to air so that ;; when identifying features from linelist do ;; not need to convert values in linelist IF xstis THEN BEGIN w2000 = where(ww GE 2000, nw2000) IF nw2000 GT 1 THEN specwav(w2000) = specwav(w2000)/1.0003071 ENDIF ; establish the box !x.title='!4k!3 (!3' + STRING(197B) + '!X)' !y.title='A!D!4k!3!N' IF KEYWORD_SET(kone) THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.3,charthick=1.0, /nodata plotpos = fltarr(4,4) plotpos(*,3) = [0.1,0.06,0.99,0.255] plotpos(*,2) = [0.1,0.305,0.99,0.50] plotpos(*,1) = [0.1,0.55,0.99,0.745] plotpos(*,0) = [0.1,0.795,0.99,0.99] IF KEYWORD_SET(kfour) AND xall EQ 0 THEN BEGIN imod = i mod 4 IF imod EQ 0 THEN erase IF imod EQ 0 THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod),xtitle=' ',xticklen=0.05,yticklen=0.02 IF ( imod EQ 1 OR imod EQ 2) THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod), /noerase, xtitle=' ',xticklen=0.05,yticklen=0.02 IF imod EQ 3 THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod), /noerase,xticklen=0.05,yticklen=0.02 ENDIF plotpos = fltarr(4,5) plotpos(*,4) = [0.04,0.03,0.97,0.20] plotpos(*,3) = [0.04,0.23,0.97,0.40] plotpos(*,2) = [0.04,0.43,0.97,0.60] plotpos(*,1) = [0.04,0.63,0.97,0.80] plotpos(*,0) = [0.04,0.83,0.97,1.00] plotpos(1,*) = (plotpos(1,*)+0.022)/1.043 plotpos(3,*) = (plotpos(3,*)+0.022)/1.043 IF KEYWORD_SET(kfive) AND xall EQ 0 THEN BEGIN imod = i mod 5 IF imod EQ 0 THEN erase IF imod EQ 0 THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod),xtitle=' ',xticklen=0.05,yticklen=0.02 IF ( imod GE 1 AND imod LE 3) THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod), /noerase, xtitle=' ',xticklen=0.05,yticklen=0.02 IF imod EQ 4 THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod), /noerase,xticklen=0.05,yticklen=0.02 ENDIF plotpos = fltarr(4,6) plotpos(*,5) = [0.04,0.03,0.98,0.18] plotpos(*,4) = [0.04,0.21,0.98,0.34] plotpos(*,3) = [0.04,0.37,0.98,0.50] plotpos(*,2) = [0.04,0.53,0.98,0.66] plotpos(*,1) = [0.04,0.69,0.98,0.82] plotpos(*,0) = [0.04,0.85,0.98,0.98] IF KEYWORD_SET(ksix) AND xall EQ 0 THEN BEGIN imod = i mod 6 IF imod EQ 0 THEN erase IF imod EQ 0 THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod),xtitle=' ',xticklen=0.05,yticklen=0.02 IF ( imod GE 1 AND imod LE 4) THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod), /noerase, xtitle=' ',xticklen=0.05,yticklen=0.02 IF imod EQ 5 THEN plot, [x1,x2],[y1,y1], $ xrange=[x1,x2], yrange=[y1,y2],xstyle=1, ystyle=1,charsize=1.0,charthick=1.0, $ /nodata , pos = plotpos(*,imod), /noerase,xticklen=0.05,yticklen=0.02 ENDIF ; plot the spectrum and the continuum oplot, specwav, specflx, psym=10, thick=2 ; plot the spectrum oplot, [xlo,xhi], [1,1], linestyle = 1 ; mark the continuum IF xstis THEN specwav = spec0 IF xft THEN $ oplot, specwav, specfit, color=colors.red,thick=2 ; overplot the fit IF xft THEN $ print, 'overplotting current best fit in red' IF xspc THEN $ oplot, spc_wav, spc_flx, color=colors.blue,thick=2 ; overplot the comparison IF xspc THEN $ print, 'overplotting comparison spectrum in blue' ;; mark individual lines ;determine lines that are greater than the threshold ;this is done element at a time wavt=fltarr(1) labt=strarr(1) mult=strarr(1) wavt2=fltarr(1) labt2=strarr(1) mult2=strarr(1) ik1=0 ik2=0 FOR zel = 0,n_elements(atoms)-1 DO BEGIN ;DETERMINE IONS GREATER THAN THRESH IF xew EQ 2 THEN $ markion, atoms(zel), x1, x2, listn, listwav, listion, listlab, listmlt, listew, -thresh, wl, lab, mlt IF xew EQ 1 THEN $ markion, atoms(zel), x1, x2, listn, listwav, listion, listlab, listmlt, listfr, thresh, wl, lab, mlt IF xew EQ 0 THEN $ markion, atoms(zel), x1, x2, listn, listwav, listion, listlab, listmlt, listfr, 0.0, wl, lab, mlt ;MAKE ARRAY OF ELEMENTS TO MARK if atoms(zel) lt 21 then begin wavt = [wavt,wl] labt = [labt,lab] mult = [mult,mlt] endif if atoms(zel) gt 20 then begin wavt2 = [wavt2,wl] labt2 = [labt2,lab] mult2 = [mult2,mlt] endif ENDFOR nwt=n_elements(wavt) nw2=n_elements(wavt2) ; print,'nlines:',nwt,nw2 ; debugger ;sort lines that met the thresh criteria into wavelength order IF (nwt GT 1) then BEGIN xx=sort(wavt) wavt=wavt(xx) labt=labt(xx) mult=mult(xx) ENDIF IF (nw2 GT 1) then BEGIN xx=sort(wavt2) wavt2=wavt2(xx) labt2=labt2(xx) mult2=mult2(xx) ENDIF ;TIDY UP ARRAYS (REMOVE ALL 0 VALUED ELEMENTS) iwt = where(wavt2 GT 0, nwt) IF (nwt GT 0) then BEGIN pwav2=wavt2(iwt) plab2=strcompress(labt2(iwt),/remove_all) pmlt2=mult2(iwt) ;OPLOT LINE IDENTIFICATIONS Z > 20 (BLUE) IF xmlt THEN plab2 = plab2+pmlt2 idls_oplot,specwav,specflx,specfit, x1,x2,y1,y2,pwav2,plab2,extend=xxtd,lcharsize=idsize ENDIF ;TIDY UP ARRAYS (REMOVE ALL 0 VALUED ELEMENTS) iwt = where(wavt GT 0, nwt) IF (nwt GT 0) then BEGIN pwav=wavt(iwt) plab=labt(iwt) plab=strcompress(labt(iwt),/remove_all) pmlt=mult(iwt) ;PLOT LINE IDENTIFICATIONS Z < 21 (RED) IF xmlt THEN plab = plab+pmlt idls_plot,specwav,specflx,specfit, x1,x2,y1,y2,pwav,plab,extend=xxtd,lcharsize=idsize ENDIF ;; mark separate spectra IF xall THEN BEGIN FOR zel = 2,28 DO BEGIN IF xew EQ 2 THEN $ markeqwid, zel, x1, x2, listn, listwav, listion, listlab, listmlt, $ listgf, listew IF xew EQ 1 THEN $ markeqwid, zel, x1, x2, listn, listwav, listion, listlab, listmlt, $ listgf, listfr ENDFOR ENDIF ;; expand the plot ; oplot, [x1, x2], [0.4,0.4] ; oplot, [x1, x2], [1.6,1.6], linestyle=1 ; oplot, specwav, 5 * specflx - 5 + 1.6, clip = [x1,1.1,x2,1.7] ; IF fit EQ 1 THEN $ ; oplot, specwav, 5 * specfit - 5 + 1.6, clip = [x1,1.1,x2,1.7], color=colors.red, thick=2 ;; check for continuation IF i LT nx-1 AND NOT xcp THEN $ read,prompt='Continue (-/n/q) ? >',format='(a)', cont ENDIF ENDFOR IF xcp EQ 1 THEN BEGIN print,' ' print,' Output: idlines plot written to file: ',plotfile print,' ' device,/close set_plot,'X' dev='X' ENDIF END