# How can I make arrows in Patches by using these codes?

8 views (last 30 days)
Chris on 7 Dec 2021
Edited: Matt J on 7 Dec 2021
<Code 1>
%version that explicitly performs the inversion instead of just calling
%lsqlin
%set up some distances for a 1-dimensional dipping fault problem
maxx = 50e3; %meters
dx = 1e3; %meters
[X,Y] = meshgrid(0:dx:maxx*2,-maxx:dx:maxx); %make grid of X and Y values
xt = 0; %location of top of trench
strike = 0; %Fault strikes pure north-south, perpendicular to our profile
dip = 30; %degrees, angle of subduction zone from horizontal
L = 50e3; %Length of fault (along strike)
z1 = 0e3; %depth of top of locked zone, in meters
z0 = 50e3; %depth of bottom of locked zone, in meters
W = (z0-z1)/sind(dip); %Down-dip width of fault
Wp = 10; %Number of fault patches along width (downdip)
Lp = 10; %Number of fault patches along length(along strike)
npatch = Wp*Lp; %Number of fault patches
x1 = xt+z1*cosd(dip); %horiz. location of top of locked zone
x0 = xt+z0*cosd(dip); %horiz. location of bottom of locked zone
nu = 0.25; %Poisson's ratio
fault_type = 2; %dip slip fault
%Here we are starting by inverting for coseismic slip, since that is easier
%than interseismic.
dL = L/Lp;
dW = W/Wp;
lengths = (-L/2+dL/2):dL:(L/2-dL/2);
depths = dW*sind(dip):dW*sind(dip):(W*sind(dip));
[z0p,tmpL] = meshgrid(depths,lengths); %position of center base (downdip side) of each patch in depth and along-strike coords
Xp=tmpL*sind(strike)+z0p/tand(dip)*cosd(strike);
Yp=tmpL*cosd(strike)-z0p/tand(dip)*sind(strike);
dX1=dL/2*sind(strike);
dY1=dL/2*cosd(strike);
dX2=dW*cosd(dip)*cosd(strike);
dY2=dW*cosd(dip)*sind(strike);
dZ=dW*sind(dip);
id = 0;
faultx = [];
faulty = [];
faultz = [];
for i=1:npatch
faults(i,:) = [Xp(i) Yp(i) z0p(i) dL dW strike dip];
faultx = [faultx;Xp(i)-dX1-dX2 Xp(i)-dX1 Xp(i)+dX1 Xp(i)+dX1-dX2 Xp(i)-dX1-dX2];
faulty = [faulty;Yp(i)-dY1+dY2 Yp(i)-dY1 Yp(i)+dY1 Yp(i)+dY1+dY2 Yp(i)-dY1+dY2];
faultz = [faultz;z0p(i)-dZ z0p(i) z0p(i) z0p(i)-dZ z0p(i)-dZ];
end
np=size(X(:),1)
greenY = zeros(np,npatch);
greenZ = zeros(np,npatch);
greenX = zeros(np,npatch);
smooth=smoother(Lp,Wp,1,1);
%make greens
for i=1:npatch
x0 = faults(i,1);
y0 = faults(i,2);
z0 = faults(i,3);
L = faults(i,4);
W = faults(i,5);
strike = faults(i,6);
dip = faults(i,7);
greenX(:,i) = ux;
greenY(:,i) = uy;
greenZ(:,i) = uz;
end
G = [greenX;greenY]; %like what we would use if we only had horizontal GPS data
trueslip=11; %true slip on one patch in the earthquake
slip=zeros(npatch,1); %now need a vector, one for each patch.
slip(23)=trueslip; %only put slip on 5th fault patch
noise = 0.05; %the standard deviation of the noise we add to the data
data_nonoise = G*slip;
data = data_nonoise+noise*randn(size(data_nonoise)); %add normally-distributed (Gaussian) noise to the vertical component of deformation.
%%% solve for best-fit rate
%%% The uz matrix is our model of what ONE unit of slip looks like - so
%%% this inverse problem is actually pretty simple. We'll still use a
%%% built-in matlab routine to solve it. lsqlin stands for "least-squares-linear
%Now we add smoothing - we will start with just "damping" rather than
%spatial smoothing. The variable "lambda" controls the weight/strength of
%smoothing applied to the solution
smoo = eye(npatch); % "eye" makes a diagonal matrix of ones
lambda = 0.5;
Ga = [G;lambda*smoo]; %Here we "augment" or add to the G matrix.
%model = lsqlin(uz,uz_data);
invmatrix = inv(Ga'*Ga)*G';
model = invmatrix*data;
prediction = G*model;
figure
subplot(3,2,1)
patch(faultx'/1e3,faulty'/1e3,faultz'/1e3,slip)
set(gca,'zdir','reverse')
view(20,35)
axis image
grid on
xlabel('East (km)')
ylabel('North (km)')
zlabel('Depth (km)')
title('input slip')
colorbar
subplot(3,2,2)
patch(faultx'/1e3,faulty'/1e3,faultz'/1e3,model)
set(gca,'zdir','reverse')
view(20,35)
axis image
grid on
xlabel('East (km)')
ylabel('North (km)')
zlabel('Depth (km)')
title('inferred slip')
colorbar
subplot(3,2,3)
xlabel('East (km)')
ylabel('North (km)')
colorbar
title('X-component of deformation, input data')
hold on
plot(faultx'/1e3,faulty'/1e3,'k')
cx=caxis;
subplot(3,2,4)
xlabel('East (km)')
ylabel('North (km)')
colorbar
title('Y-component of deformation, input data')
hold on
plot(faultx'/1e3,faulty'/1e3,'k')
cy=caxis;
subplot(3,2,5)
xlabel('East (km)')
ylabel('North (km)')
colorbar
title('X-component of deformation, predicted data')
hold on
plot(faultx'/1e3,faulty'/1e3,'k')
caxis(cx)
subplot(3,2,6)
xlabel('East (km)')
ylabel('North (km)')
colorbar
title('Y-component of deformation, predicted data')
hold on
plot(faultx'/1e3,faulty'/1e3,'k')
caxis(cy)
<Code 2>
clear all
Eqdate=datetime('18Jan23', 'Format', 'yyMMMdd', 'InputFormat', 'yyMMMdd');
files=dir('gps*');
for i=1:length(files)
disp(files(i).name);
%make a quick figure for each site
figure
plot(Date,Vertical)
hold on
yax=ylim;
plot([Eqdate Eqdate],yax,'r') %magenta vertical line at time of earthquake
title(files(i).name)
X(i)=East(1);%set overall location of point to first E/N point (used in calc_okada)
Y(i)=North(1);
%the following approach does not work if you have an earthquake in the
%middle of a huge gap in your data! Then you will see the coseismic +
%interseismic deformation.
before = find(Date<Eqdate); %indices of points in time at least one day before the earthquake
after = find(Date>Eqdate); %indices of points in time at least one day after the earthquake
if(~or(isempty(before),isempty(after))) %only true if both before and after are defined
%find single day before and after
before = before(end);
after = after(1);
GPSdx(i)=East(after)-East(before);
GPSdy(i)=North(after)-North(before);
GPSdz(i)=Vertical(after)-Vertical(before);
else
GPSdx(i)=NaN;
GPSdy(i)=NaN;
GPSdz(i)=NaN;
end
end
figure
quiver(X,Y,GPSdy,GPSdz,'k')
title('EQ displacements-dy,dz')