;
; 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