function out = ratestate_3D(t,x,flag,sigrate,d_c,A,B,sigma,eta,nps,G)

%%       [out1, out2, out3] = rate_state(t,x,flag,const)
%
% S-file version of fault-constitutive law
%       t = times
%       x = model vector%
%	x(1) is velocity
%	x(2) is state  variable
%	x(3) is shear stress
%
%  unwrap constants
	
    v     = x(1:nps);
    theta = x(nps+1:2*nps);
	tau   = x(2*nps+1:end);
   
     %return rates
    out=zeros(1,3*nps);
     out(nps+1:2*nps) = 1-theta.*v./d_c;
     out(2*nps+1:3*nps) = G*v+sigrate;
     out(1:nps) =  (out(2*nps+1:3*nps)'./sigma - B.*out(nps+1:2*nps)'./theta ) ./(eta./sigma + A./v);
     
     out = out';