# from gaussian fit for fossil fuel production
# derived from 'proven reserves' = -2 sig
# derived from 'optimistic reserves' = +1 sig
co2_gt_ff_future_mu = 6100 # gigatonnes CO2
co2_gt_ff_future_sig = 1700 # gigatonnes CO2

# airborne fraction is ~50% emission
# 60%: http://www.ipcc.ch/publications_and_data/ar4/wg1/en/ch2s2-3.html#2-3-1
# 55%: http://www.ipcc.ch/publications_and_data/ar4/wg1/en/figure-7-4.html
# 45%: http://www.skepticalscience.com/news.php?n=112
abfrac = 0.5
 
# Past Fossil Fuel Consumption
# CDIAC FAQ
# http://cdiac.ornl.gov/pns/faq.html
# C from fossil fuels 1850-2000: 282 GtC
# C02 from fossil fuels 1850-2000: 282 * 3.67 ~= 1000 GtCO2
# leads to ~ 65 ppm CO2 (50% airborne frac)
co2_gt_ff_past = 1000

# Past Non-Fossil Fuel CO2 Components
# CDIAC FAQ
# http://cdiac.ornl.gov/pns/faq.html
# Land use is 154GtC and cement is 6GtC
# other = 160 * 3.67 ~ 600 gigatones
# leads to ~ 35 ppm
co2_gt_other_past = 600
co2_ppm_preindustrial = 280
co2_ppm_other = 35

# see archer 2012 Real Climate
# http://www.realclimate.org/index.php/archives/2012/01/much-ado-about-methane/
# assumes 1000 Gt contribution from CH4 (does airborne frac apply here? if not, 134 ppm)
# guess a worst-case 1000 Gton from hydrates
# after 2000 Gton C (*3.67 = CO2)
co2_ppm_methane = 0
#co2_ppm_methane = 65 
#co2_ppm_methane = 200 #

# CO2 ppm per kg
# CDIAC FAQ
# http://cdiac.ornl.gov/pns/faq.html
# 3.67CO2/C *  2.13 GtC/PPM
gt_ppm = 7.81
 
 # grab default graphic parameters
def.par = par(no.readonly = TRUE) 

# flag for controlling image output
graphic_output = TRUE

# -----------------------------------------------------------------------------
# Gaussian for Future Fossil Fuel Gigatonnes Availability in Gigatonnes CO2
# -----------------------------------------------------------------------------
mu = co2_gt_ff_future_mu
sig = co2_gt_ff_future_sig
co2_gt_ff_future=0:(2*mu) # x scale, gigatonnes C02

# gaussian
co2_gt_ff_future_pdf=(1/(sig*sqrt(2*pi)))*exp(-0.5*((co2_gt_ff_future-mu)/sig)^2)

# -----------------------------------------------------------------------------
# Gaussian for Resulting Atmospheric CO2 Concentration from Fossil Fuel
# (past and future)
# -----------------------------------------------------------------------------
mu = co2_gt_ff_future_mu
sig = co2_gt_ff_future_sig
mu_index = mu

# calculate CO2 PPM for all fossil fuels, past and future
co2_ppm_ff_all = (co2_gt_ff_future + co2_gt_ff_past) * abfrac / gt_ppm

# find mu and sig using fact that this is same distribution as FF Gt CO2 above
co2_ppm_ff_all_mu = co2_ppm_ff_all[1+as.integer(mu_index)]
co2_ppm_ff_all_sig = co2_ppm_ff_all[1+as.integer(mu+sig)]-co2_ppm_ff_all[1+as.integer(mu)] 
mu = co2_ppm_ff_all_mu
sig = co2_ppm_ff_all_sig

# the pdf is the same as previous
co2_ppm_ff_all_pdf = co2_gt_ff_future_pdf 

# -----------------------------------------------------------------------------
# Gaussian for Cumulative Atmospheric CO2 Concentration (all sources)
# -----------------------------------------------------------------------------
# change in forcing / 280 for preindustral / 35 for land use
co2_ppm_all = co2_ppm_ff_all + co2_ppm_preindustrial + co2_ppm_other + co2_ppm_methane

dF = 5.35 * log((co2_ppm_all)/(co2_ppm_preindustrial + co2_ppm_other))
# ecs/3.7
lam15 = 1.5/3.7
lam20 = 2.0/3.7
lam25 = 2.5/3.7
lam30 = 3.0/3.7
lam35 = 3.5/3.7
lam40 = 4.0/3.7
lam45 = 4.5/3.7
# calculate equilibrium change in temp as function of CO2 forcings and ecs
dT15 = lam15 * dF
dT20 = lam20 * dF
dT25 = lam25 * dF
dT30 = lam30 * dF
dT35 = lam35 * dF
dT40 = lam40 * dF
dT45 = lam45 * dF
 
# -----------------------------------------------------------------------------
# ECS distribution: Poisson distribution
# -----------------------------------------------------------------------------
# mean = climate sensitivity = 3C
lambda = 2.5 # ECS mean = 2.5
xmax = 14 # long enough tail to go to zero
xx = 0:xmax # for Possion, use integers
xx2 = xx/2 # scale the Poisson back down to 0:7
yy = lambda^xx * exp(-lambda) / factorial(xx) # Poisson
Ts_poisson = approx(xx2,yy,n=71) # need more points, so interpolate to dx = 0.1
shift = 17 # shift peak back to ECS = 3
yy2 = rep(0,71) 
yy2[shift:71] = Ts_poisson$y[1:(71-shift+1)]
Ts_poisson$y = yy2 / 5 # shift complete
sum(Ts_poisson$y[Ts_poisson$x>=2 & Ts_poisson$x<=4.5])
#[1] 0.8526035

# -----------------------------------------------------------------------------
# ECS distribution: Gamma distribution
# -----------------------------------------------------------------------------
# mean = climate sensitivity ~= 3C
#alpha = 12			# width of dist, higher is narrower
#beta = alpha/3.25 	# shift, higher is right shift
alpha = 9			# width of dist, higher is narrower
beta = alpha/2.8 	# shift, higher is right shift
x = seq(0.1,7,by=0.1) # x axis
y = beta^alpha/gamma(alpha) * x^(alpha-1) * exp(-beta*x)
y = y/sum(y) # normalized
Ts_gamma = list(x=x,y=y)
plot(Ts_gamma)
abline(v=2.5)
sum(Ts_gamma$y[Ts_gamma$x>=1.5&Ts_gamma$x<=4.0])

yscale = max(Ts_gamma$y)*1.6
sum(Ts_gamma$y[x>=2&x<=4.5])/sum(Ts_gamma$y) 
#[1] 0.8459356

# -----------------------------------------------------------------------------
# Combined PDF: ecs x co2
# Combined VAL: dT = lambda * dF, lambda = ecs/3.7, dF = 5.35log(CO2b/CO2a)
# -----------------------------------------------------------------------------

# probabilty distro 
# sum(Ts$y[1:21])      # .0637
# sum(Ts$y[22:46])     # .8405
# 1 - sum(Ts$y[1:46])  # .0958

# sample co2 concentrations to get vector same length as Ts_gamma
n = length(Ts_gamma$y)
m = length(co2_ppm_all)
i = seq(1,m,as.integer(m/(n-1)))

# these PDFs were calculated above
CO2_PDF = co2_ppm_ff_all_pdf[i]  # gamma fit to Meinshausen
CO2_PDF=CO2_PDF/sum(CO2_PDF) # normalize
DF_PDF=CO2_PDF
LAMBDA_PDF = Ts_gamma$y			# gaussian

# these values were calculated above
CO2_VAL = co2_ppm_all[i]
DF_VAL = 5.35 * log(CO2_VAL/(co2_ppm_preindustrial + co2_ppm_other))
LAMBDA_VAL = Ts_gamma$x / 3.7 

# calculate 2D probability matrix
DT_PDF = matrix(LAMBDA_PDF,n,n) %*% diag(DF_PDF)
DT_PDF = DT_PDF/sum(DT_PDF) # normalized
#image(DT_PDF)

# calculate 2D temperature response
DT_VAL = matrix(LAMBDA_VAL,n,n) %*% (diag(DF_VAL))
#image(DT_VAL)

# calculate probability that CO2 > 1000 ppm or warming > 6C
b = matrix(CO2_VAL>1000,n,n) + matrix(DT_VAL>6,n,n) > 0
sum(DT_PDF[b])
# [1] 0.1003253 without methane
# [1] 0.3168074 with methane

# calculate probability that warming > 4C
b = matrix(DT_VAL>4,n,n)
sum(DT_PDF[b])
# [1] 0.4979181 without methane
# [1] 0.7712913 with methane 

# pdf for resulting warming
s = as.integer(max(DT_VAL))*10+1
dt_pdf = rep(NA,s-1)
for (i in 1:s) {
	b = matrix(DT_VAL>=(i/10),n,n) & matrix(DT_VAL<((i+1)/10),n,n) > 0
	dt_pdf[i+1]=sum(DT_PDF[b])
}
dt_pdf[is.na(dt_pdf)]=0
#dt_pdf

# this is a hack
# run this whole script twice
# just to this point to initialize these two
# once with methane = 1000 gt (65 ppm) and once with methane = 0
if (co2_ppm_methane > 0 ) {
	dt_pdf_methane = dt_pdf
} else {
	dt_pdf_no_methane = dt_pdf
}



# -----------------------------------------------------------------------------
# Plot Gaussian for Future Fossil Fuel Gigatonnes Availability in Gigatonnes CO2
# -----------------------------------------------------------------------------
mu = co2_gt_ff_future_mu
sig = co2_gt_ff_future_sig
y = co2_gt_ff_future_pdf
x = co2_gt_ff_future
xtck = as.integer(seq(mu-4*sig,mu+4*sig,sig/2)) # x-axis ticks
xoffset = 250 # for placing text

if (graphic_output) png("img/boe-ff-emit-gaussian.png",w=600,h=400)
	par(oma=c(0,0,0,0))
	plot(y~x,type="l",ylab="",yaxt="n",xlab="gigatonnes of CO2",xaxt="n")
	axis(1, at=xtck,labels=xtck,las=2)
	abline(v=mu+sig, col="blue",lty=3,lwd=1.5)  # +1 sig
	abline(v=mu-sig, col="blue",lty=3,lwd=1.5)  # -1 sig
	abline(v=mu+2*sig, col="red",lty=3,lwd=1.5) # +2 sig
	abline(v=mu-2*sig, col="red",lty=3,lwd=1.5) # -2 sig
	abline(v=mu, col="green",lty=3,lwd=1.5) # mean
	# this is the proven reserves
	points(x[as.integer(1+mu-2*sig)],y[as.integer(1+mu-2*sig)],col="red",cex=1.5)
	text(x[as.integer(1+mu-2*sig)]-xoffset,y[as.integer(1+mu-2*sig)]*1.2,"a")
	# this is the mean derived
	points(x[1+as.integer(mu)],y[1+as.integer(mu)],col="green",cex=1.5)
	text(x[1+as.integer(mu)]-xoffset,y[1+as.integer(mu)]*1.015,"c")
	# this is the 'most optimistic' from OGJ (unconfirmed)
	points(x[1+as.integer(mu+sig)],y[1+as.integer(mu+sig)],col="blue",cex=1.5)
	text(x[1+as.integer(mu+sig)]-xoffset,y[1+as.integer(mu+sig)]*1.015,"b")
	text(x[as.integer(2*mu*0.8)],y[as.integer(mu*0.8)],paste("sigma = \n",as.integer(sig)," gigatonnes CO2",sep=""),cex=0.8)
	title(main="Future CO2 Emissions From Fossil Fuels")
	legend("topleft", legend = c("Scenarios:", "(a) Proven", "(b) Optimistic","(c) Most Likely"),
			bty="n", xjust = 1, yjust = 1, cex=0.9)
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")
if (graphic_output) dev.off()

# -----------------------------------------------------------------------------
# Plot Gaussian for Resulting Atmospheric CO2 Concentration from Fossil Fuel
# (past and future)
# -----------------------------------------------------------------------------
mu = co2_ppm_ff_all_mu
sig = co2_ppm_ff_all_sig
y = co2_ppm_ff_all_pdf
x = co2_ppm_ff_all
xtck = as.integer(seq(mu-4*sig,mu+4*sig,sig/2)) # x-axis ticks
xoffset = 20 # for placing text
xscale=length(x)/(x[length(x)]-x[1])

if (graphic_output) png("img/boe-co2-ppm-gaussian.png",w=600,h=400)
	plot(y~x,type="l",ylab="",yaxt="n",xlab="PPM of CO2",xaxt="n")
	axis(1, at=xtck,labels=xtck,las=2)
	abline(v=mu+sig, col="blue",lty=3,lwd=1.5)  # +1 sig
	abline(v=mu-sig, col="blue",lty=3,lwd=1.5)  # -1 sig
	abline(v=mu+2*sig, col="red",lty=3,lwd=1.5) # +2 sig
	abline(v=mu-2*sig, col="red",lty=3,lwd=1.5) # -2 sig
	abline(v=mu, col="green",lty=3,lwd=1.5) # mean
	# this is the proven reserves
	i = round((mu-2*sig-x[1])*xscale) + 1
	points(x[i],y[i],col="red",cex=1.5)
	text(x[i]-xoffset,y[i]*1.2,"a")
	# this is the mean derived
	i = round((mu-x[1])*xscale) + 1
	points(x[i],y[i],col="green",cex=1.5)
	text(x[i]-xoffset,y[i]*1.015,"c")
	# this is the 'most optimistic' from OGJ (unconfirmed)
	i = round((mu+sig-x[1])*xscale) + 1
	points(x[i],y[i],col="blue",cex=1.5)
	text(x[i]-xoffset,y[i]*1.015,"b")
	i = round((mu+3*sig-x[1])*xscale)
	text(x[i],y[round(length(x)/2)*0.8],paste("sigma = \n",round(sig)," PPM CO2",sep=""),cex=0.8)
	title(main="CO2 Concentration due to Fossil Fuel Sources\n(Past and Future, Airborne Fraction)")
	legend("topleft", legend = c("Scenarios:", "(a) Proven", "(b) Optimistic","(c) Most Likely"),
			bty="n", xjust = 1, yjust = 1, cex=0.9)
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")
if (graphic_output) dev.off()


# -----------------------------------------------------------------------------
# Plot Temperature Response for PDF of CO2 and 3 key ECS (2,3,4.5)
# (CO2 from all sources)
# -----------------------------------------------------------------------------
mu = co2_ppm_ff_all_mu + co2_ppm_preindustrial + co2_ppm_other + co2_ppm_methane
sig = co2_ppm_ff_all_sig
y = rep(0,length(co2_ppm_all))
x = co2_ppm_all
xtck = as.integer(seq(mu-4*sig,mu+4*sig,sig/2)) # x-axis ticks
xoffset = 20 # for placing text
xscale=length(x)/(x[length(x)]-x[1])
 
if (graphic_output) png("img/boe-co2-2sig-ecs-1sig-gaussian.png",w=600,h=400)
	par(oma=c(0,0,0,0))
	plot(y~x, type="l",ylim=c(1,8),ylab="dT (deg C)",xaxt="n",
					xlab="CO2 PPM Total (all sources, future and past)", lty=1)
	axis(1, at=xtck,labels=xtck,las=2)
	lines(x, dT45, lty=1) 
	lines(x, dT3, lty=2) 
	lines(x, dT2, lty=3)
	abline(h=6,col="red",lwd=2) # danger zone
	abline(v=1000,col="red",lwd=2) # danger zone
	abline(v=mu-2*sig,col="red",lty=3) # -2 sig
	abline(v=mu+2*sig,col="red",lty=3) # +2 sig
	abline(v=mu-sig,col="blue",lty=3) # -1 sig
	abline(v=mu+sig,col="blue",lty=3) # +1 sig
	abline(v=mu,col="green",lty=3) # max
	text(mu+3*sig,2.2,paste("sigma = \n",round(sig)," PPM CO2"),cex=0.8)
	title(main="Temperature Response for \nThree Equlibrium Climate Sensitivities (1-sigma)\nand Probablistic CO2 Concentration (1-,2-sigma)", 
		cex.sub=0.7,cex.main=1)
	legend("topleft", legend = c("Ts = 4.5C", "Ts = 3.0C","Ts = 2.0C"),
			bty="n", xjust = 1, yjust = 1, cex=0.8, lty=c(1,2,3))
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")
if (graphic_output) dev.off()

# -----------------------------------------------------------------------------
# Plot Poisson Distribution of ECS
# (with and without overlay on Meinshausen 2009
# -----------------------------------------------------------------------------
# the size of the png is set to work with meinshausen overlay
xxtck = 0:7 # labels for the chart
if (graphic_output) png("img/boe-ecs-pdf-poisson.png",w=512,h=295)
	plot(Ts_poisson, type="l",yaxt="n",ylab="",xaxt="n",xlab="T (deg C)",
		xaxs="i",yaxs="i",lwd=4,col="red",ylim=c(0,max(Ts_poisson$y)*1.6))
	axis(1, at=xxtck,labels=xxtck,)
	abline(v=2,col="blue",lty=3,lwd=2)
	abline(v=4.5,col="blue",lty=3,lwd=2)
	abline(v=3,col="green",lty=3,lwd=2)
	title(main="Poisson PDF for Equilbrium Climate Sensitivity")
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")
if (graphic_output) dev.off()

library("rimage")
meinshausenECSpdf <- read.jpeg("img/meinshausen-fig3-cropped.jpg")
if (graphic_output) png("img/boe-ecs-poisson-meinshausen.png",w=512,h=295)
	plot(meinshausenECSpdf)
	par(new=T)
	plot(Ts_poisson, type="l",yaxt="n",ylab="",xaxt="n",xlab="T (deg C)",
		xaxs="i",yaxs="i",lwd=4,col="red",ylim=c(0,max(Ts_poisson$y)*1.6))
	axis(1, at=xxtck,labels=xxtck,)
	abline(v=2,col="blue",lty=3,lwd=2)
	abline(v=4.5,col="blue",lty=3,lwd=2)
	abline(v=3,col="green",lty=3,lwd=2)
	title(main="Poisson PDF for Equilbrium Climate Sensitivity\nOverlay on Meinshausen 2009 Fig 3a")
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")
if (graphic_output) dev.off()


# -----------------------------------------------------------------------------
# Plot Gamma Distribution of ECS
# (with and without overlay on Meinshausen 2009
# -----------------------------------------------------------------------------
# the size of the png is set to work with meinshausen overlay
xxtck = 0:7 # labels for the chart
if (graphic_output) png("img/boe-ecs-pdf-gamma.png",w=512,h=295)
	plot(Ts_gamma, type="l",yaxt="n",ylab="",xaxt="n",xlab="T (deg C)",
		xaxs="i",yaxs="i",lwd=4,col="red",ylim=c(0,yscale))
	axis(1, at=xxtck,labels=xxtck,)
	abline(v=2,col="blue",lty=3,lwd=2)
	abline(v=4.5,col="blue",lty=3,lwd=2)
	abline(v=3,col="green",lty=3,lwd=2)
	title(main="Gamma PDF for Equilbrium Climate Sensitivity")
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")
if (graphic_output) dev.off()

meinshausenECSpdf <- read.jpeg("img/meinshausen-fig3-cropped.jpg")
if (graphic_output) png("img/boe-ecs-gamma-meinshausen.png",w=512,h=295)
	plot(meinshausenECSpdf); par(new=T)
	plot(Ts_gamma,type="l",yaxt="n",ylab="",xlab="T (deg C)",
		xaxs="i",yaxs="i",lwd=4,col="red",ylim=c(0,yscale),xlim=c(0,7))
	axis(1, at=xxtck,labels=xxtck,)
	abline(v=2,col="blue",lty=3,lwd=2)
	abline(v=4.5,col="blue",lty=3,lwd=2)
	abline(v=3,col="green",lty=3,lwd=2)
	title(main="Gamma PDF for Equilbrium Climate Sensitivity\nOverlain on Meinshausen 2009 Fig 3a")
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")
if (graphic_output) dev.off()

# -----------------------------------------------------------------------------
# Plot Poisson -v- Gamma
# -----------------------------------------------------------------------------
if (graphic_output) png("img/boe-ecs-gamma-poisson.png",w=512,h=295)
	plot(Ts_poisson, type="l",yaxt="n",ylab="",xaxt="n",xlab="T (deg C)",
		xaxs="i",yaxs="i",lwd=4,col="red",ylim=c(0,yscale*.8))
	lines(Ts_gamma,lwd=4,col="blue")
	axis(1, at=xxtck,labels=xxtck,)
	abline(v=2,col="blue",lty=3,lwd=2)
	abline(v=4.5,col="blue",lty=3,lwd=2)
	abline(v=3,col="green",lty=3,lwd=2)
	title(main="Poisson -v- Gamma ECS estimates")
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")
if (graphic_output) dev.off()


# -----------------------------------------------------------------------------
# Plot Combined Graphs for ECS/CO2 combined PDFs and dT
# -----------------------------------------------------------------------------
# contours for combined pdf
t=list(x=Ts_gamma$x,y=CO2_VAL,z=DT_VAL)
a=list(x=Ts_gamma$x,y=CO2_VAL,z=DT_PDF)
az = sort(a$z,decreasing=T)
azcumsum = az
for (i in 1:length(az)) {
	azcumsum[i] = sum(az[1:i])
}
azsig2 = az[sum(azcumsum<=.95)]
azsig1 = az[sum(azcumsum<=.84)]
az50 = az[sum(azcumsum<=.50)]


#================================
if (graphic_output) png("img/boe-pdf-ecs-ppm.png",w=600,h=400)
	layout(matrix(c(0,1,7,2,0,3,5,0,6,4),2,5,byrow=TRUE), c(1,7,1,7,1), c(1,7),TRUE) 
	par(mar=c(0,0,0,0),oma = c(1,0.2,3,0.2))
	#layout.show(nf)
		plot(Ts_gamma,type="l",ylab="",yaxt="n",xlab="",xaxt="n")
			text(0.5,max(Ts_gamma$y)*0.8,"a")
			text(5.5,max(Ts_gamma$y)*0.8,"ECS")
			text(5.5,max(Ts_gamma$y)*0.4,"PDF")
		plot(Ts_gamma,type="l",ylab="",yaxt="n",xlab="",xaxt="n")
			text(0.5,max(Ts_gamma$y)*0.8,"a")
			text(5.5,max(Ts_gamma$y)*0.8,"ECS")
			text(5.5,max(Ts_gamma$y)*0.4,"PDF")
		plot(-CO2_PDF,CO2_VAL,type="l",ylab="",yaxt="n",xlab="",xaxt="n")
			text(-max(CO2_PDF)*0.8,CO2_VAL[60],"b")
			text(-max(CO2_PDF)*0.5,CO2_VAL[12],"CO2")
			text(-max(CO2_PDF)*0.5,CO2_VAL[7],"PPM")
			text(-max(CO2_PDF)*0.5,CO2_VAL[2],"PDF")
		plot(CO2_PDF,CO2_VAL,type="l",ylab="",yaxt="n",xlab="",xaxt="n")
			text(max(CO2_PDF)*0.8,CO2_VAL[60],"b")
			text(max(CO2_PDF)*0.5,CO2_VAL[12],"CO2")
			text(max(CO2_PDF)*0.5,CO2_VAL[7],"PPM")
			text(max(CO2_PDF)*0.5,CO2_VAL[2],"PDF")
	image(a$z,xaxt="n",yaxt="n")
		axis(1, at=a$x[seq(1,70,10)]/a$x[70],labels=as.integer(a$x[seq(1,70,10)]))
		axis(3, at=a$x[seq(1,70,10)]/a$x[70],labels=rep(" ",7))
		axis(4, at=seq(0,1,by=.1),labels=as.integer(a$y[seq(1,70,6.9)]),las=2,cex=0.7)
	b.full=contourLines(a$z,levels=c(az50,azsig1,azsig2))
	line.type=c(1,2,3)
	line.width=c(0.8,0.8,0.8)
	for(c in 1:length(b.full)){
		  ##This loop is used in case the density is multimodal or if the desired contour
		  ##  extends outside the plotting region
		  b=list("x"=as.vector(unlist(b.full[[c]][2])),"y"=as.vector(unlist(b.full[[c]][3])))
		  
		  ##plot desired contour
		  line.dat<-xspline(b,open=TRUE,draw=FALSE)
		  lines(line.dat,lty=line.type[c],lwd=line.width[c])
	}
	 par(new=T)
	 plot.new()
	 text(0,1,"c",col="black")

	 
	#=================================
	# contour lines on integer dT
	ncon = as.integer(max(t$z))-as.integer(min(t$z))
	d.full=contourLines(t$z,levels=seq(0,ncon))

	# chart
	#def.par <- par(no.readonly = TRUE)
	image(t$z,xaxt="n",yaxt="n",col = heat.colors(ncon))
		axis(1, at=a$x[seq(10,70,10)]/a$x[70],labels=as.integer(a$x[seq(10,70,10)]))
		axis(3, at=a$x[seq(10,70,10)]/a$x[70],labels=rep(" ",7))
		axis(2, at=seq(0,1,by=.1),labels=rep(" ",11))
	for(c in 1:length(d.full)){
		  ##This loop is used in case the density is multimodal or if the desired contour
		  ##  extends outside the plotting region
		  d=list("x"=as.vector(unlist(d.full[[c]][2])),"y"=as.vector(unlist(d.full[[c]][3])))
		  
		  ##plot desired contour
		  line.dat<-xspline(d,open=TRUE,draw=FALSE)
		  lines(line.dat,lty=2,col="green")
	}
	for(c in 1:length(b.full)){
		  ##This loop is used in case the density is multimodal or if the desired contour
		  ##  extends outside the plotting region
		  b=list("x"=as.vector(unlist(b.full[[c]][2])),"y"=as.vector(unlist(b.full[[c]][3])))
		  
		  ##plot desired contour
		  line.dat<-xspline(b,open=TRUE,draw=FALSE)
		  lines(line.dat,lty=line.type[c],lwd=line.width[c])
	}
	 par(new=T)
	 plot.new()
	 text(0,1,"d",col="black")
	 # title
	 mtext("Combined PDF and Temperature Response\nFor Probablistic Distributions of ECS and Fossil Fuel Emissions", side=3,outer = TRUE)
	 mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), side=1, adj=1, cex=0.7, col="dark gray",outer= TRUE)
	 plot.new()
	 text(0.5,0.7,"CO2",col="black")
	 text(0.5,0.3,"PPM",col="black")
if (graphic_output) dev.off()
par(def.par)#- reset to default

#===============================================================
#
#===============================================================
dt_pdf_methane = dt_pdf_methane[2:length(dt_pdf_methane)]
dt_pdf_no_methane = dt_pdf_no_methane[2:length(dt_pdf_no_methane)]
xx= seq(0.2, (length(dt_pdf_methane)+1)/10 , by=0.1)
x= seq(2, length(dt_pdf_methane)+1)

xx[dt_pdf_methane==max(dt_pdf_methane)]
# 4.8
xx[dt_pdf_no_methane==max(dt_pdf_no_methane)]
# 3.9

alpha = 11.2;beta = 0.21
y_ch4 = (beta)^alpha/gamma(alpha) * (x)^(alpha-1) * exp(-(beta)*(x))
plot(dt_pdf_methane,xaxt="n");lines(y_ch4)
sum((y_ch4-dt_pdf_methane)^2)
mu_ch4=alpha/beta
#[1] 53.33333
t_ch4 = x[y_ch4==max(y_ch4)] -1
# 48

xxx = x[1:length(dt_pdf_no_methane)]
alpha = 10.1;beta = 0.24
y_no = (beta)^alpha/gamma(alpha) * (xxx)^(alpha-1) * exp(-(beta)*(xxx))
plot(dt_pdf_no_methane~xxx);lines(y_no)
sum((y_no-dt_pdf_no_methane)^2)
mu_no=alpha/beta
#[1] 42.08333
t_no = xxx[y_no==max(y_no)] -1
# 37

if (graphic_output) png("img/boe-pdf-for-two-with-methane.png",w=600,h=400)
if (is.vector(dt_pdf_methane) & is.vector(dt_pdf_no_methane)) {
	plot(dt_pdf_no_methane~xxx, col="blue", type="l",xaxt="n",ylab="",xlab="warming (degC)",xaxs="i",yaxs="i")
	axis(1, at=seq(10,130,by=10),labels=seq(10,130,by=10)/10)
	lines(dt_pdf_methane, col="red")
	lines(y_ch4,col="red",lwd=2)
	lines(y_no,col="blue",lwd=2)
	abline(v=t_ch4, col="red",lty=2)
	abline(v=t_no, col="blue",lty=2)
	title(main="Warming Scenarios\nWith and Without Methane Feedback")
	legend("topright", legend = c("No Feedbacks","Methane Feedback"), col=c("blue","red"),
			bty="n", xjust = 1, yjust = 1, cex=0.9, lty=1)
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")

}
if (graphic_output) dev.off()