# Some preamble setup and read in files

source("AntLocal-Functions.r")

# Read in proxy data
allproxies<-read.ts("AntarcticProxyMatrix.csv",sep=";",header=T)
# Columns are
# 1 Talos firn core
# 2 Law Dome DSS
# 3 Plateau Remote w infill in calib
# 4 Ind22/B4
# 5 EDML
# 6 WAIS Divide 2005A
# 7 WAIS Divide 2006A with infill from 05A in calib and from 1585
# 8 ITASE 00-01
# 9 ITASE 00-05
# 10 James Ross Is
# 11 Siple Station with infill in calib

# With targets used here, the IND series is negatively correlated and weighted sum in R
# needs positive correlation, so invert it in case someone decides to play with it.
invertIND<-allproxies[,4]*-1
allproxies[,4]<-invertIND


Proxynet.fullant<-cbind(allproxies[,c(1,2,3,5,6,7,8,9,11)])
Proxynet.ea<-cbind(allproxies[,c(1,2,3,5)])
Proxynet.wa<-cbind(allproxies[,c(6,7,8,9,10,11)])


# Steig target comes as Antarctica;WAIS+Pen;EAIS;WAIS;Pen
target<-read.ts("Ant_instrum_target.csv",sep=";",header=T)
target.fullant<-target[,1]
target.ea<-target[,3]
target.wa<-target[,2]

# The WIAS WD2006A record finishes in 1584 and is infilled for normalisation and calibration using
# a smoothed version of WD2005A above this period. This allows for standardisation, calibration, verification
# Once these operations are performed, the smoothed infill is not used in the actual reconstruction outside
# of the calibration/verficiation.

# The relevant column in the series for the WD2006 record is column 6 for the fullant proxy network
# and column 2 in the West Antarctic network


# Do the full antarctic reconstruction, noting WD06 is present and in column 6
outdata <- dorecon(Proxynet.fullant,target.fullant,T,6)
label<-"FullAnt"

# Write it out to file
recfile<-paste("recon_",label,".csv",sep="")
colhead1<-paste(label,"Recon")
colhead2<-paste(label,"Recon decadal")
write.table(outdata,recfile,sep=";",col.names=
  c("Year",colhead1,"2SE",colhead2,"2SE_decadal","REs_fullcal","REs_verif1","REs_verfi2","Mean_RE"),row.names=F)

# Do the west antarctic reconstruction, noting WD06 is present and in column 2
outdata <- dorecon(Proxynet.wa,target.wa,T,2)
label<-"WestAnt"
WArecon<-outdata

# Write it out to file
recfile<-paste("recon_",label,".csv",sep="")
colhead1<-paste(label,"Recon")
colhead2<-paste(label,"Recon decadal")
write.table(outdata,recfile,sep=";",col.names=
  c("Year",colhead1,"2SE",colhead2,"2SE_decadal","REs_fullcal","REs_verif1","REs_verfi2","Mean_RE"),row.names=F)

# Do the east antarctic reconstruction, noting WD06 is absent and column is irrelevant
outdata <- dorecon(Proxynet.ea,target.ea,F,6)
label<-"EastAnt"
EArecon<-outdata

# Write it out to file
recfile<-paste("recon_",label,".csv",sep="")
colhead1<-paste(label,"Recon")
colhead2<-paste(label,"Recon decadal")
write.table(outdata,recfile,sep=";",col.names=
  c("Year",colhead1,"2SE",colhead2,"2SE_decadal","REs_fullcal","REs_verif1","REs_verfi2","Mean_RE"),row.names=F)

# Grid cells in target
# EAntarctica 4106
# Wantarctica+Peninsula 1403
EAwt<-4106/(4106+1403)
WAwt<-1403/(4106+1403)

# Area weighted reconstruction
compstart<-max(min(time(EArecon)),min(time(WArecon)))
compend<-min(max(time(EArecon)),max(time(WArecon)))

wtdsum<-window(EArecon[,c(2)],start=compstart,end=compend)*EAwt+window(WArecon[,c(2)],start=compstart,end=compend)*WAwt

outdata<-cbind(time(wtdsum),wtdsum)
write.table(outdata,"recon_EAplusWA.csv",sep=";",col.names=c("Year","WtdEA+WA"),row.names=F)

######################################################################
# end of reconstructions
###################################################################################