;+ ; NAME: ; AUX_PLONG ; PURPOSE: ; Plot age/metallicity profiles for long-slit spectra ; USAGE: ; aux_plong, , PS=ps, /NOLOSVD, /NOPOP, CHI2=chi2, HOR_AGE=hor_age, HOR_MET=hor_met, VER=ver, /NOYTIT, /NYLAB, /LOGF, XTIT=xtit, CODK=codk, IRAD=irad, GALNAME=galname, XRANGE=xrange, CZRANGE=czrange, SIGRANGE=sigrange, AGERANGE=agerange, METRANGE=metrange ; ; DESCRIPTION: ; Plot the kin+SSP profile generated with the program AUX_ALONG ; Each point is a separate '.res' file named as: '_r''.res' ; where is the generic name of the analysed frame, and ; is the position along the slit, in arcsec. ; ; ARGUMENTS: ; - input ; The prefix denoting the name (path+name) of the ; analysed frame; indentical to in AUX_ALONG ; ; KEYWORDS: ; PS ; Save the plot in a postscript file named ; /NOLOSVD ; Do not plot the kinematics ; /NOPOP ; Do not plot the age/metallicity profiles ; /CHI2 ; Plot the chi2 values ; HOR_AGE [Gyr] ; Overplot age = [,]. For example the ; mean age at 1 arcsec and 1 effective radius ; HOR_MET [dex] ; Same as hor_age but for metallicity ; VER [dimensionless] ; Overplot as a vertical line the effective radius of the galaxy, real ; /NOYTIT ; Do not plot the title of the y-axes ; /NYLAB ; Do not plot the labels of the y-axes ; /LOGF ; Plot the flux in log ; XTIT ; x-axis title, string ; CODK [dimensionless] ; Coding of the symbols for each individual profile. Useful if ; few data sets are overplotted ; 0 : red/blue ; 1 : magenta/cyan ; 2 : 'black' ; 3 : continuous line ; IRAD ; GALNAME ; Name of the galaxy, to be overplotted in the right upper ; corner, string ; XRANGE default = [0, abs(max())] ; array with the x-axis plotting range ; CZRANGE [km/s] ; Plotting range for cz ; SIGRANGE [km/s] ; Velocity dispersion range ; AGERANGE [Gyr] ; Age range ; METRANGE [dex] ; metallicity range ; ; HISTORY: ; MK: 2009/09 ; Creation, used for 2009MNRAS.396.2133K ; PhP: 2010/05/26 ; Modifications for 2011MNRAS.417.1643K ; MK: 2012/10/26 ; Documentation ;- ; CATEGORY: ; AUX ;; Subrouting for plot_radprof, reading the profiles function long_read, prefix nn = 0 for k=0, n_elements(prefix)-1 do begin len = strlen(prefix[k]) + 2 if k eq 0 then begin name = file_search(prefix[k]+ '_r*.res', COUNT=nf) if nf eq 0 then begin ; search the 'old' format (without 'r') name = file_search(prefix+ '_*.res', COUNT=nf) len -- endif if nf gt 0 then pnum = intarr(nf) endif else begin name = [name, file_search(prefix[k]+ '_r*.res', COUNT=nf)] if nf eq 0 then begin ; search the 'old' format (without 'r') name = [name, file_search(prefix[k]+ '_*.res', COUNT=nf)] len -- endif if nf gt 0 then pnum = [pnum, intarr(nf)+k] endelse if nf gt 0 then rr = fltarr(nf) for i = 0, nf-1 do begin ll = strlen(strmid(name[nn+i], len)) - 4 rr[i] = strmid(name[nn+i], len, ll) endfor if nf gt 0 then if n_elements(rad) eq 0 then rad = rr else rad = [rad, rr] nn += nf endfor nf = n_elements(pnum) if nf gt 0 then begin print, name[0] chisq = replicate(0d,nf) cz = replicate(0d,nf) & e_cz = replicate(0d,nf) sig = replicate(0d,nf) & e_sig = replicate(0d,nf) age = replicate(0d,nf) & e_age = replicate(0d,nf) met = replicate(0d,nf) & e_met = replicate(0d,nf) for i=0,nf-1 do begin str = uly_solut_tread(name[i], STATUS=status) ; at this point we should test if str[0].cmp[0].para is pointing to data if size(str, /TYPE) eq 8 then begin chisq[i] = str[0].chisq cz[i] = str[0].losvd[0] & sig[i] = str[0].losvd[1] e_cz[i] = str[0].e_losvd[0] & e_sig[i] = str[0].e_losvd[1] age[i] = exp((*str[0].cmp[0].para)[0].value) & met[i] = (*str[0].cmp[0].para)[1].value e_age[i] = age[i]*(*str[0].cmp[0].para)[0].error & e_met[i] = (*str[0].cmp[0].para)[1].error endif else rad[i] = !VALUES.F_NAN undefine, str endfor ; print, 'Output file: ', prefix+'.txt' ; openw, lun, prefix+'.txt', /GET_LUN ; printf, lun, "# Radius, arcsec chi2 vel, km/s +/- sigma, km/s" ; !TEXTUNIT = LUN ; forprint, rad, chisq, cz, e_cz, sig, e_sig, text=5, FORMAT='(f10.3,f10.2, 4f6.2)' ; close, lun ; free_lun, lun return, {p:pnum, rad:rad, chisq:chisq, cz:cz, e_cz:e_cz, sig:sig, e_sig:e_sig, age:age, e_age:e_age, met:met, e_met:e_met} endif else message, 'No files ' + prefix[0]+ '_r*.res could be found' return, 0 end ;--------------------------------------------------------------------------------------- pro aux_plong, nm, PS=ps, NOLOSVD=nolosvd, NOPOP=nopop, HOR_AGE=hor_age, HOR_MET=hor_met, VER=ver, NOYTIT=noytit, NYLAB=noylab, LOGF=logf, XTIT=xtit, CODK=codk, IRAD=irad, GALNAME=galname, XRANGE=xrange, CZRANGE=czrange, SIGRANGE=sigrange, AGERANGE=agerange, METRANGE=metrange, CHI2=chi2 if n_elements(nm) eq 0 then $ print, 'Usage: AUX_PLONG, name ...' ; ; ICZ : give the possibility to revert the radial axis for some profiles profile = long_read(nm) ;profile = long_read(nm+'s') ;if size(profile, /TYPE) ne 8 then profile = long_read(nm) rk = sort(profile.rad) p = profile.p[rk] + 1 rad = profile.rad[rk] chisq = profile.chisq[rk] cz = profile.cz[rk] e_cz = profile.e_cz[rk] sig = profile.sig[rk] e_sig = profile.e_sig[rk] age = profile.age[rk] e_age = profile.e_age[rk] met = profile.met[rk] e_met = profile.e_met[rk] nprof = n_elements(nm) if n_elements(irad) ge nprof then begin for k=0,nprof-1 do begin nn = where(p eq k+1, nc) if nc gt 0 then rad[nn] *= irad[k] endfor endif if n_elements(galname) eq 0 then galname = nm[0] if n_elements(xrange) eq 0 then xrange = [0, 20] psymsize = 0.7 positive = where(rad ge 0e and abs(rad) ge xrange[0] and abs(rad) le xrange[1]) negative = where(rad lt 0e and abs(rad) ge xrange[0] and abs(rad) le xrange[1]) posik = ptrarr(nprof, /ALLO) cnt_posik = intarr(nprof+2) negak = ptrarr(nprof, /ALLO) cnt_negak = intarr(nprof+2) for k = 0, nprof-1 do begin *posik[k] = where(rad ge 0e and p eq k+1, cnt) cnt_posik[k] = cnt *negak[k] = where(rad lt 0e and p eq k+1, cnt) cnt_negak[k] = cnt endfor age /= 1000. e_age /= 1000d ; Mask some bad points (hum be sure they are really bad!) bad = where(sig gt 400 or e_age eq 0, cnt, COMPL=good) if cnt gt 0 then begin cz[bad] = !VALUES.F_NAN e_cz[bad] = !VALUES.F_NAN sig[bad] = !VALUES.F_NAN e_sig[bad] = !VALUES.F_NAN age[bad] = !VALUES.F_NAN e_age[bad] = !VALUES.F_NAN met[bad] = !VALUES.F_NAN e_met[bad] = !VALUES.F_NAN endif name = file_search(nm[0]+ '_profile.fits', COUNT=nf) if nf gt 0 then begin ; there is a light profile fits_read, name, profile, hdr naxis2 = sxpar(hdr, 'NAXIS1', COUNT=count) cdelt2 = double(string(sxpar(hdr, 'CDELT2', COUNT=count))) if count eq 0 then $ cdelt2 = double(string(sxpar(hdr, 'CD2_2', COUNT=count))) if count eq 0 then cdelt2 = 1 ; default = 1 arcsec/pix peak = double(string(sxpar(hdr,'CRPIX1'))) xx = (findgen(naxis2) - peak) * cdelt2 minrad = min(double(rad), MAX=maxrad) xxx = xx[where(xx ge minrad and xx le maxrad)] profile /= max(profile[peak]) yyy = profile[where(xx ge minrad and xx le maxrad)] positive1 = where(xxx ge 0e, COMPLEMENT=negative1) if size(negative1, /N_DIM) eq 0 then negative1=0 print, 'radial range', minrad, maxrad, ' NAXIS1:', strtrim(naxis2,2), ' CDELT2:', strtrim(cdelt2,2) endif ; Read Re if n_elements(ver) eq 0 then begin nmr = nm[0]+'_1Re.fits' if file_test(nmr) eq 1 then begin fits_read, nmr, data, hdr ver = sxpar(hdr, 'EXTRAD') ; Read Re print, 'Re=', ver, ' arcsec' if ver le xrange[0] or ver ge xrange[1] then undefine, ver endif endif ; Read the SSP solutions in 1arsec and in 1 Re if n_elements(hor_met) eq 0 and n_elements(hor_age) eq 0 then begin file = nm[0]+'_1arcsec.res' if file_test(file) eq 0 then file = nm[0]+'_1arcsec.res' if file_test(file) eq 1 then begin str = uly_solut_tread(file) hor_age = exp((*str[0].cmp[0].para)[0].value)/1000. hor_met = (*str[0].cmp[0].para)[1].value hor_sig = str[0].losvd[1] heap_free, str file = nm[0]+'_1Re.res' if file_test(file) eq 0 then file = nm[0]+'_1Re.res' if file_test(file) eq 1 then begin str = uly_solut_tread(file) hor_age = [hor_age,exp((*str[0].cmp[0].para)[0].value)/1000.] hor_met = [hor_met,(*str[0].cmp[0].para)[1].value] hor_sig = [hor_sig, str[0].losvd[1]] heap_free, str endif endif endif npanel = 5 if keyword_set(nolosvd) then npanel -= 2 if keyword_set(nopop) then npanel -= 2 if n_elements(xxx) eq 0 then npanel -- if keyword_set(chi2) then npanel ++ !P.MULTI = [0, 1, npanel] redk = intarr(nprof) + fsc_color('White') bluk = intarr(nprof) + fsc_color('White') if n_elements(codk) eq 0 then begin codk = indgen(nprof) if nprof gt 2 then codk[2:*] = 2 endif for k=0,nprof-1 do begin case codk[k] of 0 : begin redk[k] = fsc_color('Red') bluk[k] = fsc_color('Blue') end 1 : begin redk[k] = fsc_color('violet') bluk[k] = fsc_color('dodgerblue') end else : begin redk[k] = fsc_color('black') bluk[k] = fsc_color('black') endelse endcase endfor if keyword_set(ps) then begin set_plot, 'ps' !p.charsize=4 !p.thick=2 !x.thick=3 !y.thick=3 !P.FONT = 0 !X.TICKLEN = 0.07 !Y.TICKLEN = 0.03 !X.MARGIN = [7., 2.] ; def [10., 3.] marg between axes and edge of plot region !Y.MARGIN = [10., 1.] xs = 10 & ys = 12 device, FILENAME=ps+'.eps', /COLOR, XS=xs, YS=ys, FONT_SIZE=6, /TIMES, /ENC endif ; if n_elements(xtit) eq 1 and not keyword_set(chi2) then xtit='R, arcsec' keep = where(rad ge -xrange[1] and rad le xrange[1]) if n_elements(xxx) gt 0 then begin ; if there is a light profile if keyword_set(flog) then begin yyy = alog10(yyy) rp = where(abs(xxx) ge 2 and abs(xxx) le 5) cp = poly_fit(abs(xxx[rp]), yyy[rp], 1) yyp = cp[0] + cp[1] * abs(xxx) yyy -= yyp yrange = [-0.1, 0.4] endif else yrange = minmax(yyy) delta = yrange[1]-yrange[0] yrange[0] -= 0.1*delta yrange[1] += 0.1*delta if keyword_set(noylab) then ytickname=replicate(' ', 30) $ else begin if keyword_set(flog) then begin if not keyword_set(noytit) then ytit='log(Flux)-exp' else undefine, ytit if not keyword_set(noylab) then ytickformat='(f4.1)' endif else begin if not keyword_set(noytit) then ytit='Flux' else undefine, ytit if not keyword_set(noylab) then ytickformat='(f4.1)' endelse endelse plot, xxx[positive1], yyy[positive1], YMARGIN=[0., 1.], $ XR=xrange, XTICKNAME=replicate(' ', 30), $ XST=3, YST=3, YR=yrange, $ YTICKNAME=ytickname, YTICKFORMAT=ytickformat, /NODATA, $ BACK=fsc_color('white'), COLOR=fsc_color('black') ytpos= float(!D.y_ch_size)/!D.x_size * 2 ; alignment of the Y-titles if n_elements(ytit) gt 0 then $ xyouts, ytpos, mean(!y.window), ytit, CHARSIZE=!p.charsize*0.5, /NORMAL, ALIGN=0.5, ORIENT=90 plotsym, 0, .6, THICK=1, /FILL oplot, xxx[positive1], yyy[positive1], PSYM=8, COLOR=fsc_color('Red') plotsym, 0, .6, THICK=1, /FILL oplot, [-1.e*xxx[negative1]], [yyy[negative1]], PSYM=8, COLOR=fsc_color('Blue') if keyword_set(ver) then begin plots, [ver, ver], [yrange], LINEST=2, COLOR=fsc_color('Royal Blue') endif xyouts, 0.91, 0.94, galname[0], CHARSIZE=!p.charsize*0.5, /NORMAL, ALIGN=1 endif if not keyword_set(nolosvd) then begin if size(negative, /N_DIM) eq 0 then negative=0 if n_elements(czrange) eq 0 then begin yrange = minmax([cz[positive],-cz[negative]],/NAN) delta = yrange[1]-yrange[0] yrange[0] -= 0.1*delta yrange[1] += 0.1*delta endif else yrange=czrange if not keyword_set(noytit) then ytit='v, km/s' else undefine, ytit if not keyword_set(noylab) then ytickformat='(f4.1)' YMARGIN= (keyword_set(nopop))?[2, 0.]:[0.,0.] plot, rad[positive], cz[positive], YMARGIN=ymargin, XR=xrange, XTICKNAME=replicate(' ', 30), XST=3, /NODATA, YST=3, YR=yrange, YTICKNAME=ytickname, BACK=fsc_color('white'), COLOR=fsc_color('black') if n_elements(ytit) gt 0 then $ xyouts, ytpos, mean(!y.window), ytit, CHARSIZE=!p.charsize*0.5, /NORMAL, ALIGN=0.5, ORIENT=90 for k=nprof-1,0,-1 do begin if cnt_posik[k] gt 0 then begin plotsym, 0, psymsize, THICK=2, /FILL rk = where(finite(cz[*posik[k]]), cnt) if cnt gt 0 then begin rk = (*posik[k])[rk] if codk[k] lt 3 then begin oploterror, rad[rk], cz[rk], e_cz[rk], PSYM=8, ERRCOLOR=redk[k], COLOR=redk[k] endif else begin ; rk = (*posik[k])[sort(rad[*posik[k]])] oplot, rad[rk], cz[rk], PSYM=0, THICK=2.0, LINES=0 endelse endif endif endfor for k=nprof-1,0,-1 do begin if cnt_negak[k] gt 0 then begin plotsym, 0, psymsize, THICK=2, /FILL rk = where(finite(cz[*negak[k]]), cnt) if cnt gt 0 then begin rk = (*negak[k])[rk] if codk[k] lt 3 then begin oploterror, -1e*rad[rk], cz[rk], e_cz[rk], PSYM=8, ERRCOLOR=bluk[k], COLOR=bluk[k] endif endif endif endfor if keyword_set(ver) then begin plots, [ver, ver], [yrange], LINEST=2, COLOR=fsc_color('Royal Blue') endif if n_elements(sigrange) eq 0 then begin yrange = minmax(sig[keep], /NAN) delta = yrange[1]-yrange[0] yrange[0] -= 0.1*delta yrange[1] += 0.1*delta endif else yrange=sigrange if not keyword_set(noytit) then ytit=TeXtoIDL('\sigma')+', km/s' else undefine, ytit if not keyword_set(noylab) then ytickformat='(f4.1)' if keyword_set(nopop) then $ plot, rad[positive], sig[positive], YMARGIN=[3,-2], XR=xrange, XST=3, XTIT=xtit, /NODATA, YST=3, YR=yrange, YTICKNAME=ytickname, BACK=fsc_color('white'), COLOR=fsc_color('black') $ else $ plot, rad[positive], sig[positive], YMARGIN=[0.,0.], XTICKNAME=replicate(' ', 30), XR=xrange, XST=3, /NODATA, YST=3, YR=yrange, YTICKNAME=ytickname, BACK=fsc_color('white'), COLOR=fsc_color('black') if n_elements(ytit) gt 0 then $ xyouts, ytpos, mean(!y.window), ytit, CHARSIZE=!p.charsize*0.5, /NORMAL, ALIGN=0.5, ORIENT=90 for k=nprof-1,0,-1 do begin if cnt_posik[k] gt 0 then begin plotsym, 0, psymsize, THICK=2, /FILL rk = where(finite(sig[*posik[k]]), cnt) if cnt gt 0 then begin rk = (*posik[k])[rk] if codk[k] lt 3 then begin oploterror, rad[rk], sig[rk], e_sig[rk], PSYM=8, ERRCOLOR=redk[k], COLOR=redk[k] endif else begin ; rk = (*posik[k])[sort(rad[*posik[k]])] oplot, rad[rk], sig[rk], PSYM=0, THICK=2.0, LINES=0 endelse endif endif endfor for k=nprof-1,0,-1 do begin if cnt_negak[k] gt 0 then begin plotsym, 0, psymsize, THICK=2, /FILL rk = where(finite(cz[*negak[k]]), cnt) if cnt gt 0 then begin rk = (*negak[k])[rk] if codk[k] lt 3 then begin oploterror, -1e*rad[rk], sig[rk], e_sig[rk], PSYM=8, ERRCOLOR=bluk[k], COLOR=bluk[k] endif endif endif endfor if keyword_set(ver) then begin plots, [ver, ver], [yrange], LINEST=2, COLOR=fsc_color('Royal Blue') endif if n_elements(hor_sig) gt 0 then begin ; if n_elements(hor_sig) gt 1 then $ ; plots, [5,xrange[1]], [hor_sig[1],hor_sig[1]], color=fsc_color('Green') ; plots, [0, 5], [hor_sig[0], hor_sig[0]], color=fsc_color('Green') endif endif if not keyword_set(nopop) then begin ;yrange = minmax(age[keep]) ;delta = yrange[1]-yrange[0] ;yrange[0] -= 0.1*delta ;yrange[1] += 0.1*delta yrange=[0.9,20.] if n_elements(agerange) gt 0 then yrange = agerange if not keyword_set(noytit) then ytit='age, Gyr' else undefine, ytit if not keyword_set(noylab) then ytickformat='(f4.1)' ytickv = tickslog(yrange, NMIN=3) yticks = n_elements(ytickv)-1 plot, rad[positive], age[positive], YMARGIN=[2., 0.], $ XTICKNAME=replicate(' ', 30), XR=xrange, XST=3, $ /NODATA, YST=3, YR=yrange , YTICKNAME=ytickname, $ YTICKFORMAT=ytickformat, /YLOG, YTICKS=yticks, YTICKV=ytickv, $ BACK=fsc_color('white'), COLOR=fsc_color('black') if n_elements(ytit) gt 0 then $ xyouts, ytpos, mean(!y.window), ytit, CHARSIZE=!p.charsize*0.5, /NORMAL, ALIGN=0.5, ORIENT=90 if keyword_set(ver) then begin plots, [ver, ver], [yrange], LINEST=2, COLOR=fsc_color('Royal Blue') endif for k=nprof-1,0,-1 do begin if cnt_posik[k] gt 0 then begin plotsym, 0, psymsize, THICK=2, /FILL if codk[k] lt 3 then begin oploterror, rad[*posik[k]], age[*posik[k]], e_age[*posik[k]], PSYM=8, ERRCOLOR=redk[k], COLOR=redk[k] endif else begin rk = (*posik[k])[sort(rad[*posik[k]])] ;; to make the polyfill we need to sort the points and to cut the ;; things outside the x-y ranges, otherwise it plots outside the graph sortrad = sort(rad[rk]) rad[rk] = abs(rad[rk[sortrad]]) & age[rk] = age[rk[sortrad]] & e_age[rk] = 0.5 agemin = age[rk] - e_age[rk] agelim = where(agemin lt yrange[0], cnt) if cnt gt 0 then agemin[agelim] = yrange[0] agemax = age[rk] + e_age[rk] agelim = where(agemax gt yrange[1], cnt) if cnt gt 0 then agemax[agelim] = yrange[1] ; polyfill,/LINE_FIL, [rad[rk], reverse(rad[rk])], [agemin, reverse(agemax)], COLOR=fsc_color('dark gray') oplot, rad[rk], age[rk], THICK=2, LINES=0 endelse endif endfor for k=nprof-1,0,-1 do begin if cnt_negak[k] gt 0 then begin plotsym, 0, psymsize, THICK=2, /FILL if codk[k] lt 3 then begin oploterror, -1e*rad[*negak[k]], age[*negak[k]], e_age[*negak[k]], PSYM=8, ERRCOLOR=bluk[k], COLOR=bluk[k] endif else begin rk = (*negak[k])[sort(rad[*negak[k]])] ;; to make the polyfill we need to sort the points and to cut the ;; things outside the x-y ranges, otherwise it plots outside the graph sortrad = sort(rad[rk]) rad[rk] = abs(rad[rk[sortrad]]) & age[rk] = age[rk[sortrad]] & e_age[rk] = 0.5 agemin = age[rk] - e_age[rk] agelim = where(agemin lt yrange[0], cnt) if cnt gt 0 then agemin[agelim] = yrange[0] agemax = age[rk] + e_age[rk] agelim = where(agemax gt yrange[1], cnt) if cnt gt 0 then agemax[agelim] = yrange[1] ; polyfill,/LINE_FIL, [rad[rk], reverse(rad[rk])], [agemin, reverse(agemax)], COLOR=fsc_color('dark gray') oplot, rad[rk], age[rk], THICK=2, LINES=2 ; oplot, -rad[rk], age[rk], PSYM=0, THICK=2.0, LINEST=2 endelse endif endfor if n_elements(hor_age) gt 0 then begin if n_elements(hor_age) gt 1 and $ hor_age[1] gt yrange[0] and hor_age[1] lt yrange[1] then $ plots, [5,xrange[1]], [hor_age[1],hor_age[1]], COLOR=fsc_color('Green') if hor_age[0] gt yrange[0] and hor_age[0] lt yrange[1] then $ plots, [0, 5], [hor_age[0], hor_age[0]], COLOR=fsc_color('Green') endif mkeep = where(finite(met[keep])) yrange = minmax((met[keep])[mkeep]) if n_elements(hor_met) ne 0 then yrange = minmax([(met[keep])[mkeep],hor_met]) delta = yrange[1]-yrange[0] yrange[0] -= 0.1*delta yrange[1] += 0.1*delta if n_elements(metrange) gt 0 then yrange = metrange if not keyword_set(noytit) then ytit='[Fe/H], dex' else undefine, ytit if not keyword_set(noylab) then ytickformat='(f4.1)' plot, rad, met, YMARGIN=[3, -2.], XR=xrange, XST=3, $ /NODATA, YST=3, YR=yrange, $ YTICKNAME=ytickname, XTIT='R, arcsec', $ BACK=fsc_color('white'), COLOR=fsc_color('black') if n_elements(ytit) gt 0 then $ xyouts, ytpos, mean(!y.window), ytit, CHARSIZE=!p.charsize*0.5, /NORMAL, ALIGN=0.5, ORIENT=90 if keyword_set(ver) then begin plots, [ver, ver], [yrange], LINEST=2, COLOR=fsc_color('Royal Blue') endif for k=nprof-1,0,-1 do begin if cnt_posik[k] gt 0 then begin plotsym, 0, psymsize, THICK=2, /FILL ; circle of size=0.3 if codk[k] lt 3 then begin oploterror, rad[*posik[k]], met[*posik[k]], e_met[*posik[k]], PSYM=8, ERRCOLOR=redk[k], COLOR=redk[k] endif else begin rk = (*posik[k])[sort(rad[*posik[k]])] ;; to make the polyfill we need to sort the points and to cut the ;; things outside the x-y ranges, otherwise it plots outside the graph sortrad = sort(rad[rk]) rad[rk] = abs(rad[rk[sortrad]]) & met[rk] = met[rk[sortrad]] & e_met[rk] = 0.07 metmin = met[rk] - e_met[rk] metlim = where(metmin lt yrange[0], cnt) if cnt gt 0 then metmin[metlim] = yrange[0] metmax = met[rk] + e_met[rk] metlim = where(metmax gt yrange[1], cnt) if cnt gt 0 then metmax[metlim] = yrange[1] ; polyfill, /LINE_FIL, [rad[rk], reverse(rad[rk])], [metmin, reverse(metmax)], COLOR=fsc_color('dark gray') oplot, rad[rk], met[rk], PSYM=0, THICK=2.0, LINES=0 endelse endif endfor for k=nprof-1,0,-1 do begin if cnt_negak[k] gt 0 then begin plotsym, 0, psymsize, THICK=2, /FILL if codk[k] lt 3 then begin oploterror, -1e*rad[*negak[k]], met[*negak[k]], e_met[*negak[k]], PSYM=8, ERRCOLOR=bluk[k], COLOR=bluk[k] endif else begin rk = (*negak[k])[sort(-rad[*negak[k]])] ;; to make the polyfill we need to sort the points and to cut the ;; things outside the x-y ranges, otherwise it plots outside the graph sortrad = sort(rad[rk]) rad[rk] = abs(rad[rk[sortrad]]) & met[rk] = met[rk[sortrad]] & e_met[rk] = 0.07 metmin = met[rk] - e_met[rk] metlim = where(metmin lt yrange[0], cnt) if cnt gt 0 then metmin[metlim] = yrange[0] metmax = met[rk] + e_met[rk] metlim = where(metmax gt yrange[1], cnt) if cnt gt 0 then metmax[metlim] = yrange[1] ; polyfill,/line_fil, [rad[rk], reverse(rad[rk])], [metmin, reverse(metmax)], COLOR=fsc_color('dark gray') oplot, rad[rk], met[rk], PSYM=0, THICK=2.0, LINEST=2 ; oplot, -rad[rk], met[rk], PSYM=0, THICK=2.0, LINEST=2 endelse endif endfor if n_elements(hor_met) ne 0 then begin if n_elements(hor_met) ge 2 and $ hor_met[1] gt yrange[0] and hor_met[1] lt yrange[1] then $ plots, [5,xrange[1]], [hor_met[1],hor_met[1]], COLOR=fsc_color('Green') if hor_met[0] gt yrange[0] and hor_met[0] lt yrange[1] then $ plots, [0, 5], [hor_met[0], hor_met[0]], COLOR=fsc_color('Green') endif endif if keyword_set(ps) then begin niceplot, /RESET Device, /CLOSE set_plot, "x" endif !P.MULTI = 0 !P.CHARSIZE = 1 end ;----------------- end -------------------