Error in MATLAB Code : BD18_Inverse Function

I am currently working on a MATLAB project involving the BD18_Inverse function for geotechnical data analysis, specifically the inverse filtering procedure. I am facing an issue with my MATLAB code and would appreciate your expertise in troubleshooting the problem.

The code executes the BD18_Inverse function, but I encounter an error that prevents the code from running successfully. The error message I receive is as follows:

```

Array indices must be positive integers or logical values.

Error in BD18_Inverse/CalculateW1 (line 79)

z50 = 1 + 2 * (C2 .* (z50_ref - 1)) .* (1 - 1 ./ (1 + ((q0t) ./ qinv(abs(z_prime)))).^0.5);

Error in BD18_Inverse (line 12)

w1 = CalculateW1(z_prime, n_m, dcone);

Error in BD18_Script (line 13)

BD18_Inverse(z, conetip, fsm, dcone);

```

Complete Code is here

function [qinv, fsinv] = BD18_Inverse(z, qm, fsm, dcone, sigma_v, Pa)
    % Initialize parameters
    n_m = 1;
    z_prime = (z - z(1)) / dcone;
    mt = 0.1;

    % Initialize qinv as qm
    qinv = qm;
    
    % Part 1: Inversion for Tip Resistance
    while true
        w1 = CalculateW1(z_prime, n_m, dcone);
        w2 = CalculateW2(qinv, n_m);
        wc = (w1 .* w2) / sum(w1 .* w2);
        
        % Calculate qtz' for smoothing
        qtz0 = mean(qinv(abs(z_prime) <= 0.05));

        for i = 1:length(z_prime)
            qtz0 = mean(qinv(abs(z_prime - z_prime(i)) <= 0.05));
            qz0t = sum(qinv .* wc);
            qinv(i) = qm(i) + (qinv(i) - qz0t) + (qtz0 - qz0t);

        end
        
        % Smoothing
        qinv = smooth(qinv, 5);
        
        % Calculate err
        err = sum(abs(qm - qinv)) / sum(qm);
        
        if err < 1e-6
            break;
        end
    end
    
    % Smoothing to remove high-frequency noise
    qinv = smooth(qinv, 5);
    
    % Part 2: Inversion for Sleeve Friction
    for i = 1:length(z)
        Q_tn_m = CalculateQt(qm, sigma_v, Pa, n_m);
        F_r_m = CalculateFr(fsm, qm, sigma_v);
        I_c_m = CalculateIc(Q_tn_m, F_r_m);
        n_m = CalculateN(I_c_m, sigma_v, Pa);

        % Check for convergence (n_m does not change significantly)
        if abs(n_m - 1) < 1e-6
            break;
        end
    end

    % Calculate Qtn for final iteration
    Q_tn_m = CalculateQt(qm, sigma_v, Pa, n_m);

    % Part 3: Interface Correction
    sharp_transitions = abs(log(qinv(2:end)) - log(qinv(1:end-1))) / (z_prime(2:end) - z_prime(1:end-1)) > mt;
    transition_zones = IdentifyTransitionZones(sharp_transitions, dcone, z_prime);
    qinv = PerformInterfaceCorrection(qinv, transition_zones);
    
    % Calculate fsinv
    Q_tn_inv = CalculateQt(qinv, sigma_v, Pa, n_m);
    F_r_inv = CalculateFr(fsm, qinv, sigma_v);
    I_c_inv = CalculateIc(Q_tn_inv, F_r_inv);
    n_inv = CalculateN(I_c_inv, sigma_v, Pa);
    
    fsinv = CalculateFsInv(F_r_inv, qinv, sigma_v);

    % Helper functions for Part 2
  function w1 = CalculateW1(z_prime, n_m, dcone)
    disp(size(z_prime)); % Add this line to check the size of z_prime
    disp(size(qinv));   % Add this line to check the size of qinv
    % Modify the code to match the dimensions
    z_prime = reshape(z_prime, size(qinv));
    z50_ref = 4.2;
    C2 = ones(size(z_prime));
    C2(z_prime >= 0) = 0;
    C2(z_prime >= 4) = 0.8;

    q0t = qinv(abs(z_prime) <= 0.05);
    z50 = 1 + (2 * (C2 .* (z50_ref - 1)) ).* (1 - 1 ./ (1 + ((q0t) ./ qinv(abs(z_prime)))).^0.5);

    w1 = 1 ./ (1 + (z_prime ./ z50).^3);
    w1(abs(z_prime) >= 4) = 0.5;
end

    function w2 = CalculateW2(qinv, n_m)
        qz0t = sum(qinv);
        w2 = sqrt(2 ./ (1 + (qz0t ./ qinv).^n_m));
    end

    function Q_tn = CalculateQt(q, sigma_v, Pa, n)
        Q_tn = ((q - sigma_v) ./ Pa) .* (Pa ./ sigma_v).^n;
    end

    function F_r = CalculateFr(fs, q, sigma_v)
        F_r = (fs ./ (q - sigma_v)) * 100;
    end

    function I_c = CalculateIc(Q_tn, F_r)
        I_c = sqrt((3.47 - log(Q_tn)).^2 + (log(F_r) + 1.22).^2);
    end

    function n = CalculateN(I_c, sigma_v, Pa)
        n = 0.381 * I_c + 0.05 * sigma_v / Pa - 0.15;
        n(n > 1) = 1;
    end

    % Helper functions for Part 3
    function transition_zones = IdentifyTransitionZones(sharp_transitions, dcone, z_prime)
        transition_zones = [];
        i = 1;

        while i <= length(sharp_transitions)
            if sharp_transitions(i) == 1
                transition_start = i;

                while i < length(sharp_transitions) && sharp_transitions(i) == 1
                    i = i + 1;
                end

                transition_end = i;
                transition_thickness = z_prime(transition_end) - z_prime(transition_start);

                if transition_thickness >= 3 * dcone
                    transition_zones = [transition_zones; transition_start, transition_end, transition_thickness];
                end
            end

            i = i + 1;
        end
    end

    function q = PerformInterfaceCorrection(q, transition_zones)
        for i = 1:size(transition_zones, 1)
            start_index = transition_zones(i, 1);
            end_index = transition_zones(i, 2);

            if q(end_index) > q(start_index)
                upper_zone = start_index + round(0.4 * transition_zones(i, 3));
                q(start_index:upper_zone) = q(upper_zone);
                
                lower_zone = upper_zone + 1;
                q(lower_zone:end_index) = q(end_index);
            else
                upper_zone = start_index + round(0.6 * transition_zones(i, 3));
                q(start_index:upper_zone) = q(upper_zone);

                lower_zone = upper_zone + 1;
                q(lower_zone:end_index) = q(end_index);
            end
        end
    end
end

I have tried several troubleshooting steps, including checking the variable sizes, confirming inputs, and debugging with smaller data subsets. However, I have not been able to identify the root cause of the issue.

To assist you in helping me, I have posted the MATLAB code in the following MATLAB Answers thread:

I kindly request your guidance and expertise in resolving this issue. Any insights, suggestions, or solutions you can provide would be greatly appreciated.

Thank you for your time and assistance.

Best regards,

Priyam Mishra

Answer 1

[SAI SRUJAN]

Hi Priyam Mishra,

I understand that you are facing an issue involving the "BD18_Inverse" function for geotechnical data analysis, specifically the inverse filtering procedure.

The error message suggests that there is an issue with array indexing in the "CalculateW1" function.

The error is being caused by the line of code provided below.

z50 = 1 + 2 * (C2 .* (z50_ref - 1)) .* (1 - 1 ./ (1 + ((q0t) ./ qinv(abs(z_prime)))).^0.5);

The variable "z_prime" is utilized for indexing the "qinv" vector. It is essential to ensure that the value of "z_prime" adheres to the conditions, which require array indices to be positive integers or logical values. Additionally, it is advised to verify the absence of any potential division by zero error.

You can refer to the following documentation to understand more about "array indexing" in MATLAB,

https://in.mathworks.com/help/matlab/math/array-indexing.html