Coffee and Cream (Exercise 7)

function coffee_model()
% Modeling of coffee cooling in cermaic mug. Cream is added at different
% times during the cooling process.
% Assumption: When cream is added, completely assimilates into the coffee
% and becomes one mass, with the same specific heat of the coffee.
    %% Coffee mug specifications
    radius = 0.04; %m
    height = 0.10; %m
    V = pi * radius^2 * height; % cylindrical cup of coffee
    k = 1.5; % W/m/K
    A = height * 2 * pi * radius + pi * radius ^ 2; % sides and bottom
    d = 0.007; % m
    m_coffee = V * 1000; % 1000kg/m^3 density of coffee
    m_mix = m_coffee;
    m_cream = 1e-4 * 1000; % 100ml = 1e-4m^3 cream, 1000kg/m^3 density
    c = 4186; % J/kg/K
    h = 100; % W / m^2 / K

    %% Initial Temperatures
    Tinit = 370; % K
    Troom = 290; % K
    Tcream = 275; % K
    Tdrinkable = 320; % K

    %% Times
    time_init = 0;
    time_final = 60 * 60;

   %% Main function call
   figure(4)
   plot_time_drinkable()
   get_time_to_add_cream()
    %% Describes instantaneous change of the internal energy of the coffee
    function res = coffee_U_change(t, U)
        % Change in coffee temperature is caused by conduction of heat
        % through the ceramic mug
        temp_coffee = U / m_mix / c;
        conduction = k * A / d * (Troom - temp_coffee);
        res = conduction;
    end
    %% Returns the internal energy of the coffee cream mixture after
    % cream is added
    function final_energy = energy_coffee_cream(U)
        % Changes the mass of the coffee to the mass of coffee cream
        % mixture
        final_temp = (m_coffee * (U / m_coffee / c) + m_cream * Tcream)...
            / (m_coffee + m_cream);
        m_mix = m_coffee + m_cream;
        final_energy = final_temp * m_mix * c;
    end

    %% Function that computes the time it takes for the coffee to reach drinkable
    % temperatures
    function res = time_to_drinkable(t_cream)
        % Set a trigger event for the ode function to stop at
        options = odeset('Events', @event_function);
        m_mix = m_coffee;
        % TE: time at which event happens
        % VE: value of output at which event happens
        [time_1, uu_1, TE_1, VE_1] = ode45(@coffee_U_change, [time_init, t_cream], ...
           Tinit * m_coffee * c, options);
        [time_2, uu_2, TE_2, VE_2] = ode45(@coffee_U_change, [t_cream, time_final], ...
          energy_coffee_cream(uu_1(end)), options);
        % Check which trigger actually stopped because it reached drinkable
        % temperature
        U_drinkable_before_cream = Tdrinkable * c * m_coffee;
        U_drinkable_after_cream = Tdrinkable * c * m_mix;
        % If the coffee has reached drinkable temperature during the
        % addition of the cream, the time it reaches drinkable temp will be
        % right after cream is added
        res = time_2(1);
        % Check to see if the trigger event is reasonably close to the
        % drinkable internal energy of coffee
        if abs(VE_1 - U_drinkable_before_cream) < 0.0001
            % Hits the drinkable temperature before cream is added
            res = TE_1;
        end
        if abs(VE_2 - U_drinkable_after_cream) < 0.0001
            % Hits the drinkable temperature after the cream is added
            res = TE_2;
        end

        %% Trigger stop when the temperature of the coffee is at the drinkable temp
        function [value, isterminal, direction] = event_function(t, U)
           value = U / m_mix / c - Tdrinkable; % When value = 0, event is triggered
           isterminal = 1; % stop function as soon as this event is reached
           direction = 0; % At any zero (when value = 0) consider this event, not
                          % only when the function is either increasing (1) or
                          % decreasing (-1)
        end
    end

    %% Plot time it takes for the coffee to become drinkable vs time at which
    % the cream was added
    % At some point, the cream is added after the coffee gets to drinking
    % temp by itself.
    function plot_time_drinkable()
        % If you add the cream after 330 s, the coffee is already at or
        % below the drinkable temperature
        t_cream = time_init + 1 : 1 : time_final - 1;
        Time = zeros(1, length(t_cream));
        for i = 1 : length(t_cream)
           Time(i) = time_to_drinkable(t_cream(i));
        end
        plot(t_cream, Time, 'b*')
        xlabel('Time at which cream is added (s)')
        ylabel('Time to cool to drinkable temp (s)')
        title('Cream Adding Time vs Coffee Cooling to Drinkable Temperature Time')
    end

    %% Find time at which ading cream will generate fastest cooling process
    % Res is the time at which it is most optimal to add cream aka when the
    % cooling process is at a minimum
    function res = get_time_to_add_cream()
        [x, fval] = fminsearch(@time_to_drinkable, time_init + 1);
        res = x;
    end

%% Plot the time curve of the coffee cooling with adding cream at different
   % times during the process
   function plot_adding_cream()
        % If the time_cream is changed, make sure that T, Temp, extra_t and
        % extra_temp's num of rows is enough.
       time_cream = time_init + 1; % : 100 : time_final - 1;
       T = zeros(1000, length(time_cream));
       Temp = zeros(1000, length(time_cream));
       % Because ode45 does not use a constant change in time, must fill in the
       % matrix with extra points to fulfill dimensional equivalence
       % Because each coffee starts at t=0 at Tinit, this is the point that is
       % duplicated in the column vector for each different cream-adding
       % simulation
       extra_t = zeros(1000, 1);
       extra_temp = Tinit * ones(1000, 1);
       hold on
       for i = 1 : length(time_cream)
          % set mass of mix to be coffee; it will be changed to mass of mix
          % later when cream is added. Must change back for new cream adding
          m_mix = m_coffee;
          t_cream = time_cream(i);
          % Calculate internal energy of the coffee before cream addition
          [t_1, u_1] = ode45(@coffee_U_change, [time_init, t_cream], ...
               Tinit * m_coffee * c);
           % Calculate cooling process after cream is added
          [t_2, u_2] = ode45(@coffee_U_change, [t_cream, time_final], ...
              energy_coffee_cream(u_1(end)));
          % Combine time series into one column vector for each cream addition
          T(:, i) = [extra_t(1 : length(T) - length([t_1; t_2])); t_1; t_2];
          Temp(:, i) = [extra_temp(1 : length(T) - length([u_1; u_2]));...
              u_1 / m_coffee / c; u_2 / m_mix / c];
       end

       % Plot each simulation on one graph
       % Possible to graph on different figures, but make sure temp_cream
       % is a reasonable amount of figures
       for i = 1 : length(T(1,:))
          hold on
          plot(T(:, i), Temp(:, i));
          axis([time_init, time_final, Troom - 5, Tinit])
       title('Temperature of Coffee in Ceramic Mug Over Time')
       xlabel('Time (seconds)')
       ylabel('Temperature (Kelvin)')
       end
    end
end