#Data recorded from a full french press at five minute intervals.  Temperatures are in
#Farenheit.
french<- structure(list(Time = c(0L, 5L, 10L, 15L, 20L, 25L, 30L, 35L, 
40L, 45L, 50L, 55L, 60L, 65L, 70L, 75L, 80L, 85L, 90L, 95L), 
    Temp = c(201.4, 195.7, 192.1, 188.9, 185.6, 181.8, 180.1, 
    177.3, 174.8, 172.4, 169.8, 167.6, 165.2, 163.2, 161, 159, 
    156.8, 154.7, 152.9, 151.1)), .Names = c("Time", "Temp"), class = "data.frame", row.names = c(NA, -20L)) 
#Add another column for Kelvin
french$kelvin<- (french$Temp + 459.67)*5/9


french.lin<-lm(Temp ~ Time, data=french)
french.quad<-lm(Temp ~ Time + I(Time^2) , data=french)

lin.model.plot<- function() {
	with(french, plot(Time, Temp ,ylim=c(140,210),xlim=c(0,100),main="Temperature Drop in a Frieling Stainless Steel French Press:\n Linear and Quadratic models", xlab="Minutes After Brewing", ylab="Temp (F)", pch=19, cex=0.6))
	abline(a=coef(french.lin)[1],b=coef(french.lin)[2], col=2, lty=2)
	lines(french$Time,predict(french.quad), col=3)
	legend(x=75,y=185, legend=c("Linear","Quadratic"),fill=c(2,3),cex=0.8)
	}

#Some useful constants for estimating heat transfer.  Weight and Area taken from Frieling's 
#Amazon product page for this model http://www.amazon.com/gp/product/B00009ADDS/
c.water<- 4181.3
kg.water<- 0.975926476 
a.water<- 5.08*22.86*pi*2*0.0001
ambient<- (67 + 459.67)*5/9
C.press<- c.water*kg.water
#Initial temp for nls model
T.O<- french$kelvin[1]

#Generic exponential function for heat transfer by convection
exp.temp<- function(time.press, h.est) {
	exp(-h.est*(a.water/C.press)*time.press)
	}

#fit a nonlinear model.  I used 1 as a starting value so nls() didn't complain.  turns out it
#wasn't a bad guess
french.nls<- nls(kelvin ~ ambient + (T.O - ambient)*exp.temp(time.press=french$Time*300, beta1), data=french, start=list(beta1=1))
	
nls.model.plot<- function(beta1) {
	with(french, plot(Time*300, kelvin, xlim=c(0,30000), ylim=c(335,375), main="Exponential Decay Model", xlab="Seconds After Brewing", ylab="Temp (K)", pch=19, cex=0.6))
	pred.temp<- with(french, ambient + (T.O - ambient)*exp.temp(time.press=french$Time*300, coef(french.nls)))
	lines(french$Time*300, pred.temp, lty=3, col=4) 
	}

#I guess there are less ugly ways to drop the first observation.  Had I thought about this from the start 
#then data=french[-1,] would have sufficed but I hard coded a few constants and had to walk back through everything.
T.1<- french$kelvin[2]
french2<- french[-1,]
time2<- french2$Time*300

french.nls2<- nls(kelvin ~ ambient + (T.1 - ambient)*exp(-beta2*(a.water/C.press)*time2), data=french2, start=list(beta2=1))

nls.model.plot2<- function(beta1, beta2) {
	with(french, plot(Time*300, kelvin, xlim=c(0,30000), ylim=c(335,375), main="Exponential Decay Model Comparison", xlab="Seconds After Brewing", ylab="Temp (K)", pch=19, cex=0.6))
	pred.temp<- with(french, ambient + (T.O - ambient)*exp.temp(time.press=french$Time*300, coef(french.nls)))
	pred.temp2<- with(french2, ambient + (T.1 - ambient)*exp(-beta2*(a.water/C.press)*time2))
	lines(french$Time*300, pred.temp, col=4) 
	lines(time2, pred.temp2, col=2)
	legend(x=22000,y=360, legend="New model",fill=2)
	}