
% matlab script to plot velocity, stress, and strain models of Imperial fault


% compute ANALYTIC velocity, strain, and stress profiles

% compute Weertman 2-D arctan  model
	
     xp=-512:1:512;
     depth=6;
     Vo=-40;
     Vanalytic=Vo/pi*(atan(xp./depth));


% calculate analytic strain and stress from 2-D arctan model 
% convert mm to radians (km) by multiplying by 10^-6;

      unit_conv=10^-6;
      shear_mod=30*10^9;
      Strain_analytic=unit_conv.*-Vo/pi*(depth./( (xp.*xp) + (depth.*depth)));
      Stress_analytic=shear_mod.*Strain_analytic;



% next load the SCEC data and the models

%  load the data file from grdtrack results

     load scec_gps.dat
     xkm=scec_gps(:,1)-512;
     ykm=scec_gps(:,2);
     vy=scec_gps(:,4);
     vys=scec_gps(:,6);
 
 
% load velocity model grd

 yv_model=grdread('ytotvel_2002_2003.grd');
 velocity=yv_model(1:2048,512:1536);
 

% load stress grd  (Dave's model), comparison to Freed,Pa/yr
% projection isn't quite right (new stress model, old velocity model), 
%    shift grid by 8 km to match with velocity grid

zstress=grdread('zstress.grd');
stress=zstress(1:2048,8:1024+8);     % in Pa/yr


% compute strain from stress (divide by shear mod)
shear_mod=30*10^9;
strain=stress./shear_mod;    % in radians/yr


% plot map view of velocity  & strain models with gps data 
 
 figure(1)
 subplot(1,3,1),imagesc(flipud(velocity))
 colorbar
 grid
 hold on
 plot(xkm,2048-ykm,'k*');
 title('velocity mm/yr')
 subplot(1,3,2),imagesc(flipud(strain),[0,3e-6]) %  plotting microrad/yr
 colorbar
 grid
 hold on
 plot(xkm,2048-ykm,'k*');
 title('strain rate microrad/yr')
 subplot(1,3,3),imagesc(flipud(stress./10^3),[0,80]) % convert to kPa/yr
 colorbar
 grid
 hold on
 plot(xkm,2048-ykm,'k*');
 title('stress rate kPa/yr')


% extract profile lines for velocity and strain models
% imperial fault location at 2048-1700 (y=348) (matlab flips grid)

 imp_loc=348;

 imp_vel_prof=velocity(imp_loc,:);
 imp_strain_prof=strain(imp_loc,:);
 imp_stress_prof=stress(imp_loc,:);

 x_iprof_map=imp_loc*ones(1,length(imp_vel_prof));
 x_iprof=1:length(imp_vel_prof);
 
 
% grab gps points around model profile (+/- 100 km)

 imp_ymax=2048-1600;
 imp_ymin=2048-1800;
 
 j=1;
   for i=1:length(xkm)
 	if imp_ymin <= ykm(i)
 	if imp_ymax >= ykm(i)
 		ykm(i);
 		igps_bin_xkm(j)=xkm(i);
 		igps_bin_ykm(j)=ykm(i);
 		igps_bin_vy(j)=vy(i);
 		igps_bin_vys(j)=vys(i);
 		j=j+1;
 	end;
 	end;
 end
  
% calculate rms (S&S velocity model - GPS velocity) for velocity shift

	for i=1:length(igps_bin_xkm);
   		SS_model_track(i)=velocity(round(igps_bin_ykm(i)),round(igps_bin_xkm(i)));	
	end

SS_diff=igps_bin_vy-SS_model_track;
SS_rms=sqrt(mean(SS_diff.*SS_diff));
 

 
% plot S&S model, 2-D arctangent analytic model 
% apply linear shift (2 mm/yr)to analytic model for common reference frame
% could also have done this to the S&S model and gps data


	figure(2)
	orient tall
	subplot(2,1,1),errorbar(igps_bin_xkm-8-512,igps_bin_vy-SS_rms,igps_bin_vys,'ko')
	hold on
	plot(x_iprof-512,Vanalytic-2,'b')
	hold on
	plot(x_iprof-8-512,imp_vel_prof,'--c')
	axis([-60 60 -25 25])
	legend('GPS data', '2-D analytic velocity model','smith-konter velocity model')
	xlabel('distance (km)')
	ylabel('horizontal velocity (mm/yr)')
	title('Imperial velocity profile')
	
	
% plot strain rate, was in rad/yr, convert to microrad/yr
	subplot(2,1,2),	plot(x_iprof-513,Strain_analytic*10^9,'b')   
	axis([-60 60 0 2500])
	xlabel('distance (km)')
	ylabel('strain rate (nanostrain/yr)')
	title('Imperial strain rate profile')
	hold on
	plot(x_iprof-13-512,imp_strain_prof*10^9,'--c')	 
	legend( '2-D analytic shear strain rate','smith-konter strain rate')

