Anyone can help me. Could you please check my code with attach file for fig 8 page 74 by using equation 2.30 page 67.

2 views (last 30 days)
clc
clear all
close all
syms k1 k2
eta0 = 1; % Maximum displacement(m)
L = 100; % Propagate length (km)
W = 100; % Source width (km)
v = 0.14; % Ruptuer velocity to obtain maximum surface amplitude(km/s)
h = 2; % Water depth (km)
g = 0.0098; % Acceleration due to gravity (km/s^2)
t1 = 50/v; % (s)
t = (t1); % (s)
k = sqrt(k1.^2+k2.^2);
w = sqrt(g.*k.*tanh(k.*h));
% A = ((sin(w*t))./(w.*cosh(k*h))).*((eta0*v)./(2*L));
A = ((cos(w*t))./(cosh(k*h))).*((eta0*v)./(2*L));
B = (1-exp(-1i*100*k1))./(1i*k1);
C = ((exp(-1i*100*k1))./(1-((50*k1)./pi).^2)).*(1i*k1).*((50./pi).^2).*((exp(1i*100*k1))-1);
D = (exp(1i*150*k2)-exp(1i*50*k2))./(1i*k2);
E = (1./(1-(((100*k2)./pi).^2))).*(1i*k2).*((100./pi).^2).*(exp(1i*50*k2)+exp(1i*150*k2));
F = (4.*sin(50*k2))./k2;
G = (exp(-1i*50*k2)-exp(-1i*150*k2))./(1i*k2);
H = (1./(1-(((100*k2)./pi).^2))).*(1i*k2).*((100./pi).^2).*(exp(-1i*150*k2)+exp(-1i*50*k2));
K1 = 0:100;
K2 = -150:150;
for x1 = 1:101;
for x2 = 1:301;
k1 = K1(x1);
k2 = K2(x2);
if (k==0)
eta_(x1,x2)= 200*eta0*v*t;
else
eta_(x1,x2)= (A.*(B-C)).*((D-E)+F+(G+H));
end
end
end
Zsym = subs(eta_);
Z = double(Zsym);
ETA = ifft2(Z)
mesh(abs(ETA))
figure(2)
mesh(real(ETA))
figure(3)
mesh(imag(ETA))
% M = 101;
% N = 301;
% Zsym = subs(eta_);
% Z = double(Zsym);
% ETA = abs(ifft2(Z));
% I = -M/2:1:M/2-1;
% J = -N/2:1:N/2-1;
% mesh(I,J,ETA)
% % X = abs(ifft2(Z,K1,K2));
% % X = ifftshift(X);
% % % F = [-K1/2:K1/2-1]/K1;
% % plot(F,X)
in this journal fig 8

Answers (0)

Categories

Find more on Mathematics in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!