#Function to deseasonalise time series through the normal quantile transform.
#The software does not allow to apply the transform back
#Arguments of the function:
#   - data: the data you wish to deseasonalise;
#   - wmin and wmax: minimu and maximum width of the deseasonalising window;
#   - step: number of observations in each year: this must be an integer!
#Output of the function: $datadest: deseasonalised time series
#                        $kolmog: kolmogorov test of normality applied to each period of the year
#                        $kolmogmax: maximum value of the kolmogorov test in order
#                         not to reject the hypothesis of normality
#			 $datasub: uncomment line 83 below in order to get in output the
#                         list containing the data in each deseasonalising window.


dest1=function(data,wmin,wmax,step)
{
kolmog=0
kolmogmax=0
ndata=length(data)
nanni=floor(length(data)/step)
if(((length(data)/step)-nanni)!=0) cat("Beware! The data do not cover an integer number of years!",fill=T)
pr1=seq(1,step)
pr2=rep(pr1,nanni)
datai=list()
rng=0
for (i in 1:step)
{
datai[[i]]=data[pr2==i]
for (ii in 1:(wmin/2))
{
iii=i+ii
if(iii>step) iii=iii-step
datai[[i]]=c(datai[[i]],data[pr2==iii])
}
for (ii in 1:(wmin/2))
{
iii=i-ii
if(iii<1) iii=iii+step
datai[[i]]=c(data[pr2==iii],datai[[i]])
}
rng[i]=range(datai[[i]])[2]-range(datai[[i]])[1]
}
minrng=min(rng)
maxrng=max(rng)
for (i in 1:step)
{
wagg=0
wper=wmin+(wmax-wmin)/(maxrng-minrng)*(rng[i]-minrng)
if(wper>wmin)
{
wagg=(wper-wmin)/2
for(ii in 1:(ceiling(wagg)))
{
iii=i+wmin/2+ii
if(iii>step) iii=iii-step
datai[[i]]=c(datai[[i]],data[pr2==iii])
}
for(ii in 1:(ceiling(wagg)))
{
iii=i-wmin/2-ii
if(iii<1) iii=iii+step
datai[[i]]=c(data[pr2==iii],datai[[i]])
}
pr10=ceiling(wagg)-wagg
pr11=floor(pr10*nanni)
datai[[i]]=datai[[i]][pr11:(length(datai[[i]])-pr11)]
}
}
datat=data
for (i in 1:step)
{
casrtux=approx(sort(datai[[i]]),qnorm(ppoints(datai[[i]])),xout=data[pr2==i])$y
pr500=ppoints(casrtux)
pr600=pnorm(sort(casrtux),0,1)
kolmog[i]=max(abs(pr600-pr500))
kolmogmax[i]=1.36/sqrt(length(casrtux))

datat[pr2==i]=casrtux

}
out=list()
out$datadest=datat
#out$datasub=datai
out$kolmog=kolmog
out$kolmogmax=kolmogmax
return(out)
}