function [allslips,alltimes]=get_ratestate(d_c,A,B,sigma,theta0,tau0,preslip,sigrate,coslip,costress,stressG,stopt)

%Note: numsteps is number of time steps(to be determined by this algorithm)
%This program assumes the slip rate on the fault before the earthquake is
%   known (preslip) -- then, either the initial state (theta0) OR the initial stress
%   (tau0) must be specified, and the other parameter is given by the
%   rate-state equations.
%Assign empty matrix to initial value that is solved for. For example,
%   tau0=[] or theta0=[]. If tau0=[] AND theta0=[], then tau0 is assumed to
%   given by the steady-state equation for stress


%% OUTPUTS:
%allslips is 2*nafterpatches x numsteps matrix of incremental slip during each
%   time step
%alltimes is "numsteps" length vector with time of each time step 

%% INPUTS:
%costress is 2*npatch x 1 vector with coseismic shear stress change (1st npatch
%   values for strike-slip, second npatch values for dip-slip component
%stressG is 2*npatch x 2*npatch matrix the relates slip (s) to
%   shear stress (tau) on each patch, tau=stressG*s
%stopt is stop-time -- calculates slip from time=0 to time=stopt
%A,B are 2*npatch x 1 vectors of rate-state friction values 
%d_c is 2*npatch x 1 vector of critical slip distance
%theta0 is 2*npatch x 1 vector of initial state values
%tau0 is 2*npatch x 1 vector of initial shear stress (before earthquake)
%preslip is 2*npatch x 1 vector of pre-earthquake slip rates (npath 
%   is number of pathces, 1st napatch values for strike-slip, second npatch
%   values for dip-slip component)
%sigma is 2*npatch x 1 vector of normal stress -- NOTE: Model depends on sigma*A and sigma*B
%   and therefore sigma cannot be determined independently of A and B
%sigrate is 2*npatch x 1 vector of far-field stressing rates



yr2sec=60*60*24*365;
nps=length(d_c);

%nominal coefficient of friction
mu_0=0.6;

% SEISMIC RADIATION DAMPING TERM
   eta = (1.0e+3*1.0e-7);       % units of MPa-yr/m
   
%pre-quake slip rate
index=preslip==0;  %slip rate cannot be zero because of log(sliprate)
preslip(index)=10^-10;
 
if isempty(tau0) & isempty(theta0)  %assume tau0 is given by steady-state equation  
    tau0 = sigma.*(mu_0 + (A-B).*log(preslip));
    theta0 = exp(tau0./(B.*sigma)-mu_0./B-A./B.*log(preslip)+log(d_c));
    else if isempty(theta0)
        theta0 = exp(tau0./(B.*sigma)-mu_0./B-A./B.*log(preslip)+log(d_c));
    end
end
    

%shear stress immediately after the earthquake
tau0=tau0+costress;
%velocity immediately after the earthquake
v0=preslip.*exp(costress./(A.*sigma));


% START AND STOP TIME
   t0      =  0;
   tf      =  stopt;             % time in years

opt = [];
x0      =   [v0; theta0; tau0];

%note ode23_kaj is a modified version of ode23 that kicks out if the slip
%rate becomes very high (essentially if there is an aftershock)
[t,x] = ode23('ratestate_3D', [t0 tf],x0,opt,sigrate,d_c,A,B,sigma,eta,nps,stressG);
 
if length(t)<size(x,1)
    t(end+1)=t(end);
end
delt=diff(t);delt(end+1)=delt(end);delt=repmat(delt,1,nps);
if length(t)<3
delt=diff(t);delt(end+1)=delt(end);delt=delt';delt=repmat(delt,1,nps);
end

sliprate=x(:,1:nps);
allslips=(delt.*sliprate)';
totalslip=sum(allslips,2);
alltimes=t';
