;
; Plots low-frequency timeseries and their differences
;
;plot,[0,1]
multi_plot,nrow=5
loadct,39
def_1color,r,g,b,1,color='red'
def_1color,r,g,b,2,color='blue'
if !d.name eq 'X' then begin
  window,ysize=900
endif else begin
  device,xoffset=2,xsize=17
endelse
;
thalf=10.
refperiod=[1881,1975]
;
; Get region info
;
restore,filename='regboxes.idlsave'
for ireg = 0 , nreg-1 do begin
  restore,filename='instradj_'+regname(ireg)+'.idlsave'
  wdens=instradj
  mknormal,wdens,x,refperiod=refperiod,refmean=refmean,refsd=refsd
  filter_cru,thalf,tsin=wdens,tslow=instrlow,tshigh=instrhi,/nan
  ;
  restore,filename='densadj_'+regname(ireg)+'.idlsave'
  mknormal,densadj,x,refperiod=refperiod,refmean=refmean,refsd=refsd
  filter_cru,thalf,tsin=densadj,tslow=denslow,tshigh=denshi,/nan
  ;
  restore,filename='rwidadj_'+regname(ireg)+'.idlsave'
  mknormal,rwidadj,x,refperiod=refperiod,refmean=refmean,refsd=refsd
  filter_cru,thalf,tsin=rwidadj,tslow=rwidlow,tshigh=rwidhi,/nan
  ;
  pause
  plot,x,instrlow,$
    xrange=[1881,1992],xstyle=1,xtitle='Year',$
    yrange=[-3,3],ystyle=1,ytitle='Low-frequency normalised anomalies',$
    title='Region: '+regname(ireg)+'  Filter: >'+$
      string(thalf,format='(I2)')+' years'
  oplot,x,denslow,thick=5
  oplot,x,rwidlow,thick=3,linestyle=1
  legend,['Apr-Sep Temperature','Density','Width'],thick=[1,5,3],$
    linestyle=[0,0,1]
  ;
  plot,x,denslow-instrlow,thick=5,$
    xrange=[1881,1992],xstyle=1,xtitle='Year',$
    yrange=[-3,3],ystyle=1,ytitle='Low-frequency normalised difference',$
    title='Region: '+regname(ireg)+'  Filter: >'+$
      string(thalf,format='(I2)')+' years'
  oplot,!x.crange,[0.,0.]
  oplot,x,rwidlow-instrlow,thick=3,linestyle=1
  legend,['Density - Temperature','Width - Temperature'],thick=[5,3],$
    linestyle=[0,1]
  ;
  ; Now fit trend lines to 40 and 60 year segments
  ;
  nlin=1992-1881+1-40
  a=fltarr(nlin)
  for ii = 0 , nlin-1 do begin
    xx=denslow-instrlow
    acoeff=linfit(findgen(40),xx(1881-1600+ii:1881-1600+ii+39))
    a(ii)=acoeff(1)
  endfor
  plot,findgen(nlin)+1881+20,a,xrange=[1881,1992],/xstyle,$
    /ystyle,yrange=[-0.06,0.06]
  oplot,!x.crange,[0.,0.]
  ;
  nlin=1992-1881+1-60
  a=fltarr(nlin)
  for ii = 0 , nlin-1 do begin
    xx=denslow-instrlow
    acoeff=linfit(findgen(60),xx(1881-1600+ii:1881-1600+ii+59))
    a(ii)=acoeff(1)
  endfor
  oplot,findgen(nlin)+1881+30,a,linestyle=2
  legend,['Slope of 40-yr segment','Slope of 60-yr segment'],linestyle=[0,2]
  ;
  ; now plot correlations (low + hi freq) of 50-yr segments
  ;
  nlin=1992-1881+1-50
  al=fltarr(nlin)
  ah=fltarr(nlin)
  for ii = 0 , nlin-1 do begin
    xx=indgen(50)+1881-1600+ii
    al(ii)=correlate(denslow(xx),instrlow(xx))
    ah(ii)=correlate(denshi(xx),instrhi(xx))
  endfor
  plot,findgen(nlin)+1881+25,al,xrange=[1881,1992],/xstyle,$
    /ystyle,yrange=[-1.,1.]
  oplot,findgen(nlin)+1881+25,ah,linestyle=2
  oplot,!x.crange,[0.,0.]
  legend,['r(low) of 50-yr segment','r(high) of 50-yr segment'],linestyle=[0,2]
  ;
  ; now plot rmse of two straight-line segments (joining at a point)
  ;
  minlen=20
  nlin=1992-1881+1-2.*minlen      ; need at least a 20-yr period for each line
  armse=fltarr(nlin)
  xx=denslow-instrlow
  for ii = 0 , nlin-1 do begin
    print,ii
    ddd=findbest(x,xx,break=ii+minlen-1)
    armse(ii)=ddd
  endfor
  plot,findgen(nlin)+1881+minlen,al,xrange=[1881,1992],/xstyle,$
    ytitle='rms error of 2-line fit'
;    /ystyle,yrange=[-1.,1.]
  ;
endfor
;
end