
% Read data from file
% ---------------------------------------------
infile = 'data_varmepump.dat';   % Name of data file
data = load(infile, '-ascii');   % Read data
time = data(:,1);                % Time
temp_hot = data(:,2);            % Temperature in hot reservoir
temp_cold = data(:,3);           % Temperature in cold reservoir
compressor_pulses = data(:,4);   % Pulses of energy from compressor

%%Plot temperature in the hot and cold reservoir
% ---------------------------------------------
figure(1)
plot(time,temp_hot,'-b');
hold on;
plot(time,temp_cold,'--r');
title(['Temperatures']);
xlabel(['Time (s)']);
ylabel(['Temperature (' char(176) 'C)']);
xlim([0 1400]);
ylim([0 45]);
legend({'Hot reservoir','Cold reservoir'},'Location','northwest')

%%Calculate effect by compressor
% ---------------------------------------------
work = cumsum(compressor_pulses)*4.5; % Convert the pulses of energy 
                                      % into total work done by compressor.
                                      % Each pulse is 4,5 kJ.
P_w = polyfit(time,work,1);           % Fit linear polynomial to curve
P_in = P_w(1);                        % Effect of compressor

%%Plot work done by compressor
% ---------------------------------------------
figure(2)
plot(time,work,'-b');
hold on;
plot(time,P_w(1)*time+P_w(2),'--r');
title(['Work done by compressor']);
xlabel(['Time (s)']);
ylabel(['Work W (kJ)']);
xlim([0 1400]);
ylim([0 450]);
legend({'Measured','Linear fit'},'Location','northwest')

% Calculate actual heat factor
% ---------------------------------------------
P = polyfit(time,temp_hot,2);  % Fit polynomial
dT_dt = P(2) + 2*P(1)*time;    % Time derivative
m = 10;                        % Mass of hot reservoir [kg]
c = 4.180;                     % Heat capacity [J/(kg K)]
P_out = m*c*dT_dt;             % Effect out
heat_factor_exp = P_out/P_in;  % Heat factor

% Calculate ideal theoretical heat factor
% (remember to convert temperature to units of Kelvin)
% ---------------------------------------------
heat_factor_th = ( temp_hot + 273.15) ./ ( temp_hot - temp_cold );

%%Plot heat factors
% ---------------------------------------------
figure(3)
plot(time,heat_factor_exp,'-b');
hold on
plot(time,heat_factor_th,'--r');
title(['Heat factor']);
xlabel(['Time (s)']);
ylabel(['Heat factor V_f']);
xlim([0 1400]);
ylim([0 30]);
legend({'Practical','Theoretical'},'Location','northeast')
