# R-script for crude examination of the analysis by Scafetta & West 2006 doi:10.1029/2006GL027142

ftp <- 'ftp://ftp.ncdc.noaa.gov/pub/data/paleo/climate_forcing/solar_variability/lean2000_irradiance.txt'

if (!file.exists("Lean2004.txt")) {
  print("Get Lean (2004) from URL")
  Lean2004 <- readLines(ftp)
  first <- grep("YEAR",Lean2004)
  a <- Lean2004[first:length(Lean2004)]
  meta.data <- Lean2004[1:(first-1)]              
  writeLines(a,"Lean2004.txt")
}
Lean2004 <- read.table("Lean2004.txt",header=TRUE)

if (!file.exists("lean1995data.txt")) {
  print("Get Lean et al 1995 from URL")
  Lean1995 <- readLines("ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/lean1995/irradiance_data.txt")
  writeLines(Lean1995,"lean1995data.txt")
  rm(Lean1995)
}
#lean.95 <- read.table("lean1995data.txt",skip=6,
#                      col.names=c("year","TSI"))
Lean1995 <- read.table("lean1995data.txt",skip=9,
                      col.names=c("year","BACKGRN","QS+11-YR ","TSI"))
Lean1995$TSI[Lean1995$TSI < 0] <- NA
tsi.2004 <- data.frame(y=Lean2004$X11yrCYCLE.BKGRND[Lean2004$YEAR>1900],
                       x=Lean2004$YEAR[Lean2004$YEAR>1900])
tsi.1995 <- data.frame(y=Lean1995$TSI[Lean1995$year>1900],
                       x=Lean1995$year[Lean1995$year>1900])
trend.Lean2004 <- lm(y~x,data=tsi.2004)
trend.Lean1995 <- lm(y~x,data=tsi.1995)
x11()
plot(Lean2004$YEAR,Lean2004$X11yrCYCLE.BKGRND,main="TSI",ylim=c(1363,1368),
     sub="Lean (1995,2004)")
grid()
points(Lean1995$year,Lean1995$TSI,col="grey")
lines(Lean2004$YEAR[Lean2004$YEAR>1900],predict(trend.Lean2004),lty=2)
lines(Lean1995$year[Lean1995$year>1900 & is.finite(Lean1995$TSI)],predict(trend.Lean1995),lty=2,col="grey")

mu19 <- c(1365.64, 1364.68,1365.52)
mu18 <- c(1365.48, 1364.44, 1365.43)
T19 <- -0.41
T18 <- -0.49

Z <- round((T19-T18)/(mu19-mu18),2)
print("Values for Z:")
print(Z)

dTSI <- c(max(predict(trend.Lean1995))-min(predict(trend.Lean1995)),
          max(predict(trend.Lean2004))-min(predict(trend.Lean2004)),
          0.75)

print("Estimated changes:")
print(round(dTSI*Z,2))

if (!file.exists("T2m.NH.txt")) {
  print("Get CRU N.H. temperature from URL")
  T2m.NH <- readLines("http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3vnh.txt")
  T2m.NH <- T2m.NH[seq(1,length(T2m.NH),by=2)]  # discard the information about coverage
  writeLines(T2m.NH,"T2m.NH.txt")
  rm(T2m.NH)
}
T2m.NH <- read.table("T2m.NH.txt")
T2m.year <- T2m.NH[,1]
T2m.am <- T2m.NH[,14]

x11()
plot(T2m.year,T2m.am,main="Northern Memisphere Temperature",ylim=c(-0.6,0.8),
     sub="reconstruction based on Lean (1995,2004), SW, and 19th & 18th centuries")
grid()
lines(tsi.1995$x,Z[1]*(tsi.1995$y-tsi.1995$y[1]) + mean(T2m.am[is.element(T2m.year,1890:1910)]),lwd=3,col="red")
lines(tsi.2004$x,Z[1]*(tsi.2004$y-tsi.2004$y[1]) + mean(T2m.am[is.element(T2m.year,1890:1910)]),lwd=3,col="blue")