Need help to correct and modify MATLAB code

I have the following code :

%===============================%
clc
close all;
clear;

% Define the given arrays
x = [0.02 0.55 1.00 1.4 1.4 2.32 2.62 ...
    3.77 4.81 5.61 6.29 6.29 6.00 ...
    6.00 7.00 7.00 7.00];
z = [5 5 5 5 5 4 4 4 4 4 4 4 2 2 2 2 2];

% Find sudden changes in depth
sudden_changes = diff(z);

% Find indices where sudden changes occur
fault_indices = find(sudden_changes ~= 0) + 1;

% Plot the top of the rock formation
plot(x, z, '-o');
hold on;

% Plot fault lines
for i = 1:length(fault_indices)
    idx = fault_indices(i);
    plot([x(idx-1), x(idx)], [z(idx-1), z(idx)], 'r--');
end

% Calculate dip angle using the tan-rule
dip_angles = atand(diff(z) ./ diff(x));

% Identify normal and reverse faults
normal_faults = dip_angles <= 90;
reverse_faults = dip_angles > 90;

% Plot normal faults
scatter(x(fault_indices(normal_faults)), z(fault_indices(normal_faults)), 'g', 'filled');

The logical indices contain a true value outside of the array bounds.

% Plot reverse faults
scatter(x(fault_indices(reverse_faults)), z(fault_indices(reverse_faults)), 'm', 'filled');

% Add labels and legend
xlabel('X-axis');
ylabel('Depth');
title('Top of the rock formation with fault lines');
legend('Top of the rock formation', 'Fault lines', 'Normal faults', 'Reverse faults');

set(gca, 'YDir', 'reverse');
grid on;
%===============================%

I need to draw the Earth's surface where the depth is constan surface and equal to zero

and need th extend the the fault lines to the Earths surface as shown in attached figure

and need to detect the type of fault from intersected fault line with Earth's surface

[Walter Roberson]

fault_indices = find(sudden_changes ~= 0) + 1;

That is variable length, only the locations where sudden_changes is non-zero

dip_angles = atand(diff(z) ./ diff(x));

That is one element shorter than z (or x)

normal_faults = dip_angles <= 90;

Same length as dip_angles, so one element shorter than z (or x)

scatter(x(fault_indices(normal_faults)), z(fault_indices(normal_faults)), 'g', 'filled');

fault_indices is short, normal_faults is full length, fault_indices indexed at normal_faults is a problem.

Answer 1

[Sandeep Mishra]

Hi Moustafa,

From the provided code snippet, it is evident that you are trying to plot the data and based on the varying values of the data, you want to determine the type of fault.

You can implement the functionalities as outlined below:

1. You can use the MATLAB 'plot' function to represent the Earth's surface, adjusting the 'LineWidth' parameter to define the line thickness as follows:

x_range = 0:max(x); 
plot(x_range, zeros(size(x_range, 2)), 'Color', 'k', 'LineWidth',5); 
ylim([0, max(z)]);

2. To extend the faults to Earth’s surface, you can apply the ‘polyfit’ function along with the ‘roots’ function inside the for loop to compute the x-intercept of the line, then plot a line up to the x-intercept value, extending the fault line to the Earth's surface:

poly_res = polyfit(x(idx-1:idx), z(idx-1:idx), 1);   
x_intercept = roots(poly_res);
plot([x(idx-1), x_intercept], [z(idx-1), 0], 'r-');

3. To detect the fault type, you can use the 'atan2d' function to avoid NaN values that occur with 'atand' function. Based on the absolute values, you can classify the fault type as follows:

dip_angles = abs(atan2d(diff(z), diff(x)));
normal_faults = dip_angles~=0 & dip_angles <= 90; 
reverse_faults = dip_angles~=0 &  dip_angles > 90; 
scatter(x(normal_faults), z(normal_faults), 'g', 'filled'); 
scatter(x(reverse_faults), z(reverse_faults), 'm', 'filled'); 

For more information, refer to the following MathWorks Documentation:

1. ‘polyfit’ function: https://www.mathworks.com/help/matlab/ref/polyfit.html
2. ‘roots’ function: https://www.mathworks.com/help/matlab/ref/roots.html
3. ‘atan2d’ function: https://www.mathworks.com/help/matlab/ref/double.atan2d.html

I hope this helps you!