Shifting the parallel lines

I need to shift the lines (as in annexed figure based on equl dip angle of lines of the following script :

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)
y = data(:, 2);    % Bouguer gravity anomaly (g in mGal)

% Constants
k = 0.00676;               
rho = -0.220;        % Density contrast (delta_rho in kg/m^3)

% Calculate inverted depth profile 
z_inv_calculated = y./ (2 * pi * k * rho);

% Set fault threshold and find fault indices based on depth changes
fault_threshold = 0.5;  % Threshold for identifying significant depth changes
fault_indices = find(abs(diff(z_inv_calculated)) > fault_threshold);

% Initialize fault locations and dip angles arrays
fault_locations = x(fault_indices);          % Automatically determined fault locations (x-coordinates)
fault_dip_angles = nan(size(fault_indices)); % Placeholder for calculated dip angles

% Calculate dip angles for each identified fault
for i = 1:length(fault_indices)
    idx = fault_indices(i);
    
    % Calculate dip angle using adjacent points
    if idx < length(x)
        fault_dip_angles(i) = atand(abs(z_inv_calculated(idx + 1) - z_inv_calculated(idx)) / (x(idx + 1) - x(idx)));
    else
        fault_dip_angles(i) = atand(abs(z_inv_calculated(idx) - z_inv_calculated(idx - 1)) / (x(idx) - x(idx - 1)));
    end
end

% Determine max depth for fault extension lines based on data range
z_max_depth = max(abs(z_inv_calculated)) * 1.2;  % Extend slightly beyond deepest depth

% Plotting the results
figure;
ax1 = subplot(2,1,1);
plot(x, y, 'r', 'LineWidth', 2);
title('Bouguer Gravity Anomaly Profile');
xlabel('Distance ');
ylabel('y_{y} ');
grid on;

ax2 = subplot(2,1,2);
plot(x, z_inv_calculated, 'b', 'LineWidth', 2);
title('Inverted Depth Profile with Fault Extensions to Surface');
xlabel('Distance ');
ylabel('z_{inv} ');
set(gca, 'YDir', 'reverse'); % Reverse y-axis for depth
grid on;
set(gca, 'YDir','reverse')
hold on;

% Extend fault planes to surface and max depth using calculated dip angles
for i = 1:length(fault_locations)
    x_fault = fault_locations(i);
    theta = fault_dip_angles(i);
    
    % Calculate extension points based on dip angle
    if theta ~= 90
        % Calculate the end points of the fault line
        x_end = x_fault - (z_max_depth / tand(theta));
        plot([x_fault, x_end], [0, z_max_depth], 'r--', 'LineWidth', 1.5); % Red dashed line
    else
        % Vertical fault (90 degrees)
        plot([x_fault, x_fault], [0, z_max_depth], 'r--', 'LineWidth', 1.5); % Red dashed line
    end
end

hold off;

% Display results
disp('Fault Analysis Results:');

Fault Analysis Results:

disp('Location (km) | Dip Angle (degrees)');

Location (km) | Dip Angle (degrees)

for i = 1:length(fault_locations)
    fprintf('%10.2f        | %5.2f\n', fault_locations(i), fault_dip_angles(i));
end

      7.00        | 63.43
     14.00        | 90.00
     23.00        | -73.30

% Set same x-axis scale
xlim_range = [min(x) max(x)];  
set(ax1, 'XLim', xlim_range);
set(ax2, 'XLim', xlim_range);
linkaxes([ax1, ax2], 'x'); % Link all subplots' x-axes

% Tabulated results (including NaN for non-fault points)
table_results = table(y, z_inv_calculated, 'VariableNames', {'y', 'z_inv_km'});

[Matt J]

And? What's wrong with the code you've posted?

[dpb]

@Moustafa, you forgot to attach the data file as well...

[Moustafa]

Dear Sir,

Sorry for forgetten attached data

please see the attached file

Answer 1

[Voss]

See below.

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)
y = data(:, 2);    % Bouguer gravity anomaly (g in mGal)

% Constants
k = 0.00676;               
rho = -0.220;        % Density contrast (delta_rho in kg/m^3)

% Calculate inverted depth profile 
z_inv_calculated = y./ (2 * pi * k * rho);

% Set fault threshold and find fault indices based on depth changes
fault_threshold = 0.5;  % Threshold for identifying significant depth changes
fault_indices = find(abs(diff(z_inv_calculated)) > fault_threshold);

% Initialize fault locations and dip angles arrays
fault_locations = x(fault_indices);          % Automatically determined fault locations (x-coordinates)
fault_dip_angles = nan(size(fault_indices)); % Placeholder for calculated dip angles

% Calculate dip angles for each identified fault
for i = 1:length(fault_indices)
    idx = fault_indices(i);
    
    % Calculate dip angle using adjacent points
    if idx < length(x)
        fault_dip_angles(i) = atand(abs(z_inv_calculated(idx + 1) - z_inv_calculated(idx)) / (x(idx + 1) - x(idx)));
    else
        fault_dip_angles(i) = atand(abs(z_inv_calculated(idx) - z_inv_calculated(idx - 1)) / (x(idx) - x(idx - 1)));
    end
end

% Determine max depth for fault extension lines based on data range
z_max_depth = max(abs(z_inv_calculated)) * 1.2;  % Extend slightly beyond deepest depth
z_min_depth = 0;

% Plotting the results
figure;
ax1 = subplot(2,1,1);
plot(x, y, 'r', 'LineWidth', 2);
title('Bouguer Gravity Anomaly Profile');
xlabel('Distance ');
ylabel('y_{y} ');
grid on;

ax2 = subplot(2,1,2);
plot(x, z_inv_calculated, 'b', 'LineWidth', 2);
title('Inverted Depth Profile with Fault Extensions to Surface');
xlabel('Distance ');
ylabel('z_{inv} ');
set(gca, 'YDir', 'reverse'); % Reverse y-axis for depth
grid on;
set(gca, 'YDir','reverse')
hold on;

tan_theta = tand(fault_dip_angles);

% Extend fault planes to surface and max depth using calculated dip angles
for i = 1:length(fault_locations)
    
    x_fault = fault_locations(i);
    theta = fault_dip_angles(i);
    
    % Calculate extension points based on dip angle
    if theta == 90
        % Vertical fault (90 degrees)
        x_start = x_fault;
        x_end = x_fault;
    else
        % Calculate the end points of the fault line
        idx = fault_indices(i);
        x_start = x_fault + (z_inv_calculated(idx)-z_min_depth)/tan_theta(i);
        x_end = x_fault + (z_inv_calculated(idx)-z_max_depth)/tan_theta(i);
    end
    
    plot([x_start, x_end], [0, z_max_depth], 'r--', 'LineWidth', 1.5); % Red dashed line
    
end

hold off;

% Display results
disp('Fault Analysis Results:');

Fault Analysis Results:

disp('Location (km) | Dip Angle (degrees)');

Location (km) | Dip Angle (degrees)

for i = 1:length(fault_locations)
    fprintf('%10.2f        | %5.2f\n', fault_locations(i), fault_dip_angles(i));
end

      7.00        | 63.43
     14.00        | 90.00
     23.00        | -73.30

% Set same x-axis scale
xlim_range = [min(x) max(x)];  
set(ax1, 'XLim', xlim_range);
set(ax2, 'XLim', xlim_range);
linkaxes([ax1, ax2], 'x'); % Link all subplots' x-axes

●

% Tabulated results (including NaN for non-fault points)
table_results = table(y, z_inv_calculated, 'VariableNames', {'y', 'z_inv_km'});

[Moustafa]

Dear Sir;

Thanks for quick response with accepted solution

[Voss]

You're welcome!