;
; Computes variance versus n for densities
; FOR ADJUSTED TIMESERIES
;
loadct,39
sc=findgen(8)*30.+40.  &  sc=sc*(!d.n_colors-1)/255.
ss=[1,4,6,1,4,6,1,4]
yr=[0.25,0.25,0.4,0.4,0.8,0.5,0.5,1.5]
multi_plot,nrow=4
if !d.name eq 'X' then begin
  initx
  window,ysize=850
endif
;
; Get details of regions & rbar
;
restore,filename='reglists.idlsave'
restore,filename='rbar_dens_full.idlsave'
;
; Do each region separately
;
for i = 0 , nreg-1 do begin
  ;
  ; Read in timeseries of densities
  ;
  restore,filename='densadj_'+regname(i)+'.idlsave'
  keeplist=where(x lt 1991,nkeep)
  y=densadj(keeplist)
  x=x(keeplist)
  ;
  ; Compute mean n (chrons), ncore (cores) and variance from overlapping 50 yrs
  ;
  nyr=n_elements(x)
  segsize=51
  nout=nyr-segsize+1
  runx=fltarr(nout)
  runvar=fltarr(nout)
  runn=fltarr(nout)
  runncore=fltarr(nout)
  for j = 0 , nyr-segsize do begin
    yseg=y(j:j+segsize-1)
    xseg=x(j:j+segsize-1)
    nseg=n(j:j+segsize-1)
    ncoreseg=neff(j:j+segsize-1)
    runx(j)=total(xseg)/total(finite(xseg))
    runn(j)=total(nseg)/total(finite(yseg))  ; only mean n where we have data
    runncore(j)=total(ncoreseg)/total(finite(yseg))
    keeplist=where(finite(yseg),nkeep)
    if nkeep lt segsize then begin
      runvar(j)=!values.f_nan
    endif else begin
      ygot=yseg(keeplist)
      dummy=moment(ygot)
      runvar(j)=dummy(1)
    endelse
  endfor
  ;
  ; Now plot the results
  ;
  maxn=max(runn,/nan,k)
  maxvar=runvar(k)
  pause
  plot,x,y,xrange=[1600,2000],yrange=[-3,3],$
    title=regname(i),ytitle='max density timeseries'
  plot,runx,runvar,xrange=[1600,2000],yrange=[0.,yr(i)],$
    ytitle='running variance'
  plot,runx,runn,xrange=[1600,2000],ytitle='no. of chronologies'
  plot,runn,runvar,xrange=[0,100],psym=1,yrange=[0,3],$
    ytitle='running variance',xtitle='no. of chronologies'
; if i eq 0 then plot,(1.+(runn-1.)*allrbar(i))/runn,runvar,$
;   xrange=[0,0.6],psym=1,yrange=[0,6], $
;   xtitle='1 / effective mean no. of chronologies',/xstyle, $
;   ytitle='running variance' $
; else oplot,(1.+(runn-1.)*allrbar(i))/runn,runvar,psym=ss(i),color=sc(i)
  ; Keep the data for n<= 50 for fitting a line
  if i eq 0 then begin
    ; For all region, n > 50 always. Just initialise the arrays
    fitx=[0.]
    fity=[0.]
  endif else begin
    keeplist=where(finite(runvar) and (runn le 50),nkeep)
    if nkeep gt 0 then begin
      fitx=[fitx,runn(keeplist)]
      fity=[fity,1./runvar(keeplist)]
    endif
  endelse
endfor
;plot,[0,1],/nodata,xstyle=5,ystyle=5
;legend,regname,textcolors=sc,colors=sc,psym=ss
;
end