      function x=dvudt(vu)
      global V AA cd At po p pa rho eta h Vuo f ve cdrag rhoa m n me
%
%    this function calculates current dvu/dt in in^3/s 
%    vu is the air volume in in^3          
%
      g=32.174;                    %gravity, ft/sec^2
      h=(V-vu)/AA;
      p=po*(Vuo/vu)^n;
      if p < pa  |  vu > V
		ve=0.0;
      else
		ve=sqrt(2.*g*(p-pa+rho*eta*h)/(12.*rho*(1.-(At/AA)^2)));       
      end
%	  
      if vu < V
		x=cd*At*12.*ve;
      else
		x=0.0;
      end
%

