%Function of the model

function dxdt = model(t,x,u,dist,par,opt) 


%% States
T_f = x(1);
T_r = x(2);

%% Disturbances
price = dist.price.value(dist.price.time <= t);
price = price(end);  

T_o = dist.T_o.value(dist.T_o.time <= t);
T_o = T_o(end);

m_in = dist.m_in.value(dist.m_in.time <= t);
m_in = m_in(end);


%% Inputs
switch 1
    case t < u(1)
        q_f = max(par.UA_fr*(opt.lb.x.T_f-T_r),0);
    case t > u(1) && t <= dist.price.time(2)
        if T_f >= opt.ub.x.T_f
            q_f = max(par.UA_fr*(opt.ub.x.T_f-T_r),0);
        else
            q_f = opt.ub.u.q_f;
        end
    case t > dist.price.time(2);
        if T_f <= opt.lb.x.T_f
            q_f = max(par.UA_fr*(opt.lb.x.T_f-T_r),0);
        else
            q_f = 0;
        end
end

switch 1    
    case t < u(2)
        q_r =max(-par.UA_fr*(T_f-opt.lb.x.T_r) + (par.UA_ro + m_in/par.m_r*par.mCp_r)*(opt.lb.x.T_r-T_o),0);
    case t > u(2) && t <= dist.price.time(2)
        if T_r >= opt.ub.x.T_r
            q_r =max(-par.UA_fr*(T_f-opt.ub.x.T_r) + (par.UA_ro + m_in/par.m_r*par.mCp_r)*(opt.ub.x.T_r-T_o),0);
        else
            q_r = opt.ub.u.q_r;
        end
    case t > dist.price.time(2)
        if T_r <= opt.lb.x.T_r
            q_r =max(-par.UA_fr*(T_f-opt.lb.x.T_r) + (par.UA_ro + m_in/par.m_r*par.mCp_r)*(opt.lb.x.T_r-T_o),0);
        else
            q_r = 0;
        end    
end

%% Model
%Floor
dxdt(1) = (q_f/par.mCp_f)-(par.UA_fr/par.mCp_f)*(T_f-T_r);

%Room
dxdt(2) = (q_r/par.mCp_r)+(par.UA_fr/par.mCp_r)*(T_f-T_r)-(par.UA_ro/par.mCp_r)*(T_r-T_o)...
    +(m_in/par.m_r)*(T_o-T_r);

%Cost function
dxdt(3) = price*(q_f + q_r)/3600;

dxdt = dxdt'*3600;
