#Run from command line using R CMD BATCH amo_regression4.r 
#Outputs to fig*.png and amo_regression4.r.Rout
#Run from inside R using source("amo_regression4.r",echo=T)

#Options
realdata=0 #0 for synthetics, 1 for real data.
nsims=10000 #Number of Monte Carlo simulations, for synthetics only.
white_noise=1 #1 for white noise, 0 for autocorrelated ARMA(p,q) noise.
p=1 #ARMA(p,q), only used if white_noise == 0.
q=0 #ARMA(p,q), only used if white_noise == 0.
smooth_amo=1 #1 to LOWESS smooth the AMO index, or 0 to skip smoothing.
#Choose a set of parameters:
#1 - Dr. Tung's new simulation.
#2 - Dumb Scientist's new simulation.
choice=1
if(choice==1){ #Dr. Tung's new simulation.
  global_noise=0.10 #Gaussian noise with this standard deviation for globe.
  atlantic_noise=0.15 #Gaussian noise w/ this st. dev. for N. Atlantic.
} else if(choice==2){ #Dumb Scientist's new simulation.
  global_noise=0.10 #Gaussian noise with this standard deviation for globe.
  atlantic_noise=0.15 #Gaussian noise w/ this st. dev. for N. Atlantic.
} else{
  print("WARNING!!! choice not recognized!")
  doesnt_exist_so_this_will_crash
}
bandlow=50 #Lower period of bandpass filter.
bandhigh=90 #Higher period of bandpass filter.
recent = 1979 #Calculate trends after this year.
xunits="year"
step=0.01 #Stepsize for boxplot of uncertainties vs. trends, in °C/decade.
textsize=1.4
titlesize=1.8
pch1=20 #Line thickness.
par(bg = "white")
desired_width=800
desired_height=800
desired_res=120

#Use real data to define synthetic human, natural influences.
t = 1856:2011 #NOAA's N. Atlantic SST starts in 1856.
n = length(t)
t_zeroed = t - t[1]
if(realdata){
  nsims=1 #Disable Monte-Carlo simulations if real data are loaded.
  qualifier = "Real"
}
if(white_noise == 0){
  library(nlme) #for generalised least squares
}
#Load HadCRUT4 for simulations too.
global = read.table("http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.2.0.0.annual_ns_avg.txt")
#global = read.table("HadCRUT.4.2.0.0.annual_ns_avg.txt")
names(global)[1] = "t"
names(global)[2] = "temp"
keep = c("t","temp")
global = global[keep]
global = global$temp[which(global$t>=t[1] & global$t<=t[n])]
#Only load N. Atlantic SST for real data.
if(realdata){
  n_atlantic = read.table("http://www.esrl.noaa.gov/psd/data/correlation/amon.us.long.mean.data",skip=1,nrows=157) #Skip partial 2013 data.
  #n_atlantic = read.table("amon.us.long.mean.data",skip=1,nrows=157) #Skip partial 2013 data.
  amo = read.table("http://www.esrl.noaa.gov/psd/data/correlation/amon.us.long.data",skip=1,nrows=157) #Skip partial 2013 data.
  #amo = read.table("amon.us.long.data",skip=1,nrows=157) #Skip partial 2013 data.
#Average monthly N. Atlantic SST and monthly AMO index values to get annual means.
  for(i in 1:length(amo[,1])){
    n_atlantic[i,2] = mean(as.numeric(n_atlantic[i,2:length(n_atlantic[i,])]))
    amo[i,2] = mean(as.numeric(amo[i,2:length(amo[i,])]))
  }
  names(n_atlantic)[1] = "t"
  names(n_atlantic)[2] = "temp"
  names(amo)[1] = "t"
  names(amo)[2] = "temp"
  n_atlantic = n_atlantic[keep]
  amo = amo[keep]
  n_atlantic = n_atlantic$temp[which(n_atlantic$t>=t[1] & n_atlantic$t<=t[n])]
  amo = amo$temp[which(amo$t>=t[1] & amo$t<=t[n])]
}
sixth_power_regression = lm(global~t_zeroed+I(t_zeroed^2)+I(t_zeroed^3)+I(t_zeroed^4)+I(t_zeroed^5)+I(t_zeroed^6))
sixth_power = sixth_power_regression$fit
weights=t
for(j in 1:n){
  if(t[j]<=recent){
    weights[j]=0.01
  } else{
    weights[j]=1
  }
}
human = smooth.spline(t,sixth_power,weights,spar=0.85)$y

#Plot global temperatures and fits.
filename_num = 0
filename_num = filename_num+1
filename = sprintf("fig%02d.png",filename_num)
png(file=filename,width=desired_width,height=desired_height,res=desired_res)
title="Global temps, human influence"
yunits="°C"
ylims=c(min(c(global,sixth_power,human)),max(c(global,sixth_power,human)))
plot(t,global,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black",ylim=ylims)
points(t,sixth_power,type="o",pch=pch1,lwd=2,col="blue")
points(t,human,type="o",pch=pch1,lwd=2,col="red")
legend(t[1],ylims[2],c("Global (HadCRUT4) annual mean", "A 6th order polynomial trend","Spline-smoothed monotonic human influence"), pch = c(pch1,pch1,pch1), col = c("black","blue","red"), lty = c(1,1,1))
dev.off

#Plot derivative of synthetic human influence.
filename_num = filename_num+1
filename = sprintf("fig%02d.png",filename_num)
png(file=filename,width=desired_width,height=desired_height,res=desired_res)
title="Derivative of human influence"
yunits="°C/year"
human_diff = diff(human)/diff(t)
ylims=c(0,max(human_diff))
t_diff = seq(t[1],t[n-1])
plot(t_diff,human_diff,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black",ylim=ylims)
if(min(human_diff) >= 0){
  text(t_diff[1],ylims[2]*.9,pos=4,"Monotonic.",cex=2,col="green")
} else{
  text(t_diff[1],ylims[2]*.9,pos=4,"NOT MONOTONIC!",cex=2,col="red")
}
dev.off

#Bandpass filter nature timeseries.
nature = global - human
nature_fft = fft(nature)
nyquist = 1/2/(t[2]-t[1])
if((length(t) %% 2) == 1){# Odd number of samples
  f = c(seq(0, nyquist, length.out=(length(t)+1)/2), seq(-nyquist, 0, length.out=(length(t)-1)/2))
} else{#Even number
  f = c(seq(0, nyquist, length.out=length(t)/2), seq(-nyquist, 0, length.out=length(t)/2))
}
for(j in 1:length(f)){
  if(abs(f[j]) < 1/bandhigh | abs(f[j]) > 1/bandlow){
    nature_fft[j] = 0.0
  }
}
nature = Re(fft(nature_fft,inverse=T)/length(nature_fft))

#Plot natural influence.
filename_num = filename_num+1
filename = sprintf("fig%02d.png",filename_num)
png(file=filename,width=desired_width,height=desired_height,res=desired_res)
title="Detrended global temps, natural influence"
yunits="°C"
plot(t,global-human,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black")
points(t,nature,type="o",pch=pch1,lwd=2,col="green")
legend(t[1],max(global-human),c("Global (HadCRUT4) annual mean - human influence", "Smoothed natural influence"), pch = c(pch1,pch1), col = c("black","green"), lty = c(1,1))
dev.off

#Subtract human and natural influences to obtain real noise.
real_noise = global - human - nature
real_noise = real_noise - mean(real_noise)
print(var(real_noise))
print(sqrt(var(real_noise)))
filename_num = filename_num+1
filename = sprintf("fig%02d.png",filename_num)
png(file=filename,width=desired_width,height=desired_height,res=desired_res)
title="Remaining real global noise"
yunits="°C"
plot(t,real_noise,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black")
dev.off
#Analyze real noise.
autonoise = arima(real_noise,order=c(p,0,q))
print(autonoise)
if(p) ar_coef = matrix(data=NA,nrow=p)
if(q) ma_coef = matrix(data=NA,nrow=q)
#Extract noise coefficients from autonoise. AR first, then MA.
if(p) for(j in 1:p){
  ar_coef[j] = autonoise$coef[j]
}
if(q) for(j in 1:q){
  ma_coef[j] = autonoise$coef[j+p]
}
if(p) print(ar_coef)
if(q) print(ma_coef)

#Calculate recent trends.
recent_t = t[which(t>recent)]
recent_human = human[which(t>recent)]
recent_human_trend = lm(recent_human~recent_t)
true_summary = summary(recent_human_trend)
true_human_trend = coef(true_summary)[2,1]
trends = matrix(data=0,nrow=nsims,ncol=1)
trenderrors = matrix(data=0,nrow=nsims,ncol=1)
quadratics = matrix(data=0,nrow=nsims,ncol=1)
quadraticerrors = matrix(data=0,nrow=nsims,ncol=1)
all_values = matrix(data=0,nrow=3,ncol=1)
error_bars = matrix(data=0,nrow=3,ncol=1)#2 sigma
global_atlantic_corr = matrix(data=0,nrow=nsims,ncol=1)
global_var = matrix(data=0,nrow=nsims,ncol=1)
atlantic_var = matrix(data=0,nrow=nsims,ncol=1)

#Run Monte Carlo simulations.
successes=0#Number of times the true post-1979 trend is within error bars.
for(i in 1:nsims){

  #Generate synthetic global and N. Atlantic autocorrelated noise
  #with the same ARMA coefficients as the real detrended global data.
  if(!white_noise & !realdata){
    if(p & !q){#AR
      global_autonoise = arima.sim(list(ar=ar_coef),n=n)
      atlantic_autonoise = arima.sim(list(ar=ar_coef),n=n)
    } else if(!p & q){#MA
      global_autonoise = arima.sim(list(ma=ma_coef),n=n)
      atlantic_autonoise = arima.sim(list(ma=ma_coef),n=n)
    } else if(p & q){#ARMA
      global_autonoise = arima.sim(list(ar=ar_coef,ma=ma_coef),n=n)
      atlantic_autonoise = arima.sim(list(ar=ar_coef,ma=ma_coef),n=n)
    } else{
      print("WARNING!!! If p and q are both 0, just use white noise!")
      doesnt_exist_so_this_will_crash
    }
    #Adjust global variance to match real global noise.
    global_autonoise = global_autonoise/sqrt(var(global_autonoise)/var(real_noise))
    #Adjust N. Atlantic variance to match real global noise times requested atlantic/global.
    atlantic_autonoise = atlantic_autonoise/sqrt(var(atlantic_autonoise)/var(real_noise))*atlantic_noise/global_noise
    if(i==1){
      filename_num = filename_num+1
      filename = sprintf("fig%02d.png",filename_num)
      png(file=filename,width=desired_width,height=desired_height,res=desired_res)
      title="Synthetic noises"
      yunits="°C"
      #Plot using same scale as real noise for easier comparison.
      ylims=c(min(real_noise),max(real_noise))
      plot(t,global_autonoise,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black",ylim=ylims)
      points(t,atlantic_autonoise,type="o",pch=pch1,lwd=2,col="red")
      legend(t[1],ylims[2],c("Synthetic global noise","Synthetic N. Atlantic SST noise"),pch = c(pch1,pch1), col = c("black","red"), lty = c(1,1))
      dev.off
      filename_num = filename_num+1
      filename = sprintf("fig%02d.png",filename_num)
      png(file=filename,width=desired_width,height=desired_height,res=desired_res)
      par(mfrow = c(2, 2))
      acf(real_noise,main="ACF of real global noise")
      pacf(real_noise,main="PACF of real global noise")
      acf(global_autonoise,main="ACF of synthetic global noise")
      pacf(global_autonoise,main="PACF of synthetic global noise")
      dev.off
    }
  }

  #Define synthetic global, N. Atlantic temperature anomalies.
  if(!realdata){
    qualifier = "Synthetic"
    if(white_noise){
      global     =     human +     nature + rnorm(t,mean=0,sd=global_noise)
      n_atlantic = 0.8*human + 1.4*nature + rnorm(t,mean=0,sd=atlantic_noise)
    } else{
      global     =     human +     nature + global_autonoise
      n_atlantic = 0.8*human + 1.4*nature + atlantic_autonoise
    }
  }
  #Remove means so real and synthetic datasets all have exactly the same mean: 0.
  global = global - mean(global)
  n_atlantic = n_atlantic - mean(n_atlantic)
  global_atlantic_corr[i] = cor(global,n_atlantic)
  global_var[i] = var(global)
  atlantic_var[i] = var(n_atlantic)
  #Plot synthetic global, N. Atlantic temperature anomalies.
  if(i==1){
    filename_num = filename_num+1
    filename = sprintf("fig%02d.png",filename_num)
    png(file=filename,width=desired_width,height=desired_height,res=desired_res)
    title="surface temperature anomalies"
    title=paste(qualifier, title, collapse = "") 
    #ylims=c(min(c(global,n_atlantic)),max(c(global,n_atlantic)))
    #Fix plot limits so real and synthetic plots are easier to compare.
    ylims=c(-0.7,0.7)
    plot(t,global,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black",ylim=ylims)
    points(t,n_atlantic,type="o",pch=pch1,lwd=2,col="red")
    if(realdata){
      legends = c("Global (HadCRUT4) annual mean", "N. Atlantic SST (NOAA) annual mean")
    }
    else{
      legends = c("Global", "N. Atlantic")
    }
    legend(t[1],ylims[2],legends, pch = c(pch1,pch1), col = c("black","red"), lty = c(1,1))
    dev.off
  }
  #Define linearly-detrended synthetic AMO (real AMO downloaded separately).
  if(!realdata){
    n_atlantic_trend = lm(n_atlantic~t)
    amo = n_atlantic - coef(summary(n_atlantic_trend))[2,1]*(t-t[1])
  }
  if(smooth_amo){
    smoothed_amo = lowess(t,amo,f=25/n)$y
    #smoothed_amo = loess(amo~t,data=data.frame(t,amo),span=1)
  }
  #Plot detrended AMO.
  if(i==1){
    filename_num = filename_num+1
    filename = sprintf("fig%02d.png",filename_num)
    png(file=filename,width=desired_width,height=desired_height,res=desired_res)
    title="linearly-detrended AMO index"
    title=paste(qualifier, title, collapse = "") 
    ylims=c(min(amo),max(amo))
    plot(t,amo,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black",ylim=ylims)
    if(smooth_amo){
      points(t,smoothed_amo,type="o",pch=pch1,lwd=2,col="blue")
      legend(t[1],ylims[2],c("Original AMO index", "Smoothed AMO index"), pch = c(pch1,pch1), col = c("black","blue"), lty = c(1,1))
    }
    dev.off
  }
  if(smooth_amo){
    amo = smoothed_amo
  }
  #Regress global temperatures against human influence and AMO.
  human_p = human - mean(human) #Remove mean so it's handled by the intercept.
  amo_p = amo - mean(amo)
  nature_p = nature - mean(nature)
  if(white_noise == 0){
    regression = gls(global~human_p+amo_p,corr=corARMA(p=p,q=q))
    regression_summary = summary(regression)
    all_values[1] = all_values[1] + regression_summary$tTable[1,1]
    all_values[2] = all_values[2] + regression_summary$tTable[2,1]
    all_values[3] = all_values[3] + regression_summary$tTable[3,1]
    error_bars[1] = error_bars[1] + 2*regression_summary$tTable[1,2]#2 sigma
    error_bars[2] = error_bars[2] + 2*regression_summary$tTable[2,2]#2 sigma
    error_bars[3] = error_bars[3] + 2*regression_summary$tTable[3,2]#2 sigma
    est_human = regression_summary$tTable[2,1]*human_p
  }
  else{
    regression = lm(global~human_p+amo_p)
    regression_summary = summary(regression)
    all_values[1] = all_values[1] + coef(regression_summary)[1,1]
    all_values[2] = all_values[2] + coef(regression_summary)[2,1]
    all_values[3] = all_values[3] + coef(regression_summary)[3,1]
    error_bars[1] = error_bars[1] + 2*coef(regression_summary)[1,2]#2 sigma
    error_bars[2] = error_bars[2] + 2*coef(regression_summary)[2,2]#2 sigma
    error_bars[3] = error_bars[3] + 2*coef(regression_summary)[3,2]#2 sigma
    est_human = coef(regression_summary)[2,1]*human_p
  }
  #print(regression_summary)
  #Plot residuals.
  if(i==1){
    filename_num = filename_num+1
    filename = sprintf("fig%02d.png",filename_num)
    png(file=filename,width=desired_width,height=desired_height,res=desired_res)
    if(smooth_amo){
      title="Residuals of fit to human, smoothed AMO"
    }
    else{
      title="Residuals of fit to human, unsmoothed AMO"
    }
    residuals = resid(regression)
    plot(t,residuals,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black")
    dev.off
  }
  #Add residual back.
  est_human2 = est_human + residuals
  recent_est_human2 = est_human2[which(t>recent)]
  if(white_noise == 0){
    recent_est_human2_trend = gls(recent_est_human2~recent_t,corr=corARMA(p=p,q=q))
    est_summary2 = summary(recent_est_human2_trend)
    est_human2_trend = est_summary2$tTable[2,1]
    est_human2_error = 2*est_summary2$tTable[2,2]#2 sigma
  }
  else{
    recent_est_human2_trend = lm(recent_est_human2~recent_t)
    est_summary2 = summary(recent_est_human2_trend)
    est_human2_trend = coef(est_summary2)[2,1]
    est_human2_error = 2*coef(est_summary2)[2,2]#2 sigma
  }
  true_blurb = sprintf("True human. Post-%d trend = %+.3f%s/decade",recent,true_human_trend*10,yunits) 
  est_blurb2 = sprintf("Est. human. Post-%d trend = %+.3f ± %.3f",recent,est_human2_trend*10,est_human2_error*10) 
  if(true_human_trend <= est_human2_trend+est_human2_error & true_human_trend >= est_human2_trend-est_human2_error) successes = successes + 1
  #Plot true and estimated human influences.
  if(i==1){
    recent_trendline = est_human2_trend*(recent_t-recent_t[1])
    recent_trendline = recent_trendline - mean(recent_trendline) + mean(recent_est_human2)
    #Shift true human so its recent mean equals the mean of the estimated 
    #trendline, for easier comparison of recent trends.
    human = human - (mean(recent_human) - mean(recent_trendline))
    filename_num = filename_num+1
    filename = sprintf("fig%02d.png",filename_num)
    png(file=filename,width=desired_width,height=desired_height,res=desired_res)
    title="True, estimated human influences"
    ylims=c(min(c(human,est_human2)),max(c(human,est_human2)))
    plot(t,human,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col="black",ylim=ylims)
    points(t,est_human2,type="o",pch=pch1,lwd=2,col="red")
    points(recent_t,recent_trendline,type="l",pch=pch1,lwd=6,lty=3,col="blue")
    lowerbound=recent_trendline-est_human2_error*recent_t
    lowerbound=lowerbound - mean(lowerbound) + mean(recent_trendline)
    points(recent_t,lowerbound,type="l",lwd=6,lty=1,col="blue")
    upperbound=recent_trendline+est_human2_error*recent_t
    upperbound=upperbound - mean(upperbound) + mean(recent_trendline)
    points(recent_t,upperbound,type="l",lwd=6,lty=1,col="blue")
    legend(t[1],ylims[2],c(true_blurb, est_blurb2), pch = c(pch1,pch1), col = c("black","red"), lty = c(1,1))
    dev.off
    #Shift true human back to original baseline, just in case it's needed later.
    human = human + (mean(recent_human) - mean(recent_trendline))
  }
  trends[i] = est_human2_trend*10
  trenderrors[i] = est_human2_error*10
  #Calculate estimated quadratic term.
  if(white_noise == 0){
    est_human2_quadratic = gls(est_human2~t+I(t^2),corr=corARMA(p=p,q=q))
    est_quadratic = summary(est_human2_quadratic)$tTable[3,1]
    est_quadraticerror = summary(est_human2_quadratic)$tTable[3,2]
  }
  else{
    est_human2_quadratic = lm(est_human2~t+I(t^2))
    est_quadratic = coef(summary(est_human2_quadratic))[3,1]
    est_quadraticerror = coef(summary(est_human2_quadratic))[3,2]
  }
  quadratics[i] = est_quadratic
  quadraticerrors[i] = est_quadraticerror
}

#Report results.
all_values = all_values/nsims
error_bars = error_bars/nsims
print(all_values)
print(error_bars)
mean(global_atlantic_corr)
mean(global_var)
mean(atlantic_var)
mean(trends)
mean(trenderrors)
true_human_quadratic = lm(human~t+I(t^2))
true_quadratic = coef(summary(true_human_quadratic))[3,1]
quadratic_blurb = sprintf("True quadratic = %.3e",true_quadratic) 
print(quadratic_blurb)
print(successes)
successes_blurb = sprintf("The true post-1979 trend was within the estimated 95%% confidence interval %.3f%% of the time.",successes/nsims*100) 
print(successes_blurb)
if(!realdata){
  library(MASS) #for fitdistr, to determine mean and 95% C.I. of distributions.
  this_dist = fitdistr(trends,"normal")
  this_estimate = this_dist$estimate[1]
  this_error = 2*this_dist$estimate[2]#2 sigma
  filename_num = filename_num+1
  filename = sprintf("fig%02d.png",filename_num)
  png(file=filename,width=desired_width,height=desired_height,res=desired_res)
  title = sprintf("Est. post-%d human trend:\n %+.3f ± %.3f%s/decade",recent,this_estimate,this_error,yunits)
  xunits = "°C/decade"
  hist(trends,density=20,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,main=title)
  dev.off
  this_dist = fitdistr(quadratics,"normal")
  this_estimate = this_dist$estimate[1]
  this_error = 2*this_dist$estimate[2]#2 sigma
  filename_num = filename_num+1
  filename = sprintf("fig%02d.png",filename_num)
  png(file=filename,width=desired_width,height=desired_height,res=desired_res)
  title = sprintf("Est. post-%d human quadratic term:\n %+.2e ± %.2e%s/year^2",t[1],this_estimate,this_error,yunits)
  xunits = "°C/year^2"
  hist(quadratics,density=20,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,main=title)
  dev.off
  this_dist = fitdistr(global_atlantic_corr,"normal")
  this_estimate = this_dist$estimate[1]
  this_error = 2*this_dist$estimate[2]#2 sigma
  filename_num = filename_num+1
  filename = sprintf("fig%02d.png",filename_num)
  png(file=filename,width=desired_width,height=desired_height,res=desired_res)
  title = sprintf("Correlation b/w global, N. Atlantic temps:\n %.3f ± %.3f",this_estimate,this_error)
  xunits = "correlation"
  hist(global_atlantic_corr,density=20,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,main=title)
  dev.off
  this_dist = fitdistr(global_var,"normal")
  this_estimate = this_dist$estimate[1]
  this_error = 2*this_dist$estimate[2]#2 sigma
  filename_num = filename_num+1
  filename = sprintf("fig%02d.png",filename_num)
  png(file=filename,width=desired_width,height=desired_height,res=desired_res)
  title = sprintf("Variance of global temps:\n %.3f ± %.3f%s^2",this_estimate,this_error,yunits)
  xunits = "°C^2"
  hist(global_var,density=20,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,main=title)
  dev.off
  this_dist = fitdistr(atlantic_var,"normal")
  this_estimate = this_dist$estimate[1]
  this_error = 2*this_dist$estimate[2]#2 sigma
  filename_num = filename_num+1
  filename = sprintf("fig%02d.png",filename_num)
  png(file=filename,width=desired_width,height=desired_height,res=desired_res)
  title = sprintf("Variance of N. Atlantic temps:\n %.3f ± %.3f%s^2",this_estimate,this_error,yunits)
  xunits = "°C^2"
  hist(atlantic_var,density=20,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,main=title)
  dev.off
  filename_num = filename_num+1
  filename = sprintf("fig%02d.png",filename_num)
  png(file=filename,width=desired_width,height=desired_height,res=desired_res)
  title = sprintf("Post-%d trend uncertainties vs. trends",recent) 
  xunits = "Trend - °C/decade"
  yunits = "95% Uncertainties - °C/decade"
  #step=(max(trends)-min(trends))/10
  groups=floor(trends/step)*step
  boxplot(trenderrors~groups,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,main=title,xlab=xunits,ylab=yunits)
  dev.off
}
print(successes_blurb)