#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) }