;
; Reads in site-by-site MXD and temperature series in 5 yr blocks, all
; correctly normalised etc.  Rotated PCA is performed to obtain the 'decline'
; signal!
;
;----------------------------------------------------
;
restore,filename='briffa_sep98_decline1.idlsave'
;  allmxd,alltemp,mxd5yr,temp5yr,nchron,nyr,nblock,statlat,statlon
;  blockyr,timey
;
nret=4
diff5yr=mxd5yr-temp5yr
;
myeof1d,diff5yr,ev,ea,lam,lampct,lamcum,nretain=nret,correlation=1
myeof1d_rot,diff5yr,ev,ea,lam,lampct,lamcum,$
  evrot,earot,lamrot,lampctrot,lamcumrot,$
  nretain=nret,correlation=1
ev=evrot
ea=earot
lam=lamrot
lampct=lampctrot
lamcum=lamcumrot
;
; Scale EOFs and PCs so that highest value of each EOF is +1
;
for i = 0 , nret-1 do begin
  dummy=max(abs(ev(*,i)),imax)
  maxval=ev(imax,i)
  ev(*,i)=ev(*,i)/maxval
  ea(*,i)=ea(*,i)*maxval
endfor
;
; Sort modes into order of descending variance explained
;
i=reverse(sort(lampct))
lam=lam(i)
lampct=lampct(i)
lamcum(0)=lampct(0)
for j = 1 , nret-1 do lamcum(j)=lamcum(j-1)+lampct(j)
ev=ev(*,i)
ea=ea(*,i)
;
; Now prepare for plotting
;
loadct,39
multi_plot,nrow=2,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=12
endelse
def_1color,20,color='red'
def_1color,21,color='deepblue'
def_1color,22,color='lgrey'
def_1color,95,color='dgreen'
cpl_usersym,/circle,/fill
;
map=def_map(/npolar)  &  map.limit(0)=25.
coast=def_coast(/get_device) & coast.double=1
labels=def_labels(/off)
;
; Plot variance stats first
;
if 2 eq 3 then begin
plot,findgen(nret)+1.,lam,xtitle='Mode',$
  ytitle='Eigenvalue',$
  title='Correlation matrix PCA of (MXD minus TEMP)'
;
plot,findgen(nret)+1.,alog(lam),xtitle='Mode',$
  ytitle='Log of eigenvalue'
;
plot,findgen(nret+1),[!values.f_nan,lampct],xtitle='Mode',$
  ytitle='Explained variance  (%)',psym=10
;
plot,findgen(nret+1),[0.,lamcum],xtitle='Mode',$
  ytitle='Cumulative explained variance  (%)',psym=10,yrange=[0,100]
endif
;
; Now plot some of the retained EOFs and PCs
;
nret=4
;
for i = 1 , 2 do begin
  ;
  pause
  inter_boxfd,/nodata,map=map,coast=coast,labels=labels,$
  title='EOF mode #'+string(i+1,format='(I2)')+'  ('+$
    string(round(lampct(i)),format='(I3)')+' % )'
  ;
  for j = 0 , nchron-1 do begin
    onev=ev(j,i)
    if onev lt 0 then scol=21 else scol=20
    ssiz=abs(onev)*0.8       ; scaling factor to get good sizes
    ssiz=ssiz+0.1            ; minimum size!
    map_plots,statlon(j),statlat(j),psym=8,color=scol,symsize=ssiz
  endfor
  ;
  plot,blockyr,ea(*,i),/xstyle,thick=2,$
    xtitle='Year',ymargin=[10,10],xmargin=[4,4]
  oplot,blockyr,smooth(ea(*,i),5,/nan,/edge_truncate),thick=8,color=95
  ;
endfor
;
end