Coffee and Cream Problem Simulation

coffeeCream_functions.R

# T.room   # Temperature of the room
# T.coffee # Initial temperature of the coffee
# T.cream  # Temperature of the cream
# k        # Cooling rate constant
# T.d      # Desired coffee temperature

temp <- function(T.init, T.room, k, x){
  # Newton's cooling law temperature function
  return(T.room + (T.init - T.room)*exp(-k*x))
}

coffeeCream <- function(params){
  # delay = delay from time zero to add cream to coffee
  # ratio = ratio of cream to coffee (ie. 100-ratio% coffee, ratio% cream)
  
  T.coffee = params[1]
  T.d      = params[2]
  T.room   = params[3]
  T.cream  = params[4]
  ratio    = params[5]
  k        = params[6]
  t.end    = params[7]
  delay    = params[8]
  
  if(delay < 0){
    stop('coffeeCream(): delay must be >= 0')
  }
  
  # Initialize plot window with blank plot
  plot(0,bty='n',pch='',ylab='Temperature',xlab='Time', 
       ylim=c(min(32, T.cream),(T.coffee)), xlim=c(0,t.end))
  
  # Draw horizontal lines for temperatures and vertical grid lines
  abline(h=T.room, lty=2)
  abline(h=T.d, lty=5, col='red')
  abline(h=T.cream, col='blue')
  abline(v=seq(0,t.end, 10), col='grey', lty=3)
  
  # Draw cooling curves
  if(delay == 0){
  
    # if delay == 0, then cream added immediately
    segments(x0=0, y0=T.coffee, x1=0, y1=(1-ratio)*T.coffee+ratio*T.cream,
          col='blue', lwd=2)
    
    # Draw cooling curve from new initial mixed temperature
    curve(temp((1-ratio)*T.coffee+ratio*T.cream, T.room, k, x), 
          from=0, to=t.end, add=TRUE)
          
  } else if(delay > 0){
    # If delay > 0, cream is added after some cooling period
    
    # Calculate temperature of coffee and then mixture after delay period
    tempAtDelay <- temp(T.coffee, T.room, k, delay)
    tempAfterCream <- (1-ratio)*tempAtDelay+ratio*T.cream
    
    # Draw cooling curve up to delay
    curve(temp(T.coffee, T.room, k, x), from=0, to=delay, add=TRUE)
    
    # Add cream
    segments(x0=delay, y0=tempAtDelay,
             x1=delay, y1=tempAfterCream, col='blue', lwd=2)
    
    # Allow to cool further
    curve(temp(tempAfterCream, T.room, k, x-delay), 
          from=delay, to=t.end, add=TRUE)
    
    # Draw reference curve of immediate mixing
    curve(temp((1-ratio)*T.coffee+ratio*T.cream, T.room, k, x), 
          from=0, to=t.end, add=TRUE, col=402, lty=3)
  }
}

server.R

library(shiny)

# source coffeeCream(params) function
source('coffeeCream_functions.R')

shinyServer(function(input,output) {
  sliderValues <- reactive({
    c(input$T.coffee, input$T.d, input$T.room, input$T.cream, 
      input$ratio, input$k, input$t.end, input$delay)
  })
  
  output$coffeeCreamPlot <- renderPlot({
    coffeeCream(sliderValues())
  })
})

ui.R

library(shiny)

# Define UI for dataset viewer application
shinyUI(pageWithSidebar(
  
  # Application title
  headerPanel("Coffee and Cream"),
  
  sidebarPanel(
    sliderInput('T.coffee', 'Initial Temp of Coffee:',
                min = 75, max = 175, step = 5, value = 175),
    
    sliderInput('T.d', 'Desired Temp of Coffee:',
                min = 50, max = 125, step = 5, value = 100),
    
    sliderInput('T.room', 'Room Temperature:',
                min = 60, max = 80, step = 1, value = 68),
    
    sliderInput('T.cream', 'Temp of Cream:',
                min = 32, max = 60, step = 2, value = 34),
    
    sliderInput('ratio', '% of cream in mixed liquid (5% = 0.05):',
                min = 0, max = 0.25, step = 0.05, value = 0.15),
    
    sliderInput('k', 'Rate of cooling coefficient:',
                min = 0.01, max = 0.1, step = 0.005, value = 0.08),
    
    sliderInput('t.end', 'End time:',
                min = 20, max = 80, step=10, value = 50),
    
    sliderInput('delay', 'Delay before adding coffee:',
               min=0, max=50, step=1, value=0,
               animate=animationOptions(interval=500, loop=T)),
    
    helpText('When should I add cream to my coffee?')
    ),
  
  mainPanel(
    tabsetPanel(
      tabPanel('Plot',  plotOutput('coffeeCreamPlot'))
      )
    )
))

To run this app, first make sure that you have R installed. I highly recommend you also install the excellent R Studio.

Then, open R Studio to start R. Or, at the command line (using Terminal.app on Mac or by running cmd on Windows), type R and hit enter to start R.

Once in R, make sure you have Shiny installed by running

install.packages("shiny")
library(shiny)

After loading shiny, to start up a local web sever so that you can view this shiny app, just run

runGist(9370271)

If you happen to get a download has nonzero exit status when you run the runGist command, then run the following and try again.

options('download.file.extra' = '--no-check-certificate')

As mentioned here.