function [f qc] = calcf(t,x,zc,qbck,qp,hBit,Va,VaDot,Vd,betaa,betad,M,...
    theta1,theta2,rhod,rhoa,g)
%Function calculates right hand side f = [f1;f2;f3] of system equations
%x = [pc;pp;qbit] 
%?NS - 2009
%Modified: DEH - 2010

% get states
pc = x(1);
pp = x(2);
qbit = x(3);

% calc qc
qc = 0.0025*zc*sqrt(2/rhoa*(pc-0));

% calc f
f1 = (betaa/Va) * (qbit - qc + qbck + VaDot);
f2 = (betad/Vd) * (qp - qbit);

% project
if qbit>0
    f3 = (1/M) * (pp-pc-theta1*qbit-theta2*abs(qbit)*qbit +...
        (rhod-rhoa)*g*hBit);
else
    f3 = max((1/M) * (pp-pc+(rhod-rhoa)*g*hBit),0);
end

f = [f1;f2;f3];
end