##VARIATION WITH "CORRECT' DATA

library(fields)
library(signal)
library(xlsReadWrite)
 	source("d:/climate/scripts/mann.2008/utilities.v6.txt")
  	period_MBH=seq(1800,200,-100);M=length(period_MBH);M
  	use0="pairwise.complete.obs"	

	url="http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/cps-validation.xls"
	download.file(url,"temp.xls",mode="wb")
	cps.info=scan("http://www.climateaudit.info/data/mann.2008/statistics/cps.info.dat",what="")

	load("d:/climate/scripts/mann.2008/notes/Real_tornetrask.tab")
	load("d:/climate/scripts/mann.2008/notes/Notilj.tab")


#GDD (file="d:/climate/images/2010/mann08/aug8_figure1.png",type="png",h=480, w=560)
layout(1) 
	par(mar=c(3,4,2,1))
	stat.mbh=get.stat.mbh(paste(Target$hemi[j],"full",Target$target[j],sep="_"))
	plot(stat.mbh$century, stat.mbh$earlym_RE,type="l",col=2,xlab="Year AD",ylab="Mannian RE",lty=3,ylim=c(-1,1))
	lines(stat.mbh$century, stat.mbh$latem_RE,lty=3,col=4)
	A=Notilj$stat[[j]]
	points(A$century, A$earlym_RE,pch=3,col=2)
	points(A$century, A$latem_RE,pch=3,col=4)
	B=Real_tornetrask$stat[[j]]
	points(B$century, B$earlym_RE,pch=19,col=2,cex=.8)
	points(B$century, B$latem_RE,pch=19,col=4,cex=.8)
	title("Mannian RE Statistics: realdata")
	legend("bottomright",lty=c(3,-1,-1),pch=c(-1,3,19),merge=TRUE,legend=c("M08", "Notilj","realdata") )
	legend("bottomleft",fill=c(2,4),legend=c("earlym","latem") )
	
#dev.off()


#GDD (file="d:/climate/images/2010/mann08/aug8_figure2.png",type="png",h=480, w=560)
layout(1) 
j=1
	par(mar=c(3,4,2,1))
	stat.mbh=get.stat.mbh(paste(Target$hemi[j],"full",Target$target[j],sep="_"))
	plot(stat.mbh$century, stat.mbh$avg_RE,type="l",col=1,xlab="Year AD",ylab="Mannian RE",lty=3,ylim=c(-1,1))
	A=Notilj$stat[[j]]
	points(A$century, A$avg_RE,pch=3,col=3)
	B=Real_tornetrask$stat[[j]]
	points(B$century, B$avg_RE,pch=19,col=2,cex=.8)
	abline(h=.36,lty=3)
	lines(stat.mbh$century, stat.mbh$adj_RE)
	title("M08 RE: the realdata version")
#dev.off()









GDD (file="d:/climate/images/2010/mann08/aug8_figure1a.gif",type="gif",h=560, w=480)

layout(array(1:4,dim=c(2,2)) ,heights=c(1.1,1.3),widths=c(1.25,1.1) )

par(mar=c(0,4,2,0))	
	case="earlym";hemi="NH"
	X= cps[["CRU"]][[9]]
	instr.smooth=Basics$instr.smooth
	plot(X[[case]][,"NH"],type="l",xlim=c(900,2050),xaxs="i",ylim=c(-.9,1) ,xlab="",ylab="deg C",axes=FALSE)
	axis(side=1,labels=FALSE); axis(side=2); box(); abline(h=0,lty=3)
	lines(instr.smooth[,"CRU_NH"],lwd=1,col=2)
	lines(X[[case]][,"NH"])
	text(900,.9,paste(case,hemi,A[hemi,paste(case,"RE",sep="_")]),font=2,pos=4)
	
par(mar=c(3,4,0,0))	
	case="latem";hemi="NH"
	plot(X[[case]][,"NH"],type="l",xlim=c(900,2050),xaxs="i",ylim=c(-.9,1) ,xlab="",ylab="deg C",axes=FALSE)
	axis(side=1); axis(side=2,labels=FALSE); box(); abline(h=0,lty=3)
	lines(instr.smooth[,"CRU_NH"],lwd=1,col=2)
	text(900,.9,paste(case,hemi,A[hemi,paste(case,"RE",sep="_")]),font=2,pos=4)
	lines(X[[case]][,"NH"])

par(mar=c(0,0,2,1))	
	case="earlym";hemi="SH"
	plot(X[[case]][,"SH"],type="l",xlim=c(900,2050),xaxs="i",ylim=c(-.9,1) ,xlab="",ylab="",axes=FALSE)
	axis(side=1,labels=FALSE); axis(side=2,labels=FALSE); box(); abline(h=0,lty=3)
	lines(instr.smooth[,"CRU_SH"],lwd=1,col=2)
	text(900,.9,paste(case,hemi,A[hemi,paste(case,"RE",sep="_")]),font=2,pos=4)
	lines(X[[case]][,"SH"])
	mtext(side=3,"AD1000 Mann et al 2008 'Standard'",font=2,cex=1.2,line=0.5)

par(mar=c(3,0,0,1))	
	case="latem";hemi="SH"	
	plot(X[[case]][,"SH"],type="l",xlim=c(900,2050),xaxs="i",ylim=c(-.9,1) ,xlab="",ylab="",axes=FALSE)
	axis(side=1); axis(side=2,labels=FALSE); box(); abline(h=0,lty=3)
	lines(instr.smooth[,"CRU_SH"],lwd=1,col=2)
	text(900,.9,paste(case,hemi,A[hemi,paste(case,"RE",sep="_")]),font=2,pos=4)
	lines(X[[case]][,"SH"])

dev.off()

	load("d:/climate/scripts/mann.2008/notes/Notilj.tab")
	case="earlym";hemi="NH"
	X= cps[["CRU"]][[9]]
	case="latem";hemi="NH"
	plot(X[[case]][,"NH"],type="l",xlim=c(900,2050),xaxs="i",ylim=c(-.9,1) ,xlab="",ylab="deg C",axes=FALSE)
	axis(side=1); axis(side=2,labels=FALSE); box(); abline(h=0,lty=3)
	lines(instr.smooth[,"CRU_NH"],lwd=1,col=2)
	text(900,.9,paste(case,hemi,A[hemi,paste(case,"RE",sep="_")]),font=2,pos=4)

	lines(Notilj$cps[["CRU"]][[9]][[case]][,"NH"],col=2,lty=3)

		
	par(mar=c(3,4,0,0))	
	case="latem";hemi="NH"
	plot(X[[case]][,"NH"],type="l",xlim=c(900,2050),xaxs="i",ylim=c(-.9,1) ,xlab="",ylab="deg C",axes=FALSE)
	axis(side=1); axis(side=2,labels=FALSE); box(); abline(h=0,lty=3)
	lines(instr.smooth[,"CRU_NH"],lwd=1,col=2)
	text(900,.9,paste(case,hemi,A[hemi,paste(case,"RE",sep="_")]),font=2,pos=4)
	lines(X[[case]][,"NH"])



k=9
case="earlym"	
par(mar=c(0,1,2,1))	
	X=emulation[[k]][[case]]$unique
	Y=emulation[[k]][[case]]$proxybase
	Z=emulation[[k]][[case]]$proxydata
	Y$long=details$long[Y$number]
	X$flip= Z$flip[match(X$unique,Z$unique)]
	X$col=2; X$col[X$lat<0]=3;X$col[X$flip]=4
	plot(X$long,X$lat,col=X$col,axes=FALSE, cex=6*sqrt(X$weight),xlab="",ylab="",pch=19,xlim=c(-175,175),ylim=c(-70,85))
	text(-180,80,case,font=2,pos=4,cex=1.2)
	axis(side=1,labels=FALSE); axis(side=2);box(); abline(h=0,lty=3)
		world(add=TRUE,col=5)#"limegreen")
	points(Y$long[!Y$pass],Y$lat[!Y$pass],pch="+",col=1)


case="latem"	
par(mar=c(3,1,0,1))	
	X=emulation[[k]][[case]]$unique
	Y=emulation[[k]][[case]]$proxybase
	Z=emulation[[k]][[case]]$proxydata
	Y$long=details$long[Y$number]
	X$flip= Z$flip[match(X$unique,Z$unique)]
	X$col=2; X$col[X$lat<0]=3;X$col[X$flip]=4
	plot(X$long,X$lat,col=X$col,axes=FALSE, cex=6*sqrt(X$weight),xlab="",ylab="",pch=19,xlim=c(-175,175),ylim=c(-70,85))
	text(-180,80,case,font=2,pos=4,cex=1.2)
	axis(side=1); axis(side=2);box(); abline(h=0,lty=3)
	world(add=TRUE,col=5)#"limegreen")
		points(Y$long[!Y$pass],Y$lat[!Y$pass],pch="+",col=1)

dev.off()	

dev.off()