function cool_off (initial_temp, timestep, end_time)
%
% Compute the temperature of an object as a function of time.
%   The object begins at a temperature of "initial_temp"
%   and cools toward an ultimate temperature of absolute zero.
%   Keep going until we reach the time "end_time."
%
% We use Newton's law of cooling.
%
%        dT / dt  =   - k*T
%
%   where "k" is the coefficient of cooling, and is
%   equal to the reciprocal of the time constant.
%   We set the time constant in the program's initialization
%   stage.  The calculations employ Euler's method.
%
% The program prints to the screen one line per timestep, giving
%   the time (in seconds) and the temperature (in Kelvin), like this:
%
%          t    91.000  temp    23.328
%
%
% Usage:  cool_off (initial_temp, timestep, end_time)
%
%    where
%           initial_temp       is starting temperature of the
%                                  object (Kelvin)
%
%           timestep           is the timestep to use in the integration
%                                  (seconds)
%
%           end_time           is the point at which we stop the
%                                  calculations (seconds)
%
%
% MWR 11/29/2011

% set this to 1 to see diagnostic messages
debug = 0;

% the temperature of the object (Kelvin)
temp = initial_temp;
% the timestep for integration (seconds)
dt = timestep;
% time since the object started cooling (seconds)
t = 0;
% time constant for cooling (seconds)
time_constant = 10.0;



if (debug > 0)
    fprintf(1, 'initial_temp is %f \n', initial_temp);
    fprintf(1, 'timestep is %f \n', timestep);
    fprintf(1, 'end_time is %f \n', end_time);
end

while (t <= end_time)
   fprintf(1, 't %9.3f  temp %9.3f \n', t, temp);

   % how much will the object cool during this coming timestep?
   dtemp = - (1.0 / time_constant) * temp;
   if (debug > 0)
       fprintf(1, '  dT  %9.3f \n', dtemp);
   end

   % use Euler's method to compute the new temperature
   new_temp = temp + dtemp;

   % update variables for the next iteration
   temp = new_temp;
   t = t + dt;


end

Published with MATLAB® 7.13