#SCRIPT JUNE 29
#requires R-package png
 # install.packages("png")
#IDEAL GAS LAW RELATIONS  pV=nRT
pascal= 6894.76 #pascals per 1 psi
R= 8.31441 #gas constant J K-1 mol-1 
Vf= .004237 #m3 volume of football 4237 cm3
cent =function(x) (x-32)*5/9
fahr = function(x) 9/5*x +32
coltp=13.0
bias=.38 #between Logo and Nonlogo Gauge

Nx=function(targetp,pregame=71) (6894.76*(14.7+targetp))*Vf/(R*(273+cent(pregame)))
 #function to calculate moles from target pressure in psig and pregame temperature in deg 
 Nc= Nx(coltp,pregame=71);Nc # 0.3302896

Pf=  function(tem, n=Nc) (n*R*(273+cent(tem))/Vf)/6894.76 - 14.7   #in psi 
  #function to calculate pressure in psig from temperature in deg F and moles 
  Pf(48,n=Nx(coltp,pregame=71)) #11.798

Tf= function(P, n=Nc) fahr( (Vf* ((P+14.7)*6894.76))/(n*R)  -273 )  #in psi 
  #function to calcuate temperature in deg F from pressure in psig and moles
  Tf(11.798,n=Nc)#47.98409  
 
#CONVERSION FORMULAS TO MASTER GAUGE
###################################
logof=  function(x)  (x+.2846)/1.05   # Logo to Master: Master = Logo + 0.2836 psi/1.05
  logof(c(11.49,12.74))#11.214 12.404
nonlogof= function(x) (x+.1444)/1.015 #and Master = Non-Logo + 0.1444 ps /1.015
  nonlogof(c(11.11,12.33)) #11.088 12.290
logoi=function(x) 1.05*x -.2846 #Logo from Master
nonlogoi=function(x) 1.015*x -.1444 #Nonlogo from master

##INPUT FIT TO DIGITIZED TRANSIENTS IN FIGURES 25 (Non-Logo) and 27 (Logo) 
B=read.csv("http://www.climateaudit.info/data/football/transient_coef.csv")
B$pregame=rep(c(71,67),each=4)
B$gauge=factor(B$gauge,levels=c("Logo","Non-Logo","Master"))

setwd("d:/2015/football")
foot=read.csv("http://www.climateaudit.info/data/football/football.csv")
sim=read.csv("http://www.climateaudit.info/data/football/sim_table14.csv",nrow=6 )

#FIG27 WITH TRANSIENT CALCULATIONS
dest="d:/temp/temp.dat"
o=F #controls output file
ind=seq(0,14,.2)
library(png)
download.file("http://www.climateaudit.info/data/football/football_fig27.PNG",dest,mode="wb")
img27=readPNG(dest) ##"football_fig27.png"

#SHOW FIG27 ANNOTATED FOR BLOG POST
if(o) png("figure27_annotated.png",w=640,h=560)
par(mar=c(0,5,2,5))
if(o) par(list(font.lab=2,font.axis=2,cex.axis=1.5,cex.lab=1.5))
plot(0,type="n",xlim=c(-.25,14),xaxs="i",ylim=c(10.25,12.8),axes=F,ylab="")
rasterImage(img27,-3.22,10.18,17.51,12.92)
title("Re-stated Figure 27")
mtext(font=2,line=-.45,"67 Deg F Logo Gauge Initialization",cex=.9)
axis(side=2,at=seq(10.6,12.8,.2),las=1)
mtext(side=2, "PSIG (Master Gauge)",line=3.5,font=2)
axis(side=4,at=logof(seq(11,13,.2)), labels=seq(11,13,.2), las=1)
mtext(side=4, "PSIG (Conversion to Logo Scale)",line=3.5,font=2)

#Colt
k=7
f=function(x) B$A[k]+B$B[k]*exp(-B$C[k]*x)
lines(ind,f(ind),col=4,lwd=1)
f(0) # 11.90
g=function(x)  (Pf (48,n=Nx(x,pregame=67))- f(0))^2    #46.12
P0=optimize(g, c( 12,14))$minimum;P0  #12.898
logoi(P0) # 13.25807
points(0, Pf (48,n=Nx(logof(13),pregame=67)),pch="+",col=4,cex=2)
delta=f(0)- Pf (48,n=Nx(logof(13),pregame=67));delta # 0.2369
f=function(x) B$A[k]+B$B[k]*exp(-B$C[k]*x)-delta
lines(ind,f(ind),col=4,lwd=4)
lines(ind,f(ind)-.45,col=4,lwd=2,lty=2)
points(8,logof(mean(foot$Logo[12:15])),pch=19,cex=2,col=4)
text(3.59,11.94,col=4,pos=4,font=2,"13.0 psi")
 text(3.48,12.2,col="steelblue",pos=2,font=2,"13.26 psi")

#Patriot
k=5
f=function(x) B$A[k]+B$B[k]*exp(-B$C[k]*x)
lines(ind,f(ind),col=2,lwd=1)
f(0) # 11.49
g=function(x)  (Pf (48,n=Nx(x,pregame=67))- f(0))^2    #46.12
P1=optimize(g, c( 12,14))$minimum ;P1 #12.4702
P1+bias #  12.8502
logoi(P1) # 12.80911
Pf (48,n=Nx(logof(12.5),pregame=67)) # 11.20
delta=f(0)- Pf (48,n=Nx(logof(12.5),pregame=67)) ;delta# 0.2837
points(0, Pf (48,n=Nx(logof(12.5),pregame=67)),pch="+",col=2,cex=2)
f=function(x) B$A[k]+B$B[k]*exp(-B$C[k]*x)-delta
lines(ind,f(ind),col=2,lwd=4)
lines(ind,f(ind)-.45,col=2,lwd=2,lty=2)
points(3.75,logof(mean(foot$Logo[1:11])),pch=19,cex=2,col=2)
text(3.2,11.49,col=2,pos=4,font=2,"12.5 psi")
text(2.7,11.73,col="red4",pos=4,font=2,"12.81 psi")

if(o) dev.off()

##FIGURE 25 ANNOTATION
download.file("http://www.climateaudit.info/data/football/football_fig25.PNG",dest,mode="wb")
img25=readPNG(dest) ##"football_fig25.png"

if(o) png("figure25_annotated.png",w=640,h=560)
par(mar=c(0,5,2,5))
if(o) par(list(font.axis=2,font.lab=2,cex.axis=1.5,cex.lab=1.5))
plot(0,type="n",xlim=c(-0.12,14),xaxs="i",ylim=c(10.25,12.9),axes=F,ylab="")
rasterImage(img25,-3.23,10.225,17.61,12.92)
title("Figure 25 Annotated")
mtext(side=3,font=2,"71 deg F Non-Logo Initialization",line=-1)
axis(side=2,at=seq(10.6,12.8,.2),las=1)
mtext(side=2, "PSIG (Master Gauge)",line=3.5,font=2)
#axis(side=2,at=seq(10.6,12.8,.2)-del25,labels= seq(10.6,12.8,.2),las=1)
axis(side=4,at=nonlogof(seq(10.6,12.8,.2)), labels=seq(10.6,12.8,.2), las=1)
mtext(side=4, "PSIG (Conversion to Non-Logo Scale)",line=3.5,font=2)

#Colt
k=3
f=function(x) B$A[k]+B$B[k]*exp(-B$C[k]*x)
lines(ind,f(ind),col=4,lwd=1)
f(0) # 11.85
g=function(x)  (Pf (48,n=Nx(x,pregame=71))- f(0))^2    #46.12
P0=optimize(g, c( 12,14))$minimum;P0  #13.05
nonlogoi(P0) #13.10165
P1=Pf (48,n=Nx(nonlogof(13),pregame=71)) #11.75114
delta=f(0)- P1 ;delta#0.09580259
points(0, Pf (48,n=Nx(nonlogof(13),pregame=71)),pch="+",col=4,cex=2)
f=function(x) B$A[k]+B$B[k]*exp(-B$C[k]*x)-delta
f(0) #11.63
lines(ind,f(ind),col=4,lwd=4)
lines(ind,f(ind)-.45,col=4,lwd=2,lty=2)
points(8,nonlogof(mean(foot$Nonlogo[12:15])),pch=19,cex=2,col=4)
text(3.09,11.99,col=4,pos=4,font=2,"13.0 psi")
 text(2.98,12.12,col="steelblue",pos=2,font=2,"13.1 psi")

#Patriot
k=1
f=function(x) B$A[k]+B$B[k]*exp(-B$C[k]*x)
lines(ind,f(ind),col=2,lwd=1)
f(0) # 11.42815
g=function(x)  (Pf (48,n=Nx(x,pregame=71))- f(0))^2    #46.12
P0=optimize(g, c( 12,14))$minimum;P0  #12.61252
nonlogoi(P0) #12.66
P1=Pf (48,n=Nx(nonlogof(12.5),pregame=71));P1 #11.27989
delta=f(0)- P1 ;delta#0.1482604
points(0, Pf (48,n=Nx(nonlogof(12.5),pregame=71)),pch="+",col=2,cex=2)
f=function(x) nonlogof(12.5) - (nonlogof(12.5)- Pf (48,n=Nx(nonlogof(12.5),pregame=71)))*exp(-.115*x)
f(0) #11.63
f=function(x) B$A[k]+B$B[k]*exp(-B$C[k]*x)-delta
lines(ind,f(ind),col=2,lwd=4)
lines(ind,f(ind)-.45,col=2,lwd=2,lty=2)
points(3.75,nonlogof(mean(foot$Nonlogo[1:11])),pch=19,cex=2,col=2)
text(3.09,11.53,col=2,pos=4,font=2,"12.5 psi")
 text(3,11.72,col="red4",pos=2,font=2,"12.66 psi")

if(o) dev.off()