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

4 vues (au cours des 30 derniers jours)
Chris
Chris le 7 Déc 2021
Modifié(e) : Matt J le 7 Déc 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);
[ux,uy,uz] = calc_okada(1,X(:)-x0,Y(:)-y0,.25,dip,z0,dL,dW,2,strike);
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)
pcolor(X/1e3,Y/1e3,reshape(data(1:np),size(X))),shading flat
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)
pcolor(X/1e3,Y/1e3,reshape(data(np+1:end),size(X))),shading flat
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)
pcolor(X/1e3,Y/1e3,reshape(prediction(1:np),size(X))),shading flat
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)
pcolor(X/1e3,Y/1e3,reshape(prediction(np+1:end),size(X))),shading flat
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);
[Date,East,North,Vertical] = readGPSfile(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')

Réponses (0)

Catégories

En savoir plus sur Earthquake Engineering dans Help Center et File Exchange

Produits

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by