function [sys,x0,str,ts] = Padova_model(t,x,u,flag,parameters)
% Attenzione: x(1)=stomach1 (mg)
%             x(2)=stomach2 (mg)
%             x(3)=intestine (mg)
%             x(4)=plasma glucose + insulin-independent tissues (mg/kg)
%             x(5)=insulin-dependent tissues (mg/kg)
%             x(6)=plasma insulin (pmol/kg)
%             x(7)=insulin action, X  (pmol/L)
%             x(8)=delay compartment for insulin action on glucose production (pmol/L)
%             x(9)=insulin action on glucose production, Id (pmol/L)
%             x(10)=insulin in the liver (pmol/kg)
%             x(11)=1st compartment of subcutaneous insulin (pmol/kg)
%             x(12)=2nd compartment of subcutaneous insulin (pmol/kg)
%             x(13)=subcutaneous glucose (mg/kg)
% Sudden change in dynamics
% switch parameters.val
%     case 2 %Test data
%         if t<(4*24*60)
%             % Dynmaic Mode B
%             parameters.p2u = parameters.p2u*4;
%             parameters.ki = parameters.ki*2;
%             parameters.k1 = parameters.k1*2;
%             parameters.kp2 = parameters.kp2;%*0.5;
%             parameters.kp3 = parameters.kp3*0.2;
%         end
%     case 1 % Training Data
%         if t>(4*24*60)
%             % Dynmaic Mode B
%             parameters.p2u = parameters.p2u*4;
%             parameters.ki = parameters.ki*2;
%             parameters.k1 = parameters.k1*2;
%             parameters.kp2 = parameters.kp2;%*0.5;
%             parameters.kp3 = parameters.kp3*0.2;
%         end        
% end
% parameters.k1 = parameters.k1*(1+0.2*randn(1));
% parameters.k2 = parameters.k2*(1+0.2*randn(1));
% parameters.p2m = parameters.k1*(1+0.2*randn(1));
% parameters.ki = parameters.ki*(1+0.2*randn(1));
% %parameters.m1 = parameters.m1*(1+0.2*randn(1));
% parameters.m30 = parameters.m30*(1+0.2*randn(1));
% parameters.p2u = parameters.p2u*(1+0.2*randn(1));
% parameters.ksc = parameters.ksc*(1+0.2*randn(1));

switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
   [sys,x0,str,ts]=mdlInitializeSizes(parameters);

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u,parameters);
   
  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3,                                               
    sys = mdlOutputs(t,x,u,parameters);    

  %%%%%%%%%%%%%%%%%%%
  % Unhandled flags %
  %%%%%%%%%%%%%%%%%%%
  case { 2, 4, 9 },
    sys = [];

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);

end

%end sfundsc1

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts]=mdlInitializeSizes(parameters)
sizes = simsizes;

sizes.NumContStates  = 13;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 21;
sizes.NumInputs      = 3;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;



sys = simsizes(sizes);
x0=parameters.x0;
str = [];
ts  = [0 0]; 

% end mdlInitializeSizes

%
%=======================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=======================================================================
%
function sys=mdlDerivatives(t,x,u,parameters)
% ABSORPTION
%glucose in the stomach
qsto=x(1)+x(2);
parameters.dosekempt=u(3)+parameters.dose(find(parameters.Tdose<=t,1,'last'));
%Stomach solid
sys(1)=-parameters.kmax*x(1)+u(1);
if parameters.dosekempt>0
    aa=5/2/(1-parameters.b)/parameters.dosekempt;
    cc=5/2/parameters.d/parameters.dosekempt;
    kgut=parameters.kmin+(parameters.kmax-parameters.kmin)/2*(tanh(aa*(qsto-parameters.b*parameters.dosekempt))-tanh(cc*(qsto-parameters.d*parameters.dosekempt))+2);
%     if t>=parameters.Tdose(1)
%         figure(3)
%         plot(parameters.dosekempt, parameters.kmax,'*')
%         hold on
%         plot(qsto, kgut,'o')
%         %pause(0.01)
%     end;
else
    kgut=parameters.kmax;
end;
%stomach liquid
sys(2)=parameters.kmax*x(1)-x(2)*kgut;
%intestine
sys(3)=kgut*x(2)-parameters.kabs*x(3);

%Rate of appearance
Rat=parameters.f*parameters.kabs*x(3)/parameters.BW;
%Glucose Production
EGPt=parameters.kp1-parameters.kp2*x(4)-parameters.kp3*x(9);
%Glucose Utilization
Uiit=parameters.Fsnc;

% RENAL EXCRETION ;
if x(4)>parameters.ke2
     Et=parameters.ke1*(x(4)-parameters.ke2); %equazione 27
else Et=0;
end

% GLUCOSE KINETICS; 
sys(4)=max(EGPt,0)+Rat-Uiit-Et-parameters.k1*x(4)+parameters.k2*x(5);

    Vmt=parameters.Vm0+parameters.Vmx*x(7);
    Kmt=parameters.Km0; %
    Uidt=Vmt*x(5)/(Kmt+x(5));
    sys(5)=-Uidt+parameters.k1*x(4)-parameters.k2*x(5);

Gt=x(4)/parameters.Vg;

% INSULIN KINETICS
sys(6)=-(parameters.m2+parameters.m4)*x(6)+parameters.m1*x(10)+parameters.ka1*x(11)+parameters.ka2*x(12);
It=x(6)/parameters.Vi;

% INSULIN ACTION ON GLUCOSE UTILIZATION
sys(7)=(-parameters.p2u*x(7)+parameters.p2u*(It-parameters.Ib));

%NSULIN ACTION ON PRODUCTION
sys(8)=-parameters.ki*(x(8)-It);
sys(9)=-parameters.ki*(x(9)-x(8));


%INSULIN IN THE LIVER (pmol/kg)
sys(10)=(-(parameters.m1+parameters.m30)*x(10)+parameters.m2*x(6));

%(CSUBCUTANEOUS INSULIN KINETICS Modified by Chiara following Volund)(INPUT 2 IS INSULIN) 
sys(11)=(u(2)-(parameters.ka1+parameters.kd)*x(11)); % I used the parameters.u2ss as external input, thus NOW your controller has to provide the entire pump
sys(12)=(parameters.kd*x(11)-parameters.ka2*x(12));

%SUBCUTANEOUS GLUCOSE 
sys(13)=((-parameters.ksc*x(13)+parameters.ksc*x(4))); 


% end mdlDerivatives

%
%=======================================================================
% mdlOutputs
% Return the output vector for the S-function
%=======================================================================
%
function sys = mdlOutputs(t,x,u,parameters)
sys(1)=x(13)/parameters.Vg;%subcutanoeous glucose
sys(2)=x(4)/parameters.Vg;%plasma glucose
sys(3)=parameters.kp1-parameters.kp2*x(4)-parameters.kp3*x(9);%Production
sys(4)=parameters.f*parameters.kabs*x(3)/parameters.BW;%Rate of Appearance
Vmt=parameters.Vm0+parameters.Vmx*x(7);
Kmt=parameters.Km0; %
Uidt=Vmt*x(5)/(Kmt+x(5));
sys(5)=parameters.Fsnc+Uidt;%Utilization
sys(6)=x(11)*parameters.ka1+x(12)*parameters.ka2;% Rate of appearance dell'insulina
sys(7) =x(6)/parameters.Vi;%plasma Insulin
sys(8:20)=x(1:13);
sys(21)=x(1)+x(2);%Qsto
