function dy = hps(~,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%   This function calculates derivations of HPS states from the linear    %
% mathematical model.                                                     %
%                                                                         %
%   INPUT:                                                                %
%       y  = [hLH hRH hLL hRL pH pL oL oR pp]                             %
%   OUTPUT:                                                               %
%       dy = [dhLH dhRH dhLL dhRL dpH dpL duL duR dpp]                    %
%                                                                         %
%                                                 Last update: 18.05.2010 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global uLmin uLmax uRmin uRmax aLV aRV sLV sRV suL suR auL auR
global alfaLH alfaRH alfaLL alfaRL sLH sRH sLL sRL SL SR 
global betaH betaL gamaH gamaL s0H s0L V0 H0
global pA rok g m R cp T

dy=zeros(9,1);
y(7)=(y(7)-uLmin)/(uLmax-uLmin);
y(8)=(y(8)-uRmin)/(uRmax-uRmin);

for i=1:8
    if (y(i)<0)
            y(i)=0;
    end
end

for i=7:8
    if (y(i)>1)
            y(i)=1;
    end
end

if (y(9)>6e4)
        y(9)=6e4;
elseif (y(9)<(y(5)-pA))
    y(9)=y(5)-pA;
end

QLV=(aLV*sLV*y(7)+auL*suL)*sqrt(2*rok*(y(9)+pA-y(5)));
QRV=(aRV*sRV*y(8)+auR*suR)*sqrt(2*rok*(y(9)+pA-y(5)));

dQLH=2*((y(5)-y(6))/rok+g*y(1));
if (dQLH<0)
    dQLH=0;   
end

dQRH=2*((y(5)-y(6))/rok+g*y(2));
if (dQRH<0)
    dQRH=0;   
end

dQLL=2*((y(6)-pA)/rok+g*y(3));
if (dQLL<0)
    dQLL=0;
end

dQRL=2*((y(6)-pA)/rok+g*y(4));
if (dQRL<0)
    dQRL=0;     
end

QLH=alfaLH*sLH*rok*sqrt(dQLH);
QRH=alfaRH*sRH*rok*sqrt(dQRH);
QLL=alfaLL*sLL*rok*sqrt(dQLL);
QRL=alfaRL*sRL*rok*sqrt(dQRL);

dy(1,1)=1/(rok*SL)*(QLV-QLH);
dy(2,1)=1/(rok*SR)*(QRV-QRH);
dy(3,1)=1/(rok*SL)*(QLH-QLL);
dy(4,1)=1/(rok*SR)*(QRH-QRL);


if (y(5)>=pA)
    Q0H=betaH*s0H/R*sqrt(2*cp/T*y(5)*pA*((pA/y(5))^(2/m-1)-(pA/y(5))^(1/m)));
else
    Q0H=-betaH*s0H/R*sqrt(2*cp/T*y(5)*pA*((y(5)/pA)^(2/m-1)-(y(5)/pA)^(1/m)));  
end

VH=gamaH*V0+SL*(H0-y(1))+SR*(H0-y(2));
dy(5,1)=y(5)/VH*(SL*dy(1)+SR*dy(2))-R*T/VH*Q0H;

if (y(6)>=pA)
    Q0L=betaL*s0L/R*sqrt(2*cp/T*y(6)*pA*((pA/y(6))^(2/m-1)-(pA/y(6))^(1/m)));
else
    Q0L=-betaL*s0L/R*sqrt(2*cp/T*y(6)*pA*((y(6)/pA)^(2/m-1)-(y(6)/pA)^(1/m)));
end

VL=gamaL*V0+SL*(H0-y(3))+SR*(H0-y(4));
dy(6,1)=y(6)/VL*(SL*dy(3)+SR*dy(4))-R*T/VL*Q0L;