function [theta, s] = mass_eval()

s_mass = 0;
t_amb = 20;
dt = 1.0;
Td = 10;
N = 100;

P_max = 3000;
P = 1000;
theta = t_amb;


for n=1:N,
	duty = mod(n*dt/Td*P_max, P_max);
	heat = (duty < P);
	Pp = P_max * heat;
	[theta, s_mass] = mass(s_mass, dt, 4190, 20, 0.1, Pp, 10, theta-t_amb);
	theta_(n) = theta;
	Pp_(n) = Pp;
end

t = (0:N-1)*dt;

subplot (2,1,1)
plot(t, theta_); grid;
subplot (2,1,2)
plot(t, Pp_); grid;

endfunction