How to plot for this code?
2 views (last 30 days)
Show older comments
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));
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
if (k1 == 0)
eta_(x1,x2)= (A.*(100-C)).*((D-E)+F+(G-H));
else
eta_(x1,x2)= (A.*(B-C)).*((D-E)+F+(G-H));
end
if (k2 == 0)
eta_(x1,x2)= (A.*(B-C)).*((100-E)+200+(100-H));
else
eta_(x1,x2)= (A.*(B-C)).*((D-E)+F+(G-H));
end
end
end
end
ETA = ifft2(eta_)
X = abs(ifft2(ETA,K1,K2));
X = ifftshift(X);
F = [-K1/2:K1/2-1]/K1;
plot(F,X)
0 Comments
Answers (1)
Image Analyst
on 21 Feb 2021
Since youi're calline fft2(), you're dealing with images not 1-D signals. So call imshow() instead of plot(), unless you want to display some cross sections (profiles) through the image.
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!