How to shift lines to their correction positions (I need to correct figure)

I need to correct figure (3) in the following script to be similar in the attached file (with automatic way)

clc; close all; clear;
% Load the Excel data
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename, 'Sheet1');
% Extract relevant columns
x = data(:, 1); % Distance (x in km)
Field = data(:, 2); % Earth filed (Field in unit)
%==============================================================%
% % Input number of layers and densities
% num_blocks = input('Enter the number of blocks: ');
% block_densities = zeros(1, num_blocks);
% for i = 1:num_blocks
%     block_densities(i) = input(['Enter the density of block ', num2str(i), ' (kg/m^3): ']);
% end
% EDIT to run code here
num_blocks = 5;
block_densities = [2.40, 2.45, 2.50, 2.55, 2.60];

%==============================================================%
% Constants
G = 0.00676; 
Lower_density = 2.67; % in kg/m^3
%==============================================================%
% Calculate inverted depth profile for each layer
z_inv = zeros(length(x), num_blocks);
for i = 1:num_blocks
    density_contrast = block_densities(i) - Lower_density;
    if density_contrast ~= 0
        z_inv(:, i) = Field ./ (2 * pi * G * density_contrast);
    else
        z_inv(:, i) = NaN; % Avoid division by zero
    end
end
%==============================================================%
% Compute vertical gradient (VG) of inverted depth (clean)
VG = diff(z_inv(:, 1)) ./ diff(x);
%==============================================================%
% Set fault threshold and find f indices based on d changes
f_threshold = 0.5; % Threshold for identifying significant d changes
f_indices = find(abs(diff(z_inv(:, 1))) > f_threshold);
%==============================================================%
% Initialize f locations and dip arrays
%==============================================================%
f_locations = x(f_indices);  % Automatically determined f locations
f_dip_angles = nan(size(f_indices)); % Placeholder for calculated dip 
% Calculate dip  for each identified f
for i = 1:length(f_indices)
    idx = f_indices(i);
    if idx < length(x)
        f_dip_angles(i) = atand(abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / (x(idx + 1) - x(idx)));
    else
        f_dip_angles(i) = atand(abs(z_inv(idx, 1) - z_inv(idx - 1, 1)) / (x(idx) - x(idx - 1)));
    end
end
%==============================================================%
%                    Displacement of faults 
%==============================================================%
D_faults = zeros(size(f_dip_angles));
for i = 1:length(f_indices)
    idx = f_indices(i);
    dip_angle_rad = deg2rad(f_dip_angles(i)); % Convert dip to radians
    D_faults(i) = abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / sin(dip_angle_rad);
end
% Assign displacement values
D1 = D_faults(1); % NF displacemen
D2 = D_faults(2); % VF displacement
D3 = D_faults(3); % RF displacement
%==============================================================%
%                 Processing Data for Interpretation
%==============================================================%
A = [x Field z_inv];  % New Data Obtained
col_names = {'x', 'Field'};
for i = 1:num_blocks
    col_names{end+1} = ['z', num2str(i)];
end
dataM = array2table(A, 'VariableNames', col_names);
t1 = dataM;
[nr, nc] = size(t1);
t1_bottoms = t1;
for jj = 3:nc
    for ii = 1:nr-1
        if t1_bottoms{ii, jj} ~= t1_bottoms{ii+1, jj}
            t1_bottoms{ii, jj} = NaN;
        end
    end
end
%==============================================================%
%                     Identifying NaN rows
%==============================================================%
nans = isnan(t1_bottoms{:, 3:end});
nan_rows = find(any(nans, 2));
xc = t1_bottoms{nan_rows, 1}; % Corrected x-coordinates
yc = zeros(numel(nan_rows), 1); % y-coordinates for NaN rows
for ii = 1:numel(nan_rows)
    idx = find(~nans(nan_rows(ii), :), 1, 'last');
    if isempty(idx)
        yc(ii) = 0;
    else
        yc(ii) = t1_bottoms{nan_rows(ii), idx+2};
    end
end
%==============================================================%
%                 Plot f Interpretation
%==============================================================%
figure(1)
plot(A(:, 1), A(:, 3:end))
hold on
grid on
set(gca, 'YDir', 'reverse')
xlabel('Distance (km)'); 
ylabel('D (km)'); 
title('Interpretation of profile data model')

%==============================================================%
figure(2)
hold on
plot(t1_bottoms{:, 1}, t1_bottoms{:, 3:end}, 'LineWidth', 1)
set(gca, 'YDir', 'reverse')
box on
grid on
xlabel('Distance (km)'); 
ylabel('Ds (km)'); 
title('New interpretation of profile data')

%==============================================================%
% Plot the interpreted d profiles
figure(3)
hold on
% Plot the interpreted d profiles
plot(t1_bottoms{:, 1}, t1_bottoms{:, 3:end}, 'LineWidth', 1)
yl = get(gca, 'YLim');  % Get Y-axis limits
% Define f locations and corresponding dip 
f_locations = [7.00, 14.00, 23.00];  
d_angles = [58.47, 90.00, -69.79];    
for ii = 1:numel(f_locations)
    % Find the nearest x index for each fault location
    [~, idx] = min(abs(t1_bottoms{:, 1} - f_locations(ii)));
    
    % Get the starting x and y coordinates (fault starts at the surface)
    x_f = t1_bottoms{idx, 1};  
    y_f = yl(1); % Start at the surface (0 km depth)
    
    % Check if the dip angle is 90° (vf)
    if d_angles(ii) == 90
        x_r = [x_f x_f]; % Vertical line
        y_r = [yl(1), yl(2)];    % From surface to de limit
    else
        % Convert dip to slope (m = tan(angle))
        m = tand(d_angles(ii));             
        
        % Define the x range for fault line
        
        % x_r = linspace(x_f - 5, x_f + 5, 100); % Extend 5 km on each side
        x_r = x;
        y_r = y_f - m * (x_r - x_f);  % Line equation  (+ m)
        
        % Clip y_range within the plot limits
        y_r(y_r > yl(2)) = yl(2);
        y_r(y_r < yl(1)) = yl(1);
        
    end
    
    % Plot the fault lines in black (matching the image)
    plot(x_r, y_r, 'k', 'LineWidth', 3)
    
    % Display dip angles as text near the faults
    %  text(x_f, y_f + 1, sprintf('\\theta = %.2f°', d_angles(ii)), ...
    %    'Color', 'k', 'FontSize', 10, 'FontWeight', 'bold', 'HorizontalAlignment', 'right')
end
set(gca, 'YDir', 'reverse')  % Reverse Y-axis for depth representation
box on
grid on
xlabel('Distance (km)'); 
ylabel('D (km)'); 
title('New Interpretation of Profile Data')

●

%==============================================================%
% Plotting results
figure;
subplot(3,1,1);
plot(x, Field, 'r', 'LineWidth', 2);
title('Field Profile');
xlabel('Distance (km)');
ylabel('Field (unit)');
grid on;
subplot(3,1,2);
plot(x(1:end-1), VG, 'b', 'LineWidth', 1.5);
xlabel('Distance (km)');
ylabel('VG (munit/km)');
title('VG Gradient');
grid on;
subplot(3,1,3);
hold on;
for i = 1:num_blocks
    plot(x, z_inv(:, i), 'LineWidth', 2);
end
title('Inverted D Profile for Each Block');
xlabel('Distance (km)');
ylabel('D (km)');
set(gca, 'YDir', 'reverse');
grid on;

%==============================================================%
% Display results 
fprintf('F Analysis Results:\n');

F Analysis Results:

fprintf('Location (km) | Dip  (degrees) | D_faults (km)\n');

Location (km) | Dip  (degrees) | D_faults (km)

for i = 1:length(f_locations)
    fprintf('%10.2f | %17.2f | %10.2f\n', f_locations(i), f_dip_angles(i), D_faults(i));
end

      7.00 |             58.47 |       2.87
     14.00 |             90.00 |       4.48
     23.00 |            -69.79 |      -4.34

% Ensure both variables are column vectors of the same size
f_locations  = f_locations(:);   % Convert to column vector
f_dip_angles = f_dip_angles(:); % Convert to column vector
% Concatenate and write to Excel
xlswrite('f_analysis_results.xlsx', [f_locations, f_dip_angles]);

Warning: Unable to write to Excel format, attempting to write file to CSV format.

% Save results to Excel (only if fs are detected)
if ~isempty(f_locations)
    xlswrite('f_analysis_results.xlsx', [f_locations, f_dip_angles]);
else
    warning('No significant fs detected.');
end

Warning: Unable to write to Excel format, attempting to write file to CSV format.

[Mathieu NOE]

hello

your code is asking for some numeric inputs - what should we use to get the correct graph ?

% Input number of layers and densities
num_blocks = input('Enter the number of blocks: ');
block_densities = zeros(1, num_blocks);
for i = 1:num_blocks
    block_densities(i) = input(['Enter the density of block ', num2str(i), ' (kg/m^3): ']);
end

[Moustafa]

Put the block_densities =5

block_densities(i) = [2.40, 2.45, 2.50, 2.55, 2.60];

[Mathieu NOE]

hello again

so ok we have now a working code but we lack some explanations about what is the logic to change your lines positions in your plot ? is the code wrong somewhere / what is the problem ? can you explain a bit more ?

[Moustafa]

The script treats a group of rock layers that were subjected to geological faults (breaking or rift) through plans that are fully represented in Figure (3) from the left to the right as a normal plan, a vertical plan, and a reverse plan. These plans were named, depending on the dip angle of the plans and upthrown and downthrown on both sides of the rift ... And when drawing the plans , they shifted from the right places with the exception of the vertical level in its right place. This is my problem

[Mathieu NOE]

hello again

it's a bit annoying that you created so many variables that basically are related to the same data. It does not help understand how to dela with your problem

there is one first point which is not clear to me is wheter this values need to be upadated (and how..) , based on the other data (A) :

% Define f locations and corresponding dip

f_locations = [7.00, 14.00, 23.00];

d_angles = [58.47, 90.00, -69.79];

because it seems the wide black lines are not properly positionned

here a much simplified a shortened code to work on :

clc; close all; clear;
% Load the Excel data
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename, 'Sheet1');
% Extract relevant columns
x = data(:, 1); % Distance (x in km)
Field = data(:, 2); % Earth filed (Field in unit)
%==============================================================%
% % Input number of layers and densities
% num_blocks = input('Enter the number of blocks: ');
% block_densities = zeros(1, num_blocks);
% for i = 1:num_blocks
%     block_densities(i) = input(['Enter the density of block ', num2str(i), ' (kg/m^3): ']);
% end
% EDIT to run code here
num_blocks = 5;
block_densities = [2.40, 2.45, 2.50, 2.55, 2.60];
%==============================================================%
% Constants
G = 0.00676; 
Lower_density = 2.67; % in kg/m^3
%==============================================================%
% Calculate inverted depth profile for each layer
z_inv = zeros(length(x), num_blocks);
for i = 1:num_blocks
    density_contrast = block_densities(i) - Lower_density;
    if density_contrast ~= 0
        z_inv(:, i) = Field ./ (2 * pi * G * density_contrast);
    else
        z_inv(:, i) = NaN; % Avoid division by zero
    end
end
%==============================================================%
% Compute vertical gradient (VG) of inverted depth (clean)
VG = diff(z_inv(:, 1)) ./ diff(x);
%==============================================================%
% Set fault threshold and find f indices based on d changes
f_threshold = 0.5; % Threshold for identifying significant d changes
f_indices = find(abs(diff(z_inv(:, 1))) > f_threshold);
%==============================================================%
% Initialize f locations and dip arrays
%==============================================================%
f_locations = x(f_indices);  % Automatically determined f locations
f_dip_angles = nan(size(f_indices)); % Placeholder for calculated dip 
% Calculate dip  for each identified f
for i = 1:length(f_indices)
    idx = f_indices(i);
    if idx < length(x)
        f_dip_angles(i) = atand(abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / (x(idx + 1) - x(idx)));
    else
        f_dip_angles(i) = atand(abs(z_inv(idx, 1) - z_inv(idx - 1, 1)) / (x(idx) - x(idx - 1)));
    end
end
%==============================================================%
%                    Displacement of faults 
%==============================================================%
D_faults = zeros(size(f_dip_angles));
for i = 1:length(f_indices)
    idx = f_indices(i);
    dip_angle_rad = deg2rad(f_dip_angles(i)); % Convert dip to radians
    D_faults(i) = abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / sin(dip_angle_rad);
end
% Assign displacement values
D1 = D_faults(1); % NF displacemen
D2 = D_faults(2); % VF displacement
D3 = D_faults(3); % RF displacement
%==============================================================%
%                 Processing Data for Interpretation
%==============================================================%
A = [x Field z_inv];  % New Data Obtained
col_names = {'x', 'Field'};
for i = 1:num_blocks
    col_names{end+1} = ['z', num2str(i)];
end
dataM = array2table(A, 'VariableNames', col_names);
t1 = dataM;
[nr, nc] = size(t1);
t1_bottoms = t1;
for jj = 3:nc
    for ii = 1:nr-1
        if t1_bottoms{ii, jj} ~= t1_bottoms{ii+1, jj}
            t1_bottoms{ii, jj} = NaN;
        end
    end
end
%==============================================================%
%                     Identifying NaN rows
%==============================================================%
nans = isnan(t1_bottoms{:, 3:end});
nan_rows = find(any(nans, 2));
xc = t1_bottoms{nan_rows, 1}; % Corrected x-coordinates
yc = zeros(numel(nan_rows), 1); % y-coordinates for NaN rows
for ii = 1:numel(nan_rows)
    idx = find(~nans(nan_rows(ii), :), 1, 'last');
    if isempty(idx)
        yc(ii) = 0;
    else
        yc(ii) = t1_bottoms{nan_rows(ii), idx+2};
    end
end

%==============================================================%
% Plot the interpreted d profiles
figure(3)
hold on
yl = [0 ceil(max(z_inv(:)))];  % Get Y-axis limits

% Define f locations and corresponding dip 
% ARE THESE VALUES CORRECT ??
f_locations = [7.00, 14.00, 23.00];  
d_angles = [58.47, 90.00, -69.79];    

for ii = 1:numel(f_locations)
    % Find the nearest x index for each fault location
    [~, idx] = min(abs(x - f_locations(ii)));
    
    % Get the starting x and y coordinates (fault starts at the surface)
    x_f = t1_bottoms{idx, 1};  
    y_f = yl(1); % Start at the surface (0 km depth)
    
    % Check if the dip angle is 90° (vf)
    if d_angles(ii) == 90
        x_r = [x_f x_f]; % Vertical line
        y_r = [yl(1), yl(2)];    % From surface to de limit
    else
        % Convert dip to slope (m = tan(angle))
        m = tand(d_angles(ii));             
        
        % Define the x range for fault line
        x_r = x;
        y_r = y_f - m * (x_r - x_f);  % Line equation  (+ m)
        
        % Clip y_range within the plot limits
        y_r(y_r > yl(2)) = yl(2);
        y_r(y_r < yl(1)) = yl(1);
        
    end
    
    % Plot the fault lines in black (matching the image)
    plot(x_r, y_r, 'k', 'LineWidth', 3)
    
end

% find indexes of start / end groups of data
zz = sum(z_inv,2);

plot(x,z_inv, 'LineWidth', 1)

set(gca, 'YDir', 'reverse')  % Reverse Y-axis for depth representation
box on
grid on
xlabel('Distance (km)'); 
ylabel('D (km)'); 
title('New Interpretation of Profile Data')

according to the posted image , if we compare the black and red lines => there is obviously something to do here

so it's still unclear to me if we have one or multiple corrections to do :

• correct thick angled black lines to match red lines position (how ? )
• then make sure the horizontal lines are touching the red lines

[Mathieu NOE]

just for my own fun I changed a bit the code and also the values of

f_locations = [8.50, 14.00, 21.00];

d_angles = [85, 90.00, -78];

is it to be this way ?

clc; close all; clear;
% Load the Excel data
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename, 'Sheet1');
% Extract relevant columns
x = data(:, 1); % Distance (x in km)
Field = data(:, 2); % Earth filed (Field in unit)
%==============================================================%
% % Input number of layers and densities
% num_blocks = input('Enter the number of blocks: ');
% block_densities = zeros(1, num_blocks);
% for i = 1:num_blocks
%     block_densities(i) = input(['Enter the density of block ', num2str(i), ' (kg/m^3): ']);
% end
% EDIT to run code here
num_blocks = 5;
block_densities = [2.40, 2.45, 2.50, 2.55, 2.60];
%==============================================================%
% Constants
G = 0.00676; 
Lower_density = 2.67; % in kg/m^3
%==============================================================%
% Calculate inverted depth profile for each layer
z_inv = zeros(length(x), num_blocks);
for i = 1:num_blocks
    density_contrast = block_densities(i) - Lower_density;
    if density_contrast ~= 0
        z_inv(:, i) = Field ./ (2 * pi * G * density_contrast);
    else
        z_inv(:, i) = NaN; % Avoid division by zero
    end
end
%==============================================================%
% Compute vertical gradient (VG) of inverted depth (clean)
VG = diff(z_inv(:, 1)) ./ diff(x);
%==============================================================%
% Set fault threshold and find f indices based on d changes
f_threshold = 0.5; % Threshold for identifying significant d changes
f_indices = find(abs(diff(z_inv(:, 1))) > f_threshold);
%==============================================================%
% Initialize f locations and dip arrays
%==============================================================%
f_locations = x(f_indices);  % Automatically determined f locations
f_dip_angles = nan(size(f_indices)); % Placeholder for calculated dip 
% Calculate dip  for each identified f
for i = 1:length(f_indices)
    idx = f_indices(i);
    if idx < length(x)
        f_dip_angles(i) = atand(abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / (x(idx + 1) - x(idx)));
    else
        f_dip_angles(i) = atand(abs(z_inv(idx, 1) - z_inv(idx - 1, 1)) / (x(idx) - x(idx - 1)));
    end
end
%==============================================================%
%                    Displacement of faults 
%==============================================================%
D_faults = zeros(size(f_dip_angles));
for i = 1:length(f_indices)
    idx = f_indices(i);
    dip_angle_rad = deg2rad(f_dip_angles(i)); % Convert dip to radians
    D_faults(i) = abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / sin(dip_angle_rad);
end
% Assign displacement values
D1 = D_faults(1); % NF displacemen
D2 = D_faults(2); % VF displacement
D3 = D_faults(3); % RF displacement
%==============================================================%
%                 Processing Data for Interpretation
%==============================================================%
A = [x Field z_inv];  % New Data Obtained
col_names = {'x', 'Field'};
for i = 1:num_blocks
    col_names{end+1} = ['z', num2str(i)];
end
dataM = array2table(A, 'VariableNames', col_names);
t1 = dataM;
[nr, nc] = size(t1);
t1_bottoms = t1;
for jj = 3:nc
    for ii = 1:nr-1
        if t1_bottoms{ii, jj} ~= t1_bottoms{ii+1, jj}
            t1_bottoms{ii, jj} = NaN;
        end
    end
end
%==============================================================%
%                     Identifying NaN rows
%==============================================================%
nans = isnan(t1_bottoms{:, 3:end});
nan_rows = find(any(nans, 2));
xc = t1_bottoms{nan_rows, 1}; % Corrected x-coordinates
yc = zeros(numel(nan_rows), 1); % y-coordinates for NaN rows
for ii = 1:numel(nan_rows)
    idx = find(~nans(nan_rows(ii), :), 1, 'last');
    if isempty(idx)
        yc(ii) = 0;
    else
        yc(ii) = t1_bottoms{nan_rows(ii), idx+2};
    end
end

%==============================================================%
% Plot the interpreted d profiles
figure(3)

% Plot the fault lines in black (matching the image)
plot(x, zeros(size(x)), 'k', 'LineWidth', 3)
    
hold on
yl = [0 ceil(max(z_inv(:)))];  % Get Y-axis limits

% Define f locations and corresponding dip 
% ARE THESE VALUES CORRECT ??
% f_locations = [7.00, 14.00, 23.00];  
% d_angles = [58.47, 90.00, -69.79];    

% suggestion ? 
f_locations = [8.50, 14.00, 21.00];  
d_angles = [85, 90.00, -78];    

for ii = 1:numel(f_locations)
    % Find the nearest x index for each fault location
    [~, idx] = min(abs(x - f_locations(ii)));
    
    % Get the starting x and y coordinates (fault starts at the surface)
    x_f = t1_bottoms{idx, 1};  
    y_f = yl(1); % Start at the surface (0 km depth)
    
    % Check if the dip angle is 90° (vf)
    if d_angles(ii) == 90
        x_r = [x_f x_f]; % Vertical line
        y_r = [yl(1), yl(2)];    % From surface to de limit
    else
        % Convert dip to slope (m = tan(angle))
        m = tand(d_angles(ii));             
       
        % Define the x range for fault line
        y_r = linspace(yl(1),yl(2),25);
        x_r = x_f - y_r/m;
    end
    
    % Plot the fault lines in black (matching the image)
    plot(x_r, y_r, 'r', 'LineWidth', 3)
    
end

plot(x,z_inv, 'LineWidth', 1)

set(gca, 'YDir', 'reverse')  % Reverse Y-axis for depth representation
box on
grid on
xlabel('Distance (km)'); 
ylabel('D (km)'); 
title('New Interpretation of Profile Data')

[Moustafa]

Your answer Sir will be accepted if you find automatical way to locate the fault plane as you do and correct the extension of brocken blocks with respect to the fault blane as shown in attached file

[Mathieu NOE]

hello again

yeap - the only "little" issue is to undertsand how to place the red lines as there is no clear information yet how we are suppose to place them (and then extend the horizontal lines) . I could do it with arbitrary choices but I doubt this is the right approach - especially if you wnt to use the code in "automatic" mode... so is there any publication or article that could help figure out how to solve this problem ?

Answer 1

[Mathieu NOE]

so finally and hopefully I found somehow a path to the solution

first problem was to find the best position and slopes for the thick black lines

so instead of using hard coded values :

% Define f locations and corresponding dip

f_locations = [7.00, 14.00, 23.00];

d_angles = [58.47, 90.00, -69.79];

I opted for a kind of linear fit (polyfit) to get these infos , based on data samples that are supposed to represent the start / end portion of each block

this is the result so far (plotted in red thick lines)

clc; close all; clear;
% Load the Excel data
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename, 'Sheet1');
% Extract relevant columns
x = data(:, 1); % Distance (x in km)
Field = data(:, 2); % Earth filed (Field in unit)
%==============================================================%
% % Input number of layers and densities
% num_blocks = input('Enter the number of blocks: ');
% block_densities = zeros(1, num_blocks);
% for i = 1:num_blocks
%     block_densities(i) = input(['Enter the density of block ', num2str(i), ' (kg/m^3): ']);
% end
% EDIT to run code here
num_blocks = 5;
block_densities = [2.40, 2.45, 2.50, 2.55, 2.60];
%==============================================================%
% Constants
G = 0.00676; 
Lower_density = 2.67; % in kg/m^3
%==============================================================%
% Calculate inverted depth profile for each layer
z_inv = zeros(length(x), num_blocks);
for i = 1:num_blocks
    density_contrast = block_densities(i) - Lower_density;
    if density_contrast ~= 0
        z_inv(:, i) = Field ./ (2 * pi * G * density_contrast);
    else
        z_inv(:, i) = NaN; % Avoid division by zero
    end
end
%==============================================================%
% Compute vertical gradient (VG) of inverted depth (clean)
VG = diff(z_inv(:, 1)) ./ diff(x);
%==============================================================%
% Set fault threshold and find f indices based on d changes
f_threshold = 0.5; % Threshold for identifying significant d changes
f_indices = find(abs(diff(z_inv(:, 1))) > f_threshold);
%==============================================================%
% Initialize f locations and dip arrays
%==============================================================%
f_locations = x(f_indices);  % Automatically determined f locations
f_dip_angles = nan(size(f_indices)); % Placeholder for calculated dip 
% Calculate dip  for each identified f
for i = 1:length(f_indices)
    idx = f_indices(i);
    if idx < length(x)
        f_dip_angles(i) = atand(abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / (x(idx + 1) - x(idx)));
    else
        f_dip_angles(i) = atand(abs(z_inv(idx, 1) - z_inv(idx - 1, 1)) / (x(idx) - x(idx - 1)));
    end
end
%==============================================================%
%                    Displacement of faults 
%==============================================================%
D_faults = zeros(size(f_dip_angles));
for i = 1:length(f_indices)
    idx = f_indices(i);
    dip_angle_rad = deg2rad(f_dip_angles(i)); % Convert dip to radians
    D_faults(i) = abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / sin(dip_angle_rad);
end
% Assign displacement values
D1 = D_faults(1); % NF displacemen
D2 = D_faults(2); % VF displacement
D3 = D_faults(3); % RF displacement
%==============================================================%
%                 Processing Data for Interpretation
%==============================================================%
A = [x Field z_inv];  % New Data Obtained
col_names = {'x', 'Field'};
for i = 1:num_blocks
    col_names{end+1} = ['z', num2str(i)];
end
dataM = array2table(A, 'VariableNames', col_names);
t1 = dataM;
[nr, nc] = size(t1);
t1_bottoms = t1;
for jj = 3:nc
    for ii = 1:nr-1
        if t1_bottoms{ii, jj} ~= t1_bottoms{ii+1, jj}
            t1_bottoms{ii, jj} = NaN;
        end
    end
end
%==============================================================%
%                     Identifying NaN rows
%==============================================================%
nans = isnan(t1_bottoms{:, 3:end});
nan_rows = find(any(nans, 2));
xc = t1_bottoms{nan_rows, 1}; % Corrected x-coordinates
yc = zeros(numel(nan_rows), 1); % y-coordinates for NaN rows
for ii = 1:numel(nan_rows)
    idx = find(~nans(nan_rows(ii), :), 1, 'last');
    if isempty(idx)
        yc(ii) = 0;
    else
        yc(ii) = t1_bottoms{nan_rows(ii), idx+2};
    end
end

%==============================================================%
% Plot the interpreted d profiles
figure(3)

% Plot the fault lines in black (matching the image)
plot(x, zeros(size(x)), 'k', 'LineWidth', 3) % top horizontal lines
hold on
yl = [0 ceil(max(z_inv(:)))];  % Get Y-axis limits

% fault lines are estimated with polyfit
for kk = 1:numel(f_indices)
    xa = x(f_indices(kk));
    xb = x(f_indices(kk)+1);
    ya = z_inv(f_indices(kk),:);
    yb = z_inv(f_indices(kk)+1,:);
    p = polyfit([ya yb],[xa*ones(1,5) xb*ones(1,5)],1);
    y_r = linspace(yl(1),yl(2),25);
    xf = polyval(p,y_r);
    % Plot the fault lines in red (matching the image)
    plot(xf, y_r, 'r', 'LineWidth', 3)
end


plot(x,z_inv,'-*', 'LineWidth', 1)

set(gca, 'YDir', 'reverse')  % Reverse Y-axis for depth representation
box on
grid on
xlabel('Distance (km)'); 
ylabel('D (km)'); 
title('New Interpretation of Profile Data')

●

[Mathieu NOE]

then the last challenge was to remove the points that are on the wrong side of the red lines and do the nice horizontal plot , connecting to the red lines

had to scratch my head to find a way to get this plot as I wanted , but glad to post it here :

clc; close all; clear;
% Load the Excel data
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename, 'Sheet1');
% Extract relevant columns
x = data(:, 1); % Distance (x in km)
Field = data(:, 2); % Earth filed (Field in unit)
%==============================================================%
% % Input number of layers and densities
% num_blocks = input('Enter the number of blocks: ');
% block_densities = zeros(1, num_blocks);
% for i = 1:num_blocks
%     block_densities(i) = input(['Enter the density of block ', num2str(i), ' (kg/m^3): ']);
% end
% EDIT to run code here
num_blocks = 5;
block_densities = [2.40, 2.45, 2.50, 2.55, 2.60];
%==============================================================%
% Constants
G = 0.00676; 
Lower_density = 2.67; % in kg/m^3
%==============================================================%
% Calculate inverted depth profile for each layer
z_inv = zeros(length(x), num_blocks);
for i = 1:num_blocks
    density_contrast = block_densities(i) - Lower_density;
    if density_contrast ~= 0
        z_inv(:, i) = Field ./ (2 * pi * G * density_contrast);
    else
        z_inv(:, i) = NaN; % Avoid division by zero
    end
end
%==============================================================%
% Compute vertical gradient (VG) of inverted depth (clean)
VG = diff(z_inv(:, 1)) ./ diff(x);
%==============================================================%
% Set fault threshold and find f indices based on d changes
f_threshold = 0.5; % Threshold for identifying significant d changes
f_indices = find(abs(diff(z_inv(:, 1))) > f_threshold);
%==============================================================%
% Initialize f locations and dip arrays
%==============================================================%
f_locations = x(f_indices);  % Automatically determined f locations
f_dip_angles = nan(size(f_indices)); % Placeholder for calculated dip 
% Calculate dip  for each identified f
for i = 1:length(f_indices)
    idx = f_indices(i);
    if idx < length(x)
        f_dip_angles(i) = atand(abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / (x(idx + 1) - x(idx)));
    else
        f_dip_angles(i) = atand(abs(z_inv(idx, 1) - z_inv(idx - 1, 1)) / (x(idx) - x(idx - 1)));
    end
end
%==============================================================%
%                    Displacement of faults 
%==============================================================%
D_faults = zeros(size(f_dip_angles));
for i = 1:length(f_indices)
    idx = f_indices(i);
    dip_angle_rad = deg2rad(f_dip_angles(i)); % Convert dip to radians
    D_faults(i) = abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / sin(dip_angle_rad);
end
% Assign displacement values
D1 = D_faults(1); % NF displacemen
D2 = D_faults(2); % VF displacement
D3 = D_faults(3); % RF displacement
%==============================================================%
%                 Processing Data for Interpretation
%==============================================================%
A = [x Field z_inv];  % New Data Obtained
col_names = {'x', 'Field'};
for i = 1:num_blocks
    col_names{end+1} = ['z', num2str(i)];
end
dataM = array2table(A, 'VariableNames', col_names);
t1 = dataM;
[nr, nc] = size(t1);
t1_bottoms = t1;
for jj = 3:nc
    for ii = 1:nr-1
        if t1_bottoms{ii, jj} ~= t1_bottoms{ii+1, jj}
            t1_bottoms{ii, jj} = NaN;
        end
    end
end
%==============================================================%
%                     Identifying NaN rows
%==============================================================%
nans = isnan(t1_bottoms{:, 3:end});
nan_rows = find(any(nans, 2));
xc = t1_bottoms{nan_rows, 1}; % Corrected x-coordinates
yc = zeros(numel(nan_rows), 1); % y-coordinates for NaN rows
for ii = 1:numel(nan_rows)
    idx = find(~nans(nan_rows(ii), :), 1, 'last');
    if isempty(idx)
        yc(ii) = 0;
    else
        yc(ii) = t1_bottoms{nan_rows(ii), idx+2};
    end
end

%==============================================================%
% Plot the interpreted d profiles
figure(3)

% Plot the fault lines in black (matching the image)
plot(x, zeros(size(x)), 'k', 'LineWidth', 3) % top horizontal lines
set(gca, 'YDir', 'reverse')  % Reverse Y-axis for depth representation
box on
grid on
xlabel('Distance (km)'); 
ylabel('D (km)'); 
title('New Interpretation of Profile Data')
hold on
yl = [0 ceil(max(z_inv(:)))];  % Get Y-axis limits

% fault lines are estimated with polyfit
ind_nan = [];
for kk = 1:numel(f_indices)
    xtmp = x;
    xa = x(f_indices(kk));
    xb = x(f_indices(kk)+1);
    ya = z_inv(f_indices(kk),:);
    yb = z_inv(f_indices(kk)+1,:);
    p{kk} = polyfit([ya yb],[xa*ones(1,5) xb*ones(1,5)],1);
    xf = polyval(p{kk},yl);
    % Plot the fault lines in red (matching the image)
    plot(xf, yl, 'r', 'LineWidth', 3)
    ind1 = find(xtmp>min(xf)-1, 1 );
    ind2 = find(xtmp<max(xf)+1, 1, 'last' );
    ind_nan = [ind_nan [ind1: ind2]];
end

x(ind_nan) = [];
z_inv(ind_nan,:) = [];

% complete horizontal lines datas
xnew = x;
z_inv_new = z_inv;

for kk = 1:numel(f_indices)
    xf = polyval(p{kk},yl);
    ind1(kk) = find(x>min(xf)-1, 1 );
    ind2(kk) = find(x<max(xf)+1, 1, 'last' );
    ya = z_inv(ind1(kk),:);
    yb = z_inv(ind2(kk),:);
   
    % replace ends of z_inv data to match red lines 
    % define new end values 
    xa = polyval(p{kk},ya);
    ytmp = diag(ya);
    ytmp(abs(ytmp)<eps) = NaN;
    % add new end values
    xnew = [xnew; xa(:)];
    z_inv_new = [z_inv_new; ytmp];
    % define new start values 
    xb = polyval(p{kk},yb);
    ytmp = diag(yb);
    ytmp(abs(ytmp)<eps) = NaN;
    % add new start values
    xnew = [xnew; xb(:)];
    z_inv_new = [z_inv_new; ytmp]; 
   
end

%% finally , the plot !!
% have to plot each individual segment otherwise the plot is messy
[z_inv_unic,ia,ic] = uniquetol(z_inv_new(:),1e-6);
ii = isnan(z_inv_unic);
z_inv_unic(ii) = [];
ia(ii) = [];

for kk = 1:numel(z_inv_unic)
    [r,c] = find(z_inv_new == z_inv_unic(kk));
    for m = 1:numel(c)
        plot(xnew(r),z_inv_new(r,c),'-', 'LineWidth', 1)
    end
end

●