Hi @MANYINGA MARC, I understand that you need help in plotting of reflection and transmission coefficient and to plot a dispersion fucntion as a function of frequency or incident angle by using SMM.
I have done some modifications in the code file you provided. Below is the code after modification
%% Step 1: Define the Units %%%
Pa = 1;
GPa = 1e9 * Pa;
meters = 1;
micrometers = 1e-6 * meters;
mm = 1e-3 * meters;
deg = pi / 180;
cf = 1400; % speed of sound in fluid
%%% Parameters %%%
n = 30;
H = 3.434 * mm;
h = H / 30;
rho = 1500;
rho_f = 1000;
c11 = 12.1 * GPa; c12 = 5.5 * GPa; c44 = 6.95 * GPa; c22 = 12.3 * GPa;
c13 = 5.9 * GPa; c55 = 6.21 * GPa; c33 = 132 * GPa; c23 = 6.9 * GPa; c66 = 3.32 * GPa;
n11 = 0.043 * GPa; n22 = 0.037 * GPa; n33 = 0.4 * GPa; n12 = 0.021 * GPa;
n23 = 0.001 * GPa; n13 = 0.016 * GPa; n44 = 0.02 * GPa; n55 = 0.015 * GPa; n66 = 0.009 * GPa;
f = 2.242 * 1e6;
w = 2 * pi * f;
k = w / cf;
theta = (-90:1:90) * deg; % Convert degrees to radians
Tx = zeros(size(theta));
Rx = zeros(size(theta));
Dx = zeros(size(theta));
c = cos(theta);
s = sin(theta);
for q = 1:length(theta)
T = diag([c55 - 1i * w * n55, c44 - 1i * w * n44, c33 - 1i * w * n33]);
R = zeros(3, 3);
R(1,3) = (-c13 + 1i * w * n13) * c(q);
R(2,3) = (-c23 + 1i * w * n23) * c(q);
R(3,1) = (c55 - 1i * w * n55) * c(q);
R(3,2) = (c44 - 1i * w * n44) * s(q);
Q = zeros(3, 3);
Q(1,1) = (c11 - 1i * w * n11) * (c(q))^2 - rho * (c(q))^2 + (c66 - 1i * w * n66) * (s(q))^2;
Q(1,2) = (c12 + c66 - 1i * w * (n12 + n66)) * s(q) * c(q);
Q(2,1) = (c12 + c66 - 1i * w * (n12 + n66)) * s(q) * c(q);
Q(2,2) = (c66 - 1i * w * n66) * (c(q))^2 - rho * (c(q))^2 + (c22 - 1i * w * n22) * (s(q))^2;
Q(3,3) = (-c55 + 1i * w * n55) * (c(q))^2 - rho * (c(q))^2 + (-c44 + 1i * w * n44) * (s(q))^2;
N = [(T^-1) * R', T^-1; -Q - R * (T^-1) * R', -R * (T^-1)];
[V, D] = eig(N);
A1 = V(1:3, 1:3);
A2 = V(1:3, 4:6);
B1 = V(4:6, 1:3);
B2 = V(4:6, 4:6);
EXP = diag([exp(D(1,1) * k * h), exp(D(6,6) * k * h), exp(D(4,4) * k * h)]);
KJ = [B2, B1 * EXP; B2 * EXP, B1] * ([A2, A1 * EXP; A2 * EXP, A1])^-1;
KJ11 = KJ(1:3, 1:3);
KJ12 = KJ(1:3, 4:6);
KJ21 = KJ(4:6, 1:3);
KJ22 = KJ(4:6, 4:6);
% Redefine KB11, KB12, KB21, KB22 for the iterative process
KB11 = KJ11;
KB12 = KJ12;
KB21 = KJ21;
KB22 = KJ22;
for ii = 1:10
KA11 = KJ11;
KA12 = KJ12;
KA21 = KJ21;
KA22 = KJ22;
kg = [KA11 + KA12 * ((KB11 - KA22)^-1) * KA21, -KA12 * ((KB11 - KA22)^-1) * KB12;
KB21 * ((KB11 - KA22)^-1) * KA21, KB22 - KB21 * ((KB11 - KA22)^-1) * KB12];
KJ11 = kg(1:3, 1:3);
KJ12 = kg(1:3, 4:6);
KJ21 = kg(4:6, 1:3);
KJ22 = kg(4:6, 4:6);
end
K = [KJ11, KJ12; KJ21, KJ22];
S = K^-1;
lamda = c(q) / (1i * w * rho_f * cf);
T = (2 * lamda * S(6,3)) / ((S(3,3) + lamda) * (S(6,6) - lamda) - S(6,3) * S(3,6));
Tx(q) = T;
R = (-(S(1,1) - lamda) * (S(2,2) - lamda) + S(1,2) * S(2,1)) / ((S(2,2) - lamda) * (S(1,1) + lamda) - S(1,2) * S(2,1));
Rx(q) = R;
D = ((S(3,3) + lamda) * (S(6,6) - lamda) - S(6,3) * S(3,6));
Dx(q) = D;
end
figure(1)
plot(theta / deg, abs(Tx), 'r', 'linewidth', 1)
xlabel('\theta (degrees)', 'fontsize', 14)
ylabel('T (Transmission)', 'fontsize', 14)
title('Transmission Coefficient (T)', 'fontsize', 14)
set(gcf, 'color', 'white');
figure(2)
plot(theta / deg, abs(Rx), 'b', 'linewidth', 1)
xlabel('\theta (degrees)', 'fontsize', 14)
ylabel('R (Reflection)', 'fontsize', 14)
title('Reflection Coefficient (R)', 'fontsize', 14)
set(gcf, 'color', 'white');
figure(3)
plot(theta / deg, abs(Dx), 'g', 'linewidth', 1)
xlabel('\theta (degrees)', 'fontsize', 14)
ylabel('D (Dispersion)', 'fontsize', 14)
title('Dispersion Function (D)', 'fontsize', 14)
set(gcf, 'color', 'white');
Below is the screenshot of the output:
