# This is for doing Rybski analysis of climate variables


readCRU<-function(source="../data/crutem3vgl.txt",temps=2:13,plot=T) {
f<-read.table(source,fill=TRUE)
d<-as.vector(t(as.matrix(f[seq(1,length(f[,1]),by=2),temps])))
d<-d[d!= 0]
#d<-c(d,-0,-0)
ts(d,frequency=length(temps),start=1850)
}


readGISS<-function(source="data/GISS.txt",temps=2:13,plot=T) {
f<-read.table(source,fill=TRUE,skip=10)
f<-f[f$V1 %in% 1880:2007,]
d<-as.numeric(t(as.matrix(f[seq(1,length(f[,1])),temps])))*0.01
ts(d,start=1881,frequency=1)
}

readRSS<-function(skip=3,file=UAH,col=3) {
	f<-read.table(file,fill=TRUE,skip=skip)
	ts(f[col],start=f$V1[1],frequency=12)
}

TLT<-"ftp://ftp.ssmi.com/msu/monthly_time_series/rss_monthly_msu_amsu_channel_tlt_anomalies_land_and_ocean_v03_1.txt"
TLS<-"ftp://ftp.ssmi.com/msu/monthly_time_series/rss_monthly_msu_amsu_channel_tls_anomalies_land_and_ocean_v03_0.txt"
TLTo<-"ftp://ftp.ssmi.com/msu/monthly_time_series/rss_monthly_msu_amsu_channel_tlt_anomalies_ocean_v03_1.txt"
UAH<-"data/tltglhmam_5.2.txt"
CRU<-"http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3vgl.txt"


#tlt<-readRSS(file=TLT)
#crut<-readCRU(file=,temps=14)


rho<-function(l=1,h=0.9) {
	(abs(l+1)^(2*h) + abs(l-1)^(2*h))/2 - l^(2*h)
}

sdk<-function(s,l=1,k=30,h=0.95) {
	s<-sd(s)/(k^(1-h))
	2^0.5 * s * (1-rho(l,h))^0.5
}

x<-readCRU(temps=14)
kou<-function(x,k=5,h=0.95,DS=F,lags=4) {
    cols<-c("green","red","blue","brown")
	xf<-filter(x,filter=rep(1/k,k),sides=1)
	plot(xf,ylim=c(-0.4,0.8),type="l",col="magenta",ylab="Anomaly DegC")
	for (j in 1:lags) {
	d<-NULL
	e<-NULL
	se<-NULL
	s<-sdk(na.omit(x),l=j,k=k,h=h)*2.58
	for (i in 1:(length(xf)-k*j)) {
		d<-c(d,abs(xf[i+k*j]-xf[i]))
		if (DS==T) {
			e<-c(e,abs(xf[length(xf)]-xf[i]))
			len<-length(xf)-i
			se<-c(se,sdk(na.omit(x),l=len,k=k,h=h)*2.58)
		}
	}
		if (DS==T) {
			e<-ts(e,start=tsp(x)[1],frequency=tsp(x)[3])
			lines(e,lty=3,col=cols[j+1])
			se<-ts(se,start=tsp(x)[1],frequency=tsp(x)[3])
			lines(se,lty=5,col=cols[j+1])
		}
	d<-ts(d,start=start(x)[1]+k*j/tsp(x)[3],frequency=tsp(x)[3])
	lines(c(start(x)[1],end(x)[1]),c(s,s),col=cols[j])
	lines(d,lty=2,col=cols[j])
	
	}
}