% Earthquake simulation 
% Kaj Johnson, Indiana University
% last updated 12/8/05
% CONSTITUTIVE PARAMETERS
   d_c1 = 0.01;                % critical displacement (m)
   d_c2 = 0.01;                % critical displacement (m)
   A2 =   0.0033;  %slider 2 is pulled at far vield velocity
   B2 =   0.0045;
   A1 =   0.003;  %slider 1 is pulled by slider 2
   B1 =   0.004;
   mu_0 = 0.6;                % nominal friction coefficient

% BOUNDARY CONDITIONS
   sigma1 = 250;               % normal stress (MPa)
   sigma2 = 250;               % normal stress (MPa)
   k_crit = sigma1*(B1-A1)/d_c1;   
   k1 = 0.1*abs(k_crit);           % stiffness in MPa/m
   L=10000; 
   k2=1.2*sigma2*(B2-A2)/d_c2;;
   %k2=3*10^4*.9/L;   %L is radius of circular slip patch in m
   
   v_inf  =   0.01;           % load-point velocity in m/yr
                               
% INITIAL CONDITIONS
% choose initial state to be immediately after an earthquake

   v02      =   v_inf;
   v01 = 10^-5;
   %steady state
   tau02    =    sigma2*(mu_0 + (A2-B2)*log(v02));
   theta02 = d_c2/v02;
   tau01    =    sigma1*(mu_0 + (A1-B1)*log(v01));
   theta01 = d_c1/v01;
   
%end input section
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
   
   
   % SEISMIC RADIATION DAMPING TERM
   eta = 1.0e+3*1.0e-5;       % units of MPa-yr/m


   %theta0 = exp(tau0/(B*sigma)-mu_0/B-A/B*log(v0)+log(d_c));
   
   x0      =   [v01; theta01; tau01; v02; theta02; tau02];

   const = [d_c1; d_c2; A1; B1; A2; B2; sigma1; sigma2; k1; k2; v_inf; eta];

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

  opt = [];
%  RUN ODE SOLVER
 
  [t,x] = ode23s('double_rate_state', [t0 tf],x0,opt,const);
  u = cumsum(x(:,1));  % displacement
%
%PLOT OUTPUT
%
figure
%subplot(221); plot(t,x(:,1),'b'); title('velocity1 vs time'); hold on
%subplot(222); plot(t,x(:,4),'b'); title('velocity2 vs time'); hold on
%subplot(221); plot(t,x(:,2),'b'); title('theta1 vs time'); hold on
%subplot(221); plot(t,x(:,4),'b'); title('vel2 vs time'); hold on
subplot(221); semilogx(x(:,4), x(:,6),'b'); title('phase plane plot - block2'); hold on
vel=logspace(-5,5);
%tss=sigma2*(mu_0+(A2-B2)*log(vel));
%semilogx(vel, tss,'r:')

subplot(222); plot(t,x(:,5),'b'); title('theta2 vs time'); hold on
%subplot(222); plot(t,d_c2./x(:,4),'r')
delt=diff(t);
delt(end+1)=delt(end);
%subplot(223); plot(t,cumsum(delt.*x(:,1)),'b'); title('disp1 vs time'); hold on
subplot(223); plot(t,x(:,6),'b'); title('tau2 vs time'); hold on
%tau2=sigma2*(mu_0 + (A2-B2)*log(x(:,4)));
%plot(t,tau2,'r')

subplot(224); plot(t,cumsum(delt.*x(:,4)),'b'); title('disp2 vs time'); hold on
%subplot(223); semilogx(x(:,1), x(:,3),'r'); title('phase plane plot - block1'); hold on
%subplot(224); semilogx(x(:,4), x(:,6),'r'); title('phase plane plot - block2'); hold on
%plot steady state velocity
%vel=logspace(-5,5);
%tss=sigma2*(mu_0+(A2-B2)*log(vel));
%semilogx(vel, tss,'r:')

%xlabel('log velocity');ylabel('stress')


figure
subplot(221); semilogx(x(:,1), x(:,3),'b'); title('phase plane plot - block1'); hold on
vel=logspace(-5,5);
%tss=sigma1*(mu_0+(A1-B1)*log(vel));
%semilogx(vel, tss,'r:')

subplot(222); plot(t,x(:,2),'b'); title('theta1 vs time'); hold on
delt=diff(t);
delt(end+1)=delt(end);
subplot(223); plot(t,x(:,3),'b'); title('tau1 vs time'); hold on

subplot(224); plot(t,cumsum(delt.*x(:,1)),'b'); title('disp1 vs time'); hold on

