;
; Plots low-frequency timeseries and their differences
;
;plot,[0,1]
multi_plot,nrow=3
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.
pp=[1881,1940]
refp=[ [1881,1975] , [pp] , [pp] , [pp] , [pp] , [pp] , [pp] , [pp] , [pp] , $
  [pp] , [pp] ]
;
; Get region info
;
restore,filename='regboxes.idlsave'
for ireg = 0 , nreg-1 do begin
  refperiod=reform(refp(*,ireg))
  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  Ref:'+$
      string(refperiod,format='(2I5)')
  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]
  ;
  if regname(ireg) eq 'ALL' then begin
    openw,1,regname(ireg)+'_tdwi.dat'
    printf,1,regname(ireg),format='(A10,"   Region name")'
    kl=where(x ge 1881,nk)
    printf,1,nk,format='(I4,6X,"   Number of years")'
    xout=x(kl)
    y1out=denslow(kl)
    y2out=rwidlow(kl)
    y3out=instrlow(kl)
    missout=-99.9999
    printf,1,missout,format='(F9.4,X,"   Missing value")'
    for i = 0 , nk-1 do printf,1,xout(i),y1out(i),y2out(i),y3out(i),$
      format='(I5,3F9.4)'
    close,1
  endif
  ;
  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',noclip=1
  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
  ;
  seglen=50
  nlin=1992-1881+1-seglen
  al=fltarr(nlin)
  ah=fltarr(nlin)
  for ii = 0 , nlin-1 do begin
    xx=indgen(seglen)+1881-1600+ii
    al(ii)=correlate(denslow(xx),instrlow(xx))
    ah(ii)=correlate(denshi(xx),instrhi(xx))
  endfor
  plot,findgen(nlin)+1881+seglen/2.,al,xrange=[1881,1992],/xstyle,$
    /ystyle,yrange=[-0.1,1.]
  oplot,findgen(nlin)+1881+seglen/2.,ah,linestyle=2
  oplot,!x.crange,[0.,0.]
  titst=string(seglen,format='(I3)')+'-yr segment'
  legend,['r(low) of '+titst,'r(high) of '+titst,'r(15-yr) of '+titst],$
    linestyle=[0,2,0],thick=[1,1,5]
  ;
  ; and overplot correlation of 15-yr segments (implicitly high freq.)
  ;
  seglen=15
  nlin=1992-1881+1-seglen
  a15=fltarr(nlin)
  a15h=fltarr(nlin)
  for ii = 0 , nlin-1 do begin
    xx=indgen(seglen)+1881-1600+ii
    a15(ii)=correlate(densadj(xx),instradj(xx))
    a15h(ii)=correlate(denshi(xx),instrhi(xx))
  endfor
  oplot,findgen(nlin)+1881+seglen/2.,a15,thick=5
  oplot,findgen(nlin)+1881+seglen/2.,a15h,thick=5,linestyle=1
  ;
  ; 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