I Need correct or rewrite code

Dear all;

Please I need to correct or rewrite the following code to get the depths displaed as to some extent like figure in the annexed file based on a rule you extracted from data (the Excel file of data is annexed also)

clc;
close all;
clear;

% MATLAB Code for Faulted Layer Analysis

% Constants
G = 6.67430e-3;  % Gravitational constant in m^3 kg^-1 s^-2 (scaled for this application)


%=========================================================================%
% Density for different sediment layers

rho_profile_1 = 2.430;
rho_profile_2 = 2.580;
rho_profile_3 = 2.780;
rho_profile_4 = 2.670;
rho_profile_5 = 2.800;

rho_sediments = [rho_profile_1, rho_profile_2, rho_profile_3, rho_profile_4, rho_profile_5];
N = length(rho_sediments);  % Number of sediment layers
rho_average = sum(rho_sediments) / N;  % Average density of sediments (g/cm^3)

rho_profile_6 = 2.840; % Density of basement rock (g/cm^3)


% Density contrasts for different sediment layers

rho_cont_profile_1 = rho_profile_1 - rho_profile_6;
rho_cont_profile_2 = rho_profile_2 - rho_profile_6;
rho_cont_profile_3 = rho_profile_3 - rho_profile_6;
rho_cont_profile_4 = rho_profile_4 - rho_profile_6;
rho_cont_profile_5 = rho_profile_5 - rho_profile_6;
rho_cont_profile_6 = rho_average - rho_profile_6;

% Load Bouguer gravity anomaly profile from Excel file
% Load the Excel data
filename = 'Kifl oil field_Iraq.xlsx';    % Profile of Bouguer KIFL OIL.xlsx            Kifl oil field_Iraq
data = xlsread(filename, 'Sheet1');

% Extract relevant columns
x = data(:, 1);            % Distance (x in km)
g = data(:, 2);    % Bouguer gravity anomaly (g in mGal)

%=========================================================================%

% Model
g_profile_1 = 2 * pi * G * rho_cont_profile_1 * 0.6;
g_profile_2 = 2 * pi * G * rho_cont_profile_2 * 1.0;
g_profile_3 = 2 * pi * G * rho_cont_profile_3 * 1.950;
g_profile_4 = 2 * pi * G * rho_cont_profile_4 * 2.950;
g_profile_5 = 2 * pi * G * rho_cont_profile_5 * 4.00;

g_profile_6 = 2 * pi * G * rho_cont_profile_6 * 8.00;

% Depths
z1 = abs(g-g_profile_1 / (2 * pi * G * rho_cont_profile_1));
z2 = abs(g-g_profile_2 / (2 * pi * G * rho_cont_profile_2));
z3 = abs(g-g_profile_3 / (2 * pi * G * rho_cont_profile_3));
z4 = abs(g-g_profile_4 / (2 * pi * G * rho_cont_profile_4));
z5 = abs(g-g_profile_5 / (2 * pi * G * rho_cont_profile_5));

z6 = abs(g-g_profile_6 / (2 * pi * G * rho_cont_profile_6));

% Collect depth profiles in a cell array
z_profiles = {z1, z2, z3, z4, z5, z6};
%=========================================================================%

% Apply the logical condition to eliminate depth increases
for i = 1:length(z_profiles)
    for j = 2:length(z_profiles{i})
        % Check if current depth exceeds the previous depth
        if z_profiles{i}(j) < z_profiles{i}(j-1)
            %z_profiles{i}(j) > z_profiles{i}(j-1)
            % Set depth to previous depth value to eliminate increase
            
            z_profiles{i}(j) = z_profiles{i}(j-1);
            
        end
    end
end

%=========================================================================%

% Calculate dip angles and fault types
location_km = x(2:end);  % Exclude first point to match dip angle array size
dip_angle_deg = zeros(length(location_km), 1); % Pre-allocate for dip angles
fault_type = cell(length(location_km), 1);    % Pre-allocate for fault types

for i = 2:length(x)
    x_fault = x(i) - x(i-1);
    delta_z = z_profiles{end}(i) - z_profiles{end}(i-1);
    dip_angle_deg(i-1) = atand(delta_z / x_fault);
    
    % Assign fault type based on dip angle
    if dip_angle_deg(i-1) == 90
        fault_type{i-1} = 'Vertical';
    elseif dip_angle_deg(i-1) > 30 && dip_angle_deg(i-1) < 70
        fault_type{i-1} = 'Normal';
    elseif dip_angle_deg(i-1) >= 70
        fault_type{i-1} = 'Reverse';
    else
        fault_type{i-1} = 'Normal';
    end
end

% Create fault analysis table
fault_analysis_table = table(location_km, fault_type, dip_angle_deg, ...
    'VariableNames', {'Location_km', 'Fault_Type', 'Dip_Angle_deg'});

% Display the filtered results for specific locations
locations_to_display = [8.5; 14; 21.5];
filtered_results = fault_analysis_table(ismember(fault_analysis_table.Location_km, locations_to_display), :);
if ~isempty(filtered_results)
    disp('Fault Analysis Results for specified locations:');
    disp(filtered_results);
end

% Plot Bouguer anomaly profile
subplot(2, 1, 1);
plot(x, g, 'b', 'LineWidth', 1);
title('Bouguer Anomaly');
xlabel('Distance (km)');
ylabel('Gravity (mGal)');
grid on;

% Plot inverted depth profiles and find intersections with depth = 0
subplot(2, 1, 2);
hold on;
colors = {'r', 'g', 'b', 'y', 'm'};

for i = 1:length(z_profiles) - 1
    plot(x, z_profiles{i}, 'Color', colors{i}, 'LineWidth', 1);
end

% Plot the final profile (basement depth) with black color
plot(x, z_profiles{end}, 'k', 'LineWidth', 1);
title('Inverted Depth Profiles for Faulted Layers');
xlabel('Distance (km)');
ylabel('Depth (km)');
set(gca, 'YDir', 'reverse');
grid on;

% Define the horizontal horizon value
horizon_value = 0;

% Plot the horizontal horizon line
plot(x, horizon_value * ones(size(x)), 'k--', 'LineWidth', 1);

% Define the horizontal horizon value
horizon_value = 0;

% Plot the horizontal horizon line
plot(x, horizon_value * ones(size(x)), 'k--', 'LineWidth', 1);

% Find and plot intersections with the horizontal horizon
for i = 1:length(z_profiles) - 1
    z_profile = z_profiles{i};
    crossings = find(diff(sign(z_profile - horizon_value)) ~= 0);
    
    for j = 1:length(crossings)
        idx = crossings(j);
        x_cross = x(idx) + (horizon_value - z_profile(idx)) * (x(idx+1) - x(idx)) / (z_profile(idx+1) - z_profile(idx));
        z_cross = horizon_value;
        
        plot([x_cross, x_cross], [min(z_profiles{end}), max(z_profiles{end})], 'k--', 'LineWidth', 1.5);
        
        dip_angle_label = sprintf('%.2f°', dip_angle_deg(idx));
        text(x_cross, z_cross, dip_angle_label, 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', ...
            'Color', 'red', 'FontSize', 8);
    end
end

hold off;
legend('Profile1', 'Profile2', 'Profile3', 'Profile4', 'Profile5', 'depth_basement', 'Location', 'Best');

% Plot all profiles in figure 2
figure(2)
hold on;
colors = {'r', 'g', 'b', 'y', 'm', 'k'};
for i = 1:length(z_profiles)
    plot(x, z_profiles{i}, 'Color', colors{i}, 'LineWidth', 1);
end
title('Inverted Depth Profiles for Faulted Layers');
xlabel('Distance (km)');
ylabel('Depth (km)');
set(gca, 'YDir', 'reverse');
grid on;

[dpb]

At least format your code and show us what you've gotten and where/what you think isn't what expect/want...

It would seem the figure could be created pretty easily with area given the contours versus location values.

[Moustafa]

We an Excel file and the mainly first two columns are of a digitized data x, g, where x is distance in km and g Bouguer gravity anomaly in mGal respectively .The calculated (inverted depth) z is estimated using the Bouguer gravity anomaly, expressed mathematically as: z = g / (2 * pi * G * rho_contrast)

the depths calcualted for all layers sould be interpretted to give the results similar figure as shown in annexed file

[Steven Lord]

What specific help are you looking for? "Rewrite my code" is too broad a request for MATLAB Answers IMO.

[Moustafa]

See please my annexed files with my trail

[Steven Lord]

Please don't post follow-up comments with additional information as answers. Reserve messages in the Answers section of this page for things that are suggestions that may resolve the problem. Post additional information as comments instead.

[Walter Roberson]

"See please my annexed files with my trail"

Unfortunately, that does not give us any more information about what specific help you are hoping we will provide.