Hi guys I want to write a code to calculate G11, for each station in range frequencies 0 to 15 and plot im(G11) against frequencies. I write the below code but it doesn't work

clc
clear 

alpha = 2;  % Shear wave velocity
beta = 1;  % Compressional wave velocity
stations = 0:0.01:0.1;  % Range of stations
frequencies = 0:15;  % Range of frequencies
epsilon = 1e-10;  % Small epsilon value
ro=1; %mass density

for i = 0:length (stations)
    for j = 0:length (frequencies)
        omega= 2 * pi * frequencies(j);
        q=omega/alpha;
        k=omega/beta;
        f1=(beta^2 / alpha^2) * (1 - 1i * 2 / (q * stations(i) + epsilon) - 2 / (q^2 * stations(i)^2 + epsilon)) * exp(-1i * q * stations(i))+ [1i*2/(k* stations(i)+epsilon)+2/(k^2*stations(i)^2+epsilon)]*exp(-1i*k*stations(i));
        f2=(beta^2 / alpha^2) * (1i / (q * stations(i) + epsilon) + 1 / (q^2 * stations(i)^2 + epsilon)) * exp(-1i * q * stations(i))+ [1-1i/(k*stations(i)+epsilon)-1/(k^2*stations(i)^2+epsilon)]*exp(-1i*k*stations(i));
        G11= 1/(4*pi*ro*(stations(i)+epsilon))*(f1);
        
    end
end

        
plot(frequencies, G11);

Answer 1

[Sulaymon Eshkabilov]

There are a couple of mistakes:

(1) Indices can be "0"- see loop i = 0:length(stations), j=0:length(frequencies)

(2) The last value of G11 is saved nd not all values

(3) G11 is complex values and thus, you can plot real, imaginary or absolute values of G11

(4) You need to check your formulations for f1, f2 and G11 !!!

Here is the corrected code:

alpha = 2;  % Shear wave velocity
beta = 1;  % Compressional wave velocity
stations = 0:0.01:0.1;  % Range of stations
frequencies = 0:15;  % Range of frequencies
epsilon = 1e-10;  % Small epsilon value
ro=1; %mass density

for i = 1:length(stations)
    for j = 1:length(frequencies)
        omega= 2 * pi * frequencies(j);
        q=omega/alpha;
        k=omega/beta;
        f1=(beta^2 / alpha^2) * (1 - 1i * 2 / (q * stations(i) + epsilon) - 2 / (q^2 * stations(i)^2 + epsilon)) * exp(-1i * q * stations(i))+ [1i*2/(k* stations(i)+epsilon)+2/(k^2*stations(i)^2+epsilon)]*exp(-1i*k*stations(i));
        f2=(beta^2 / alpha^2) * (1i / (q * stations(i) + epsilon) + 1 / (q^2 * stations(i)^2 + epsilon)) * exp(-1i * q * stations(i))+ [1-1i/(k*stations(i)+epsilon)-1/(k^2*stations(i)^2+epsilon)]*exp(-1i*k*stations(i));
        G11(i,j)= 1/(4*pi*ro*(stations(i)+epsilon))*(f1);   
    end
    plot(frequencies, abs(G11(i,:)))  % Absolute vaues are computed
    % plot(frequencies, real(G11(i,:))) % Real values are computed
    % plot(frequencies, real(G11(i,:))) % Imaginary values are computed
    hold all
end