% coffee_cooling_sim.m  
% Simulate coffee cooling in different mugs using Newton's Law of Cooling  
% Designed for beginners in MATLAB programming.

%% Clear workspace and command window
clear variables;
clc;

%% Common parameters
T0        = 85;       % initial coffee temperature in °C
T_room    = 22;       % ambient temperature in °C
m         = 0.35;     % mass of coffee in kg
Cp        = 4186;     % specific heat J/(kg·°C)
A         = 0.03;     % exposed surface area in m²
tspan     = [0 1800]; % simulation time span in seconds (30 minutes)

%% Heat‑transfer coefficients for different mugs
h1 = 5;   % W/(m²·°C) – "insulated mug"
h2 = 10;  % W/(m²·°C) – "standard mug"
h3 = 20;  % W/(m²·°C) – "thin mug"

%% Pre‑compute time constants τ for each
tau1 = (m * Cp) / (h1 * A);
tau2 = (m * Cp) / (h2 * A);
tau3 = (m * Cp) / (h3 * A);

fprintf('Tau values (in seconds):\n');
fprintf('   h = %.1f → τ = %.1f s (%.2f min)\n', h1, tau1, tau1/60);
fprintf('   h = %.1f → τ = %.1f s (%.2f min)\n', h2, tau2, tau2/60);
fprintf('   h = %.1f → τ = %.1f s (%.2f min)\n', h3, tau3, tau3/60);

%% Define ODE as anonymous functions for each mug
ode1 = @(t,T) -(1/tau1) * (T - T_room);
ode2 = @(t,T) -(1/tau2) * (T - T_room);
ode3 = @(t,T) -(1/tau3) * (T - T_room);

%% Solve ODEs with ode45
[t1, T1] = ode45(ode1, tspan, T0);
[t2, T2] = ode45(ode2, tspan, T0);
[t3, T3] = ode45(ode3, tspan, T0);

%% Plot results
figure;
plot(t1/60, T1, 'LineWidth', 2, 'DisplayName', 'Insulated mug (h=5)');
hold on;
plot(t2/60, T2, 'LineWidth', 2, 'DisplayName', 'Standard mug (h=10)');
plot(t3/60, T3, 'LineWidth', 2, 'DisplayName', 'Thin mug (h=20)');
hold off;

xlabel('Time (minutes)');
ylabel('Temperature (°C)');
title('Cooling of Coffee in Different Mugs');
legend('Location','SouthWest');
grid on;

%% Short written answers (you can copy into a report)
disp(' ');
disp('Written answers:');
disp('1. Which mug kept the coffee warm the longest? Why?');
disp('   Answer: The insulated mug (lowest h) kept it warm the longest because it has the largest time constant τ → slower cooling.');
disp(' ');
disp('2. What happens to τ as h increases?');
disp('   Answer: τ = (m·Cp)/(h·A). As h increases, τ decreases → faster cooling.');
disp(' ');
disp('3. How would adding a lid (reducing A) change the outcome?');
disp('   Answer: Reducing A increases τ (because denominator h·A is smaller) → slower cooling → coffee stays warmer longer.');
disp(' ');
disp('4. Suggest one way to improve the realism of this model.');
disp('   Answer: For example: include heat loss through sides and bottom of the mug (not just top surface), or let ambient temperature vary, or include radiation and convection separately.');
