function dy = hpsnr(~,y,const,dy)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%   This function calculates derivations of HPS states from the nonlinear %
% 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: 07.01.2011 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

for i=1:4
    if (y(i)>0.3499)
            y(i)=0.3499;
    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)-const(27)))
    y(9)=y(5)-const(27);
end

dpp=sqrt(2*const(28)*(y(9)+const(27)-y(5)));
QLV=(const(1)*const(3)*y(7)+const(7)*const(5))*dpp;
QRV=(const(2)*const(4)*y(8)+const(8)*const(6))*dpp;

dQLH=2*((y(5)-y(6))/const(28)+const(29)*y(1));
if (dQLH<0)
    dQLH=0;   
end

dQRH=2*((y(5)-y(6))/const(28)+const(29)*y(2));
if (dQRH<0)
    dQRH=0;   
end

dQLL=2*((y(6)-const(27))/const(28)+const(29)*y(3));
if (dQLL<0)
    dQLL=0;
end

dQRL=2*((y(6)-const(27))/const(28)+const(29)*y(4));
if (dQRL<0)
    dQRL=0;     
end

QLH=const(9)*const(13)*const(28)*sqrt(dQLH);
QRH=const(10)*const(14)*const(28)*sqrt(dQRH);
QLL=const(11)*const(15)*const(28)*sqrt(dQLL);
QRL=const(12)*const(16)*const(28)*sqrt(dQRL);


if (y(5)>=const(27))
    Q0H=const(19)*const(23)/const(31)*sqrt(2*const(32)/const(33)*y(5)*const(27)*((const(27)/y(5))^(2/const(30)-1)-(const(27)/y(5))^(1/const(30))));
else
    Q0H=-const(19)*const(23)/const(31)*sqrt(2*const(32)/const(33)*y(5)*const(27)*((y(5)/const(27))^(2/const(30)-1)-(y(5)/const(27))^(1/const(30))));  
end

if (y(6)>=const(27))
    Q0L=const(20)*const(24)/const(31)*sqrt(2*const(32)/const(33)*y(6)*const(27)*((const(27)/y(6))^(2/const(30)-1)-(const(27)/y(6))^(1/const(30))));
else
    Q0L=-const(20)*const(24)/const(31)*sqrt(2*const(32)/const(33)*y(6)*const(27)*((y(6)/const(27))^(2/const(30)-1)-(y(6)/const(27))^(1/const(30))));
end

VH=const(21)*const(25)+const(17)*(const(26)-y(1))+const(18)*(const(26)-y(2));
VL=const(22)*const(25)+const(17)*(const(26)-y(3))+const(18)*(const(26)-y(4));
prac=const(31)*const(33);


dy(1,1)=1/(const(28)*const(17))*(QLV-QLH);
dy(2,1)=1/(const(28)*const(18))*(QRV-QRH);
dy(3,1)=1/(const(28)*const(17))*(QLH-QLL);
dy(4,1)=1/(const(28)*const(18))*(QRH-QRL);
dy(5,1)=y(5)/VH*(const(17)*dy(1)+const(18)*dy(2))-prac/VH*Q0H;
dy(6,1)=y(6)/VL*(const(17)*dy(3)+const(18)*dy(4))-prac/VL*Q0L;