function [out1, out2, out3] = rate_state(t,x,flag,const)
%%       [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
	d_c1 = const(1);
	d_c2 = const(2);
	A1 = const(3); 
	B1 = const(4); 
	A2 = const(5); 
	B2 = const(6); 
	sigma1 = const(7);
    sigma2 = const(8);
    
	k1 = const(9); 
	k2 = const(10); 
	v_inf = const(11);
	eta = const(12);


    v1     = x(1);
	theta1 = x(2);
	tau1   = x(3);
    v2     = x(4);
	theta2 = x(5);
	tau2   = x(6);

     %return rates
     out1(2) = 1-theta1*v1/d_c1; %theta1
     out1(3) = k1*(v_inf - v1)-k1*(v1 - v2);  %tau1
     out1(1) =  (out1(3)/sigma1 - B1*out1(2)/theta1 ) /(eta/sigma1 + A1/v1); % v1
     
     out1(5) = 1-theta2*v2/d_c2; %theta2
     out1(6) = k2*(v1 - v2);  %tau2
     out1(4) =  (out1(6)/sigma2 - B2*out1(5)/theta2 ) /(eta/sigma2 + A2/v2); % v2
     
     out1 = out1';

