; See "http://paos.colorado.edu/research/wavelets/"
; Written January 1998 by C. Torrence
;
; Computes wavelet power spectra for long chronologies from Taimyr
; and Tornetraske and Yamal
;
iuse=2
mtit=['Tornetrask','Yamal','Taimyr']
fn=['tornad.rcs','yamal.rcs','taimyr2.rcs']
nts=n_elements(fn)
allnyr=2000
allyear=findgen(allnyr)
allts=fltarr(nts,allnyr)
;
; Repeat for each series
;
for i = 0 , nts-1 do begin
  ;
  ; First read in the appropriate data
  ;
  yrsten=intarr(2)
  openr,1,fn(i)
  readf,1,yrsten,format='(6X,2I4)'
  y1=yrsten(0) & adj1=(y1 mod 10) & y1=y1-adj1
  y2=yrsten(1) & adj2=9-(y2 mod 10) & y2=y2+adj2
  print,fn(i),yrsten,y1,y2
  nyr=y2-y1+1
  rawdat=fltarr(2,nyr)
  onedat=fltarr(2,10)
  for j = 0 , nyr-1 , 10 do begin
    readf,1,onedat,format='(10X,10(I4,I3))'
    rawdat(*,j:j+9)=onedat(*,*)
  endfor
  close,1
  ;
  ; Now sort out the missing data
  ;
  timey=findgen(nyr)+y1
  ts=reform(rawdat(0,*))
  ml=where(ts eq 9990,nmiss)
  if nmiss gt 0 then ts(ml)=!values.f_nan
  ;
  ; Now normalise relative to their entire length
  ;
  mknormal,ts
  allts(i,*)=ts(*)
  ;
endfor
;
sst=reform(allts(iuse,*))
kl=where(finite(sst),n)
sst=sst(kl)
time=allyear(kl)
;
;------------------------------------------------------ Computation
;
; normalize by standard deviation (not necessary, but makes it easier
; to compare with plot on Interactive Wavelet page, at
; "http://paos.colorado.edu/research/wavelets/plot/"
	sst = (sst - TOTAL(sst)/n)

	dt = 1.         ; annual data
	xrange = [min(time),max(time)]  ; plotting range
	pad = 1
	s0 = 2.*dt    ; this says start at a scale of 1 year
	dj = 0.25  ; this will do 4 sub-octaves per octave
	j1 = 9./dj  ; this says do 9 powers-of-two with dj sub-octaves each
	mother = 'Morlet'
        mothper=6
	recon_sst = sst   ; save an extra copy, so we don't erase original sst

; estimate lag-1 autocorrelation, for red-noise significance tests
; Note that we actually use the global wavelet spectrum (GWS)
; for the significance tests, but if you wanted to use red noise,
; here's how you could calculate it...
	lag1 = (A_CORRELATE(sst,1) + SQRT(A_CORRELATE(sst,2)))/2.
	
; Wavelet transform:
	wave = WAVELET(recon_sst,dt,PERIOD=period,SCALE=scale,S0=s0, $
                       param=mothper,$
		PAD=pad,COI=coi,DJ=dj,J=j1,MOTHER=mother,/RECON)
	power = (ABS(wave))^2  ; compute wavelet power spectrum
	global_ws = TOTAL(power,1)/n   ; global wavelet spectrum (GWS)
	J = N_ELEMENTS(scale) - 1

; Significance levels, assuming the GWS as background spectrum:
	signif = WAVE_SIGNIF(sst,dt,scale,0, $
		GWS=global_ws,SIGLVL=0.90,MOTHER=mother)
	signif = REBIN(TRANSPOSE(signif),n,J+1)  ; expand signif --> (J+1)x(N) array
	signif = power/signif   ; where ratio > 1, power is significant

; GWS significance levels:
	dof = n - scale   ; the -scale corrects for padding at edges
	global_signif = WAVE_SIGNIF(sst,dt,scale,1, $
		LAG1=lag1,DOF=dof,MOTHER=mother,CDELTA=Cdelta,PSI0=psi0)

; check total variance (Parseval's theorem) [Eqn(14)]
	scale_avg = REBIN(TRANSPOSE(scale),n,J+1)  ; expand scale-->(J+1)x(N) array
	power_norm = power/scale_avg
	variance = (MOMENT(sst))(1)
	recon_variance = dj*dt/(Cdelta*n)*TOTAL(power_norm)  ; [Eqn(14)]
	
	IF (N_ELEMENTS(recon_sst) GT 1) THEN BEGIN
		recon_variance = (MOMENT(recon_sst))(1)
; RMS of Reconstruction [Eqn(11)]
		rms_error = SQRT(TOTAL((sst - recon_sst)^2)/n)
		PRINT
		PRINT,'        ******** RECONSTRUCTION ********'
		PRINT,'original variance =',variance,' degC^2'
		PRINT,'reconstructed var =',FLOAT(recon_variance),' degC^2'
		PRINT,'Ratio = ',recon_variance/variance
		PRINT,'root-mean-square error of reconstructed sst = ',rms_error,' degC'
		PRINT
		IF (mother EQ 'DOG') THEN BEGIN
			PRINT,'Note: for better reconstruction with the DOG, you need'
			PRINT,'      to use a very small s0.'
		ENDIF
		PRINT
	ENDIF
	
; Scale-average between El Nino periods of 2--8 years
	avg = WHERE((scale GE 2) AND (scale LT 8))
	scale_avg = dj*dt/Cdelta*TOTAL(power_norm(*,avg),2)  ; [Eqn(24)]
	scaleavg_signif = WAVE_SIGNIF(sst,dt,scale,2, $
		GWS=global_ws,SIGLVL=0.90,DOF=[2,7.9],MOTHER=mother)


;------------------------------------------------------ Plotting

if !d.name eq 'X' then begin
  !P.FONT = -1
  !P.CHARSIZE = 1
  WINDOW,0,XSIZE=600,YSIZE=600
endif else begin
	DEVICE,/PORT,/INCH,XSIZE=6.5,XOFF=1,YSIZE=6,YOFF=3,/COLOR,BITS=8
	!P.FONT = 0
	!P.CHARSIZE = 0.75
endelse
!P.MULTI = 0
!X.STYLE = 1
!Y.STYLE = 1
LOADCT,39
cfac=(!d.n_colors-1)/255.

;--- Plot time series
	pos1 = [0.1,0.75,0.7,0.95]
	PLOT,time,sst,XRANGE=xrange, $
		XTITLE='Time (year)',YTITLE='(dimensionless)', $
		TITLE=mtit(iuse)+' TRW chronology', $
		POSITION=pos1
  IF (N_ELEMENTS(recon_sst) GT 1) THEN OPLOT,time,recon_sst,COLOR=144*cfac
;	XYOUTS,0.85,0.9,/NORMAL,ALIGN=0.5, $
;		'!5WAVELET ANALYSIS!X'+$
;		'!C!CC. Torrence & G.P. Compo'+$
;		'!C!Chttp://paos.colorado.edu/!Cresearch/wavelets/'

;--- Contour plot wavelet power spectrum
	yrange = [1024,2.]   ; years
	levels = [0.25,1,4,16,64]
	colors = [64,128,185,208,254]*cfac
	period2 = FIX(ALOG(period)/ALOG(2))   ; integer powers of 2 in period
	ytickv = 2.^(period2(UNIQ(period2)))  ; unique powers of 2
        kl=where((ytickv le yrange(0)) and (ytickv ge yrange(1)))
        ytickv=ytickv(kl)
	pos2 = [pos1(0),0.35,pos1(2),0.65]

	CONTOUR,power,time,period,/NOERASE,POSITION=pos2, $
		XRANGE=xrange,YRANGE=yrange,/ystyle,/ylog, $
		YTICKS=N_ELEMENTS(ytickv)-1,YTICKV=ytickv, $
		LEVELS=levels,C_COLORS=colors,/FILL, $
		XTITLE='Time (year)',YTITLE='Period (years)', $
		TITLE='b) Wavelet Power Spectrum (contours at 0.25,1,4,16,64]'
print,!y.crange
; significance contour, levels at -99 (fake) and 1 (significant)
	CONTOUR,signif,time,period,/OVERPLOT,LEVEL=1,THICK=2, $
		C_LABEL=1,C_ANNOT='90%',C_CHARSIZE=1
; cone-of-influence, anything "below" is dubious
	x = [time(0),time,MAX(time)]
	y = [MAX(period),coi,MAX(period)]
	color = !p.background
 	POLYFILL,x,y,COLOR=color,NOCLIP=0
 	POLYFILL,x,y,COLOR=color,NOCLIP=0
	PLOTS,time,coi,NOCLIP=0,THICK=1

;--- Plot global wavelet spectrum
	pos3 = [0.74,pos2(1),0.95,pos2(3)]
	blank = REPLICATE(' ',29)
	PLOT,global_ws,period,/NOERASE,POSITION=pos3, $
		THICK=2,XSTYLE=10,YSTYLE=9, $
		YRANGE=yrange,/YTYPE,YTICKLEN=-0.02, $
		XTICKS=2,XMINOR=2, $
		YTICKS=N_ELEMENTS(ytickv)-1,YTICKV=ytickv,YTICKNAME=blank, $
		XTITLE='Power',TITLE='c) Global'
	OPLOT,global_signif,period,LINES=1
	XYOUTS,2,8,'95%'
	
;--- Plot 2--8 yr scale-average time series
	pos4 = [pos1(0),0.05,pos1(2),0.25]
	PLOT,time,scale_avg,/NOERASE,POSITION=pos4, $
		XRANGE=xrange,YRANGE=[0,MAX(scale_avg)*1.25],THICK=2, $
		XTITLE='Time (year)',YTITLE='Avg variance', $
		TITLE='d) 2-8 yr Scale-average Time Series'
	OPLOT,xrange,scaleavg_signif+[0,0],LINES=1

END