clc
clear
close all

% Loops in calculation
loops = 2;                  % Number of loops
lv = 1:loops*2;             % Loop vector
even = lv(mod(lv,2)==0);    % Even values in loop vector
odd =  lv(mod(lv,2)~=0);    % Odd values in loop vector

% Conditions/constants
t_v = 15;               % Temperature in heat tank [°C]
t_b = 6.5;              % Temperature in buffer tank [°C]
t_v_lim = 0.5;

tau_int = 1;            % Calculation intervals [s]
tau_end = 600*loops;    % Calculation end time [s]

c = 4192;               % Specific heat capacity [J/(kg K)]
m_b = 4;                % Mass in buffer tank [kg]
m_v = 2;                % Mass in heat tank [kg]

P_b = -200;             % Power from heat exchanger in buffer tank [W]
P_v = 400;              % Power from dypkoger i varmetank [W]

pump = 60;              % Pump effect [L/hr]
q_m = pump/3600;        % Mass flow from pump [kg/s]

N = (tau_end/tau_int)+1;    % Number of calculations

% Start values
tv(1) = t_v;                % Temperature in heat tank
tb(1) = t_b;                % Temperature in buffer tank
tau(1) = 0;                 % Start time
tl(1) = t_v+t_v_lim;        % Upper temperature limit of heat tank
ll(1) = t_v-t_v_lim;        % Lower temperature limit of heat tank

% CALCULATION
for i = 2:N;
    tau(i) = tau(i-1)+tau_int;
    tv(i) = (((q_m*c*(tb(i-1)-tv(i-1))+P_v)*tau_int)/(m_v*c))+tv(i-1);
    tb(i) = (((q_m*c*(tv(i-1)-tb(i-1))+P_b)*tau_int)/(m_b*c))+tb(i-1);
    tl(i) = tl(1);
    ll(i) = ll(1);
        if tv(i)<=t_v-t_v_lim
            q_m = 0;
        end
        if tv(i)>=t_v+t_v_lim
            q_m = pump/3600;
        end
        if tau(i) > (odd*tau_end)/(loops*2)
            tv(i) = (((q_m*c*(tb(i-1)-tv(i-1)))*tau_int)/(m_v*c))+tv(i-1);
        end
        if tau(i) > (even*tau_end)/(loops*2)
            tv(i) = (((q_m*c*(tb(i-1)-tv(i-1))+P_v)*tau_int)/(m_v*c))+tv(i-1);
        end
end



plot(tau/60,tv,'k-',tau/60,tb,'r-',tau/60,tl,'b--',tau/60,ll,'b--')