# 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 lam2 = 2.0/3.7 lam3 = 3.0/3.7 lam45 = 4.5/3.7 # calculate equilibrium change in temp as function of CO2 forcings and ecs dT2 = lam2 * dF dT3 = lam3 * dF dT45 = lam45 * dF # ----------------------------------------------------------------------------- # ECS distribution: Poisson distribution # ----------------------------------------------------------------------------- # mean = climate sensitivity = 3C lambda = 3 # ECS mean = 3 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 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) 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()