Draw lines (representing planes) based on dipping anle of planes

In the following script I need to plot figure 3 to be corrected as in figure of the attached file depending on the dipping planes:

the number of blocks = 5

the input densiities 2.40, 2.45, 2.50, 2.55, 2.60

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

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: ');

Error using input
Support for user input is required, which is not available on this platform.

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

%==============================================================%
% 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');
fprintf('Location (km) | Dip  (degrees) | D_faults (km)\n');

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

% 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]);
% 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

Answer 1

[Moustafa]

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

%==============================================================%
% 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');
fprintf('Location (km) | Dip  (degrees) | D_faults (km)\n');

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

% 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]);



% 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

Answer 2

[Moustafa]

The attached two files again