debug my matlab program:
% Borehole Acoustic Wave Propagation Simulation
% Calculates dispersion curves and synthetic seismograms for a fluid-filled borehole
% in an elastic formation using the Newton-Raphson method and a 4x4 matrix.

% Clear workspace and figures
clear all; close all; clc;

% Parameters from the text
% Choose either fast or slow formation
formation_choice = 'fast'; % Options: 'fast' or 'slow'

switch formation_choice
    case 'fast'
        % Fast Formation
        Vp = 4000;    % P-wave velocity in formation (m/s)
        Vs = 2300;    % S-wave velocity in formation (m/s)
        rho_s = 2500; % Density of the formation (kg/m^3)
    case 'slow'
        % Slow Formation
        Vp = 3000;    % P-wave velocity in formation (m/s)
        Vs = 1354;    % S-wave velocity in formation (m/s)
        rho_s = 2500; % Density of the formation (kg/m^3)
    otherwise
        error('Invalid formation choice.');
end

% Fluid properties
Vf = 1500;    % P-wave velocity in fluid (m/s)
rho_f = 1000; % Density of the fluid (kg/m^3)

% Borehole radius
R = 0.1;      % Borehole radius (m)

% Frequency range up to 25 kHz
fmin = 100;    % Minimum frequency (Hz)
fmax = 25000;  % Maximum frequency (Hz)
df = 100;      % Frequency step (Hz)
freq = fmin:df:fmax; % Frequency vector

% Initialize variables for dispersion curves
Nfreq = length(freq);
k_mode_monopole = zeros(1,Nfreq);
k_mode_dipole = zeros(1,Nfreq);

% Tolerance and maximum iterations for Newton-Raphson method
tol = 1e-6;
max_iter = 100;

% Loop over frequencies to compute dispersion curves
disp('Calculating dispersion curves...');
for n = 1:Nfreq
    omega = 2*pi*freq(n);

    % Initial guess for wavenumber (monopole and dipole)
    k_guess_mono = omega / Vp; % Initial guess for monopole mode
    k_guess_dipo = omega / Vs; % Initial guess for dipole mode

    % Solve dispersion equation for monopole mode (n=0)
    [k_mono, success_mono] = solve_dispersion(omega, 0, k_guess_mono, Vp, Vs, Vf, rho_s, rho_f, R, tol, max_iter);
    if success_mono
        k_mode_monopole(n) = k_mono;
    else
        k_mode_monopole(n) = NaN;
    end

    % Solve dispersion equation for dipole mode (n=1)
    [k_dipo, success_dipo] = solve_dispersion(omega, 1, k_guess_dipo, Vp, Vs, Vf, rho_s, rho_f, R, tol, max_iter);
    if success_dipo
        k_mode_dipole(n) = k_dipo;
    else
        k_mode_dipole(n) = NaN;
    end
end

% Compute phase velocities
vph_mono = omega ./ real(k_mode_monopole);
vph_dipo = omega ./ real(k_mode_dipole);

% Plot dispersion curves
figure;
plot(freq/1000, vph_mono, 'b-o', 'LineWidth', 2, 'DisplayName', 'Monopole Mode (n=0)'); hold on;
plot(freq/1000, vph_dipo, 'r-s', 'LineWidth', 2, 'DisplayName', 'Dipole Mode (n=1)');
xlabel('Frequency (kHz)');
ylabel('Phase Velocity (m/s)');
title(['Dispersion Curves - ' upper(formation_choice) ' Formation']);
legend('Location', 'best');
grid on;

% Generate synthetic seismograms
disp('Calculating synthetic seismograms...');
% Time axis
dt = 1e-5;          % Time step (s)
tmax = 0.01;        % Maximum time (s)
t = 0:dt:tmax;      % Time vector
Nt = length(t);

% Source spectrum (Ricker wavelet)
f0 = 8000; % Central frequency (Hz)
S_omega = @(omega) (omega.^2 / (2*pi*f0)^2) .* exp(-(omega / (2*pi*f0)).^2);

% Receiver positions along z-axis
zmin = 0.5;   % Minimum receiver depth (m)
zmax = 2;     % Maximum receiver depth (m)
dz = 0.01;    % Depth step (m)
z = zmin:dz:zmax;
Nz = length(z);

% Initialize seismograms
seismogram_mono = zeros(Nt, Nz);
seismogram_dipo = zeros(Nt, Nz);

% Frequency vector for IFFT (including negative frequencies)
Nfft = 2^nextpow2(Nt);
freq_fft = (0:Nfft/2-1)*(1/(dt*Nfft));

% Loop over receiver positions to compute seismograms
for iz = 1:Nz
    % Distance from source
    z_rec = z(iz);

    % Frequency domain representation
    P_mono = zeros(1, Nfft/2);
    P_dipo = zeros(1, Nfft/2);

    for nfreq = 1:length(freq_fft)
        omega = 2*pi*freq_fft(nfreq);
        freq_current = freq_fft(nfreq);

        % Skip frequencies outside the calculated range
        if freq_current < fmin || freq_current > fmax
            continue;
        end

        % Find index in freq vector
        idx = find(abs(freq - freq_current) < df/2);

        if isempty(idx)
            continue;
        end

        k_mono = k_mode_monopole(idx);
        k_dipo = k_mode_dipole(idx);

        % Skip if k is NaN
        if isnan(k_mono) || isnan(k_dipo)
            continue;
        end

        % Compute source spectrum
        S = S_omega(omega);

        % Monopole mode contribution
        P_mono(nfreq) = S * exp(-1i * k_mono * z_rec);

        % Dipole mode contribution
        P_dipo(nfreq) = S * exp(-1i * k_dipo * z_rec);
    end

    % Construct full spectrum for IFFT
    P_mono_full = [P_mono, conj(fliplr(P_mono(2:end)))];
    P_dipo_full = [P_dipo, conj(fliplr(P_dipo(2:end)))];

    % Inverse Fourier Transform to time domain
    time_signal_mono = ifft(P_mono_full, 'symmetric');
    time_signal_dipo = ifft(P_dipo_full, 'symmetric');

    % Truncate or interpolate to desired time axis
    seismogram_mono(:, iz) = interp1((0:Nfft-1)*dt, time_signal_mono, t, 'linear', 0);
    seismogram_dipo(:, iz) = interp1((0:Nfft-1)*dt, time_signal_dipo, t, 'linear', 0);
end

% Plot seismograms
figure;
subplot(2,1,1);
imagesc(z, t*1e3, seismogram_mono);
xlabel('Depth (m)');
ylabel('Time (ms)');
title(['Synthetic Seismogram - Monopole Source (' upper(formation_choice) ' Formation)']);
colorbar;
set(gca, 'YDir', 'normal');

subplot(2,1,2);
imagesc(z, t*1e3, seismogram_dipo);
xlabel('Depth (m)');
ylabel('Time (ms)');
title(['Synthetic Seismogram - Dipole Source (' upper(formation_choice) ' Formation)']);
colorbar;
set(gca, 'YDir', 'normal');

disp('Simulation complete.');

% --- Function Definitions ---

% Function to solve dispersion equation using Newton-Raphson method
function [k_sol, success] = solve_dispersion(omega, n, k_guess, Vp, Vs, Vf, rho_s, rho_f, R, tol, max_iter)
    k = k_guess;
    success = false;
    for iter = 1:max_iter
        % Compute dispersion function and its derivative
        [Dk, Dk_deriv] = dispersion_function(k, omega, n, Vp, Vs, Vf, rho_s, rho_f, R);

        % Avoid division by zero
        if Dk_deriv == 0
            break;
        end

        % Update k using Newton-Raphson
        k_new = k - Dk / Dk_deriv;

        % Check for convergence
        if abs(k_new - k) < tol
            k_sol = k_new;
            success = true;
            return;
        end

        k = k_new;
    end
    % If not converged, return NaN
    k_sol = NaN;
end

% Function to compute dispersion function D(k, omega) and its derivative
function [Dk, Dk_deriv] = dispersion_function(k, omega, n, Vp, Vs, Vf, rho_s, rho_f, R)
    % Compute wavenumbers
    kp = omega / Vp;
    ks = omega / Vs;
    kf = omega / Vf;

    % Radial wavenumbers
    krp = sqrt(kp^2 - k^2);
    krs = sqrt(ks^2 - k^2);
    krf = sqrt(k^2 - kf^2);

    % Choose correct branches for square roots
    if imag(krp) < 0
        krp = -krp;
    end
    if imag(krs) < 0
        krs = -krs;
    end
    if real(krf) < 0
        krf = -krf;
    end

    % Bessel functions and their derivatives
    % Solid (Formation)
    Jn_krp_R = besselj(n, krp * R);
    Yn_krp_R = bessely(n, krp * R);
    Jn_krs_R = besselj(n, krs * R);
    Yn_krs_R = bessely(n, krs * R);

    % Derivatives using standard formulas
    Jn_prime_krp_R = -krp * R * besselj(n+1, krp * R) + n * Jn_krp_R;
    Yn_prime_krp_R = -krp * R * bessely(n+1, krp * R) + n * Yn_krp_R;
    Jn_prime_krs_R = -krs * R * besselj(n+1, krs * R) + n * Jn_krs_R;
    Yn_prime_krs_R = -krs * R * bessely(n+1, krs * R) + n * Yn_krs_R;

    % Fluid (Borehole)
    Jn_krf_R = besselj(n, krf * R);
    Jn_prime_krf_R = krf * R * besselj(n-1, krf * R) - n * Jn_krf_R;

    % Assemble the 4x4 matrix M
    M = zeros(4,4);

    % Radial displacement continuity at r = R
    M(1,1) = k_{rp} * Jn_prime_krp_R;
    M(1,2) = k_{rp} * Yn_prime_krp_R;
    M(1,3) = -\frac{n}{R} * Jn_krs_R;
    M(1,4) = -\frac{n}{R} * Yn_krs_R;

    % Radial stress continuity at r = R
    lambda = rho_s * (Vp^2 - 2 * Vs^2);
    mu = rho_s * Vs^2;
    kappa = lambda + 2 * mu; % Bulk modulus

    % Radial stress components
    % From compressional potential
    sigma_rr_phi_A = (2 * mu * krp^2) * Jn_krp_R;
    sigma_rr_phi_B = (2 * mu * krp^2) * Yn_krp_R;

    % From shear potential
    sigma_rr_psi_C = -2 * mu * k * krs * Jn_krs_R;
    sigma_rr_psi_D = -2 * mu * k * krs * Yn_krs_R;

    M(2,1) = sigma_rr_phi_A;
    M(2,2) = sigma_rr_phi_B;
    M(2,3) = sigma_rr_psi_C;
    M(2,4) = sigma_rr_psi_D;

    % Shear stress sigma_rtheta must be zero at r = R
    % Shear stress from shear potential
    sigma_rtheta_psi_C = mu * ( (n^2 / (R^2)) - k^2 ) * Jn_krs_R;
    sigma_rtheta_psi_D = mu * ( (n^2 / (R^2)) - k^2 ) * Yn_krs_R;

    M(3,1) = 0;
    M(3,2) = 0;
    M(3,3) = sigma_rtheta_psi_C;
    M(3,4) = sigma_rtheta_psi_D;

    % Shear stress sigma_rz must be zero at r = R
    % Shear stress from compressional potential
    sigma_rz_phi_A = -2 * mu * k * krp * Jn_prime_krp_R;
    sigma_rz_phi_B = -2 * mu * k * krp * Yn_prime_krp_R;

    % Shear stress from shear potential
    sigma_rz_psi_C = -2 * mu * k * krs * Jn_prime_krs_R;
    sigma_rz_psi_D = -2 * mu * k * krs * Yn_prime_krs_R;

    M(4,1) = sigma_rz_phi_A;
    M(4,2) = sigma_rz_phi_B;
    M(4,3) = sigma_rz_psi_C;
    M(4,4) = sigma_rz_psi_D;

    % Compute determinant Dk
    Dk = det(M);

    % Compute derivative of Dk with respect to k numerically
    delta_k = 1e-6;
    k_plus = k + delta_k;

    % Recompute Dk at k_plus
    [Dk_plus] = dispersion_function_no_derivative(k_plus, omega, n, Vp, Vs, Vf, rho_s, rho_f, R);

    % Compute derivative
    Dk_deriv = (Dk_plus - Dk) / delta_k;
end

% Function to compute Dk without derivative
function [Dk] = dispersion_function_no_derivative(k, omega, n, Vp, Vs, Vf, rho_s, rho_f, R)
    % Similar to dispersion_function but only computes Dk
    % Code omitted for brevity (same as dispersion_function without derivative calculations)
    % ...
    % Return Dk
    Dk = det(M);
end