; ; Reads in the output from the Matlab cross-spectral ; analysis program and plots it ; ; Prepare for plotting ; loadct,39 multi_plot,nrow=3,ncol=2,layout='large' if !d.name eq 'X' then begin window,ysize=850 !p.font=-1 endif else begin !p.font=0 device,/helvetica,/bold,font_size=9 endelse def_1color,20,color='lgrey' def_1color,21,color='white' ; spawn,'wc specout.dat',osoutput nfre=fix(osoutput)/4. nfre=nfre(0) print,nfre x=fltarr(nfre) alldat=fltarr(8,nfre) tmag=fltarr(nfre) tang=fltarr(nfre) ; openr,1,'specout.dat' readf,1,x readf,1,alldat readf,1,tmag readf,1,tang close,1 ; plot,x,alldat(0,*),$ /xstyle,xtitle='Frequency (/yr)',$ ytitle='Hugershoff Power spectral density',ylog=1,yrange=[0.01,13.],$ /ystyle r95=reform(alldat(5,*)) oplot,x,alldat(0,*)-r95 > 0.0001,linestyle=1 oplot,x,alldat(0,*)+r95,linestyle=1 ; plot,x,alldat(1,*),$ /xstyle,xtitle='Frequency (/yr)',$ ytitle='Age-banded Power spectral density',ylog=1,thick=3,yrange=[0.01,13.],$ /ystyle r95=reform(alldat(6,*)) oplot,x,alldat(1,*)-r95 > 0.,linestyle=1 oplot,x,alldat(1,*)+r95,linestyle=1 ; plot,x,alldat(0,*),$ /xstyle,xtitle='Frequency (/yr)',$ ytitle='Power spectral density',ylog=1,yrange=[0.01,13.],/ystyle oplot,x,alldat(1,*),thick=3 oplot,x,alldat(2,*),thick=5 legpos=convert_coord(0.25,0.8,/data,/to_normal) legend,['Hug','Age','Cross'],thick=[1,3,5],$ position=legpos ; plot,x,tmag,$ /xstyle,xtitle='Frequency (/yr)',$ ytitle='Transfer function magnitude' oplot,!x.crange,[1.,1.] ; plot,x,tang,$ /xstyle,xtitle='Frequency (/yr)',$ ytitle='Transfer function phase',/ystyle,yrange=[-90.,90.] oplot,!x.crange,[0.,0.] ; plot,x,alldat(4,*),$ /xstyle,xtitle='Frequency (/yr)',$ ytitle='Coherence' ; ; Now make a single plot with a split linear scales ; multi_plot,nrow=4,ncol=1 if !d.name eq 'PS' then begin !p.font=0 device,/helvetica,/bold,font_size=14 endif !y.margin=[0,0] !y.omargin=[0,4] ; pause ytop=15.5 ycut=0.5 ybot=0. y1=reform(alldat(0,*)) yr=reform(alldat(5,*)) y2=reform(alldat(1,*)) filx=[x,reverse(x)] fily=[y1-yr,reverse(y1+yr)] plot,x,y1,/nodata,$ /xstyle,xtickformat='nolabels',$ /ystyle,yrange=[ycut,ytop] polyfill,filx,fily,noclip=0,color=20 oplot,x,y1,thick=2 oplot,x,y2,thick=5 ; legend,['Hugershoff MXD spectrum and 95% range','Age Band Decomposition MXD spectrum'],$ thick=[2,5],position=[0.2,0.95] ; plot,x,y1,/nodata,$ /xstyle,xtitle='Frequency (/yr)',$ /ystyle,yrange=[ybot,ycut] polyfill,filx,fily > ybot,noclip=0,color=20 oplot,x,y1,thick=2 oplot,x,y2,thick=5 xyouts,align=0.5,charsize=0.6,0.148,0.22,'6.7' xyouts,align=0.5,charsize=0.6,0.187,0.125,'5.3' xyouts,align=0.5,charsize=0.6,0.27,0.145,'3.7' xyouts,align=0.5,charsize=0.6,0.41,0.10,'2.4' ; xyouts,-0.04,ycut,orient=90.,'Spectral density [ (!Uo!NC)!U2!N/yr ]',$ align=0.5 ; ; Now read the higher-resolution spectrum ; spawn,'wc specout256.dat',osoutput nfre=fix(osoutput)/4. nfre=nfre(0) print,nfre x=fltarr(nfre) alldat=fltarr(8,nfre) tmag=fltarr(nfre) tang=fltarr(nfre) ; openr,1,'specout256.dat' readf,1,x readf,1,alldat readf,1,tmag readf,1,tang close,1 ; !p.region=[0.3,0.75,0.92,0.88] polyfill,!p.region([0,2,2,0,0]),!p.region([1,1,3,3,1]),color=21,/normal ytop=30. ycut=0.5 ybot=0. y1=reform(alldat(0,*)) yr=reform(alldat(5,*)) y2=reform(alldat(1,*)) filx=[x,reverse(x)] fily=[y1-yr,reverse(y1+yr)] plot,x,y1,/nodata,$ /xstyle,xrange=[0,0.1],xtickformat='nolabels',$ /ystyle,yrange=[ycut,ytop] polyfill,filx,fily,noclip=0,color=20 oplot,x,y1,thick=2,noclip=0 oplot,x,y2,thick=5,noclip=0 xyouts,align=0.5,charsize=0.6,0.035,3.,'24-37' !p.region=[0.3,0.62,0.92,0.75] polyfill,!p.region([0,2,2,0,0]),!p.region([1,1,3,3,1]),color=21,/normal plot,x,y1,/nodata,$ /xstyle,xrange=[0,0.1],xtitle='Frequency (/yr)',$ /ystyle,yrange=[ybot,ycut] polyfill,filx,fily,noclip=0,color=20 oplot,x,y1,thick=2,noclip=0 oplot,x,y2,thick=5,noclip=0 axis,xaxis=0,/xstyle,xtickformat='nolabels' xyouts,-0.01,ycut,orient=90.,'Spectral density [ (!Uo!NC)!U2!N/yr ]',$ align=0.5,charsize=0.6 ; xyouts,align=0.5,charsize=0.6,0.0155,0.36,'60-80' xyouts,align=0.5,charsize=0.6,0.058,0.4,'17' xyouts,align=0.5,charsize=0.6,0.07,0.36,'14' ; !p.region=[0.,0.,0.,0.] !y.margin=[4,2] !y.omargin=[0,0] ; end