Layered Damped Soil on Elastic Rock - Site Amplification in Kramer Book

Hi,

I would like to calculate the amplification factor in multilayer damped soil on elastic rock, which is also included in Kramer's book geotechnical earthquake engineering (page 268). After making the inputs of the rock and the soil as contained in the code, I need to extract the amplification that I expected to see on the surface, but I can't find where I made a mistake. While getting the correct result in another licensed software, matlab gives very high values, especially the first two waves. Is there anyone who is working on this issue who can help me ? Thank you

Kramer Book

Code:

% Define soil layer properties
H = [20, 30, 40];          % Thickness of the layers (m)
Vs = [200, 300, 400];      % Shear wave velocities of the layers (m/s)
rho = [1800, 1900, 2000];  % Densities of the layers (kg/m^3)
xi = [0.05, 0.04, 0.03];   % Damping ratios of the layers
% Rock properties
Vs_rock = 800;             % Shear wave velocity of the rock (m/s)
rho_rock = 2200;           % Density of the rock (kg/m^3)
xi_rock = 0.02;            % Damping ratio of the rock
% Define frequency range
f = linspace(0.1, 20, 1000);  % Frequency (Hz)
omega = 2 * pi * f;          % Angular frequency (rad/s)
% Create an empty array to store amplification values
A = zeros(size(f));
% Calculate amplification for each frequency
for j = 1:length(f)
    % Complex wave number for rock
    k_rock = (omega(j) / Vs_rock) * (1 - 1i * xi_rock);
    
    % Reflection at the rock base (initial transfer matrix is identity)
    T_total = eye(2);  % Identity matrix as the starting point
    
    % Calculate wave number and transfer matrix for each soil layer
    for i = 1:length(H)
        % Complex wave number for each soil layer (using updated formula)
        k = (omega(j) / Vs(i)) * (1 - 1i * xi(i));
        
        % Transfer matrix for the current layer
        Ti = [cos(k * H(i)), sin(k * H(i)) / (rho(i) * Vs(i)); 
            -rho(i) * Vs(i) * sin(k * H(i)), cos(k * H(i))];
        
        % Update the total transfer matrix
        T_total = Ti * T_total;
    end
    
    % Apply free surface boundary condition at the surface
    A(j) = 1 / abs(T_total(1,1));
end
% Plot the frequency-amplification graph
figure;
plot(f, A, 'LineWidth', 2);
xlabel('Frequency (Hz)');
ylabel('Amplification Factor');
title('Layered Soil Amplification over Elastic Rock (with Damping)');
grid on;

Kramer Book Screenshots: