#Scafetta Figure 2. Total solar irradiance sequences and their
#maximal overlap MRA decomposition according to equation
#(2). The data for the period 1980–2000 are given
#by TSI satellite composites: (a) ACRIM [Willson and
#Mordvinov, 2003], and (b) PMOD [Fro¨hlich and Lean,
#1998]. The period 1900–1980 is covered by the TSI proxy
#reconstruction by Lean et al. [1995] that has been adjusted
#by means of a vertical shift to match the two TSI satellite
#composites in the year 1980.

#Schmidt

#‘‘mra’’ function from the waveslim
#package for R to compute the wavelet components (‘‘D =
#mra(y, wf = ‘‘la8’’, J = floor(log(length(y),2)), method =
#‘‘modwt’’)’’, taking D4(t) = D$D8 and D3(t) = D$D7.),
#which should be si

#Get Solar Data (making Acrim and pmod into monthly series)
source("http://www.climateaudit.org/scripts/solar/collation.functions.solar.txt")
	pmod=get.solar(dset="pmod")
	acrim=get.solar(dset="acrim")
	lean=get.solar(dset="lean95");tsp(lean) #1600 1999

#Make annual averages of monthly information 
  	acrim.annual= annavg(acrim)
	pmod.annual= annavg(pmod)

#Do splice a la Scafetta-West 2006
	Solar=ts.union(window(lean,start=1881),acrim.annual,pmod.annual)	
	Solar=data.frame(year=time(Solar),Solar)
	names(Solar)=c("year","lean","acrim","pmod")
	(delta.acrim= Solar$acrim[1980-1880]-Solar$lean[1980-1880]) #  -1.476112
	(delta.pmod= Solar$pmod[1980-1880]-Solar$lean[1980-1880]) #  -1.476112
	   #SW centers on 1980
	temp= Solar$year>=1980
	Solar$acrim.splice=Solar$lean+delta.acrim
	Solar$acrim.splice[temp]=Solar$acrim[temp]
	Solar$pmod.splice=Solar$lean+delta.pmod
	Solar$pmod.splice[temp]=Solar$pmod[temp]
	
#############
##wavelet.decomposition with reflection padding
###################

library(waveslim)
 
#this pads the series by reflection before wavelet decomposition and truncates back
  wavelet.decomposition<-function(x,wf0="la8",pad="reflect") {
	N<-length(x)
      #steps to interpolate missing data not relevant here
	#f<-approxfun(1:N,x)
	#y<-f(1:N)
	y=x
	temp<-!is.na(y)
	ybar<-mean(y,na.rm=TRUE)
	y<-y[temp]-ybar
	N1<-length(y)
	J0<-trunc(log(N1,2))
	pad=2^(J0-1)
	z= c( rev( y[1: pad]),y,rev( y[(N-pad+1):N]) );length(z) 
   	mod.y<-mra(z,wf0,J0,"modwt") 
	names(mod.y) #[1] "D1" "D2" "D3" "D4" "D5" "D6" "D7" "D8" "S8"
	test<-mod.y[[1]]
	for (i in 2:(J0+1)) test<-cbind(test,mod.y[[i]])
	dimnames(test)[[2]]=names(mod.y)
	test=test[(pad+1): (N+pad) ,]
	return(test)
 }

#############
# NOW Wavelet on series long enough for 2^7
#############################

	length(1881:2008)  #[1]128
 	temp=(1881:2008)-Solar$year[1]+1
	model=wavelet.decomposition(Solar$acrim.splice[temp] ,wf0="la8")
	(ybar=mean(Solar$acrim.splice[temp] ) )#  1365.152
		dim(model)  #[1] 256   8
	decomp=cbind(R2=apply(model[,c("D1","D2")],1,sum), D3=model[,"D3"],D4=model[,"D4"],S4=apply(model[,5:ncol(model)],1,sum) )

library(GDD)
#GDD(file="d:/climate/images/2009/models/acrim_splice.gif",type="gif",h=420,w=420)
	par(mar=c(3,4,2,1))
	plot(1881:2008, Solar$acrim.splice[temp],type="l", ylab="Irradiance (wm-2)")
	title("ACRIM Irradiance") 
	lines(1881:2008,decomp[,"S4"]+ybar,col=2)
	mtext(side=1,"ACRIM Spliced with Lean 1995 per Scafetta-West 2006",cex=.7,line=1.7)
	abline(v=2000,lty=3)
#dev.off()