Cross-section (Need help to correct the code or modify it )

Dear all

I need to correct or modify the following code to give the attached figure: clc

close all;
clear;

% Define the given arrays

x1 =[0, 1, 2, 3, 3.7];
z1 =[6, 6, 6, 6, 6];
x2 =[4.1, 5, 6, 7, 8, 9];
z2 =[4, 4, 4, 4, 4, 4];
x3 =[9, 10, 10, 11,12, 13, 14, 15, 16, 16.4];
z3 =[6, 6, 6, 6, 6, 6, 6, 6, 6, 6];
x4 =[15.3, 16, 17, 18, 19];
z4 =[3, 3, 3, 3, 3];

x =[x1,x2,x3,x4];
z =[z1,z2,z3,z4];


%x = [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21];
%z = [2 2 2 3 3 3 3 3 3 3 4 4 4 4 4 4 1 1 1 1 1 1];

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

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

% Initialize variables to store fault information
fault_locations = [];
fault_angles = [];
fault_types = {};

% Plot the top of the rock formation
plot(x, z, '-o');
hold on;
grid on;
% Earth's surface
earth_surface = zeros(size(x));

% Plot fault lines and gather fault information
for i = 1:length(fault_indices)
    idx = fault_indices(i);
    plot([x(idx-1), x(idx)], [z(idx-1), z(idx)], 'r--');
    
    % Calculate intersection point with Earth's surface
    intersection_x = x(idx-1) - z(idx-1) * (x(idx) - x(idx-1)) / (z(idx) - z(idx-1));
    plot(intersection_x, 0, 'ro'); % Plot intersection point
    set(gca, 'YDir', 'reverse');
    
    % Calculate dip angle using the tan-rule
    dip_angle = atand((z(idx) - z(idx-1)) / (x(idx) - x(idx-1)));
    
    % Determine fault type
    if dip_angle <= 90
        fault_type = 'Normal';
        plot(intersection_x, 0, 'go'); % Plot normal fault intersection point
    elseif dip_angle > 90
        fault_type = 'Reverse';
        plot(intersection_x, 0, 'mo'); % Plot reverse fault intersection point
    else
        fault_type = 'Strike-Slip';
    end
    
    % Store fault information
    fault_locations = [fault_locations; intersection_x];
    fault_angles = [fault_angles; dip_angle];
    fault_types{end+1} = fault_type;
end

% Display fault information in a table
if isempty(fault_indices)
    
    % Display a message if there are no faults
    disp('No fault lines detected.');
else
    % Display fault information in a table
    T = table(fault_locations, fault_angles, fault_types', ...
        'VariableNames', {'Location', 'Angle', 'Type'});
    disp(T);
end

    Location    Angle        Type   
    ________    ______    __________

       4.9      -78.69    {'Normal'}
         9          90    {'Normal'}
      14.2      69.864    {'Normal'}

Answer 1

[Walter Roberson]

atand() with one parameter returns values in the range -90 to +90, not in the range 0 to 90 to 180. No single-parameter atand() values will be greater than 90

You should consider whether you want to use the two-parameter form of atand().

And you should take into account negative results from atand()

[Moustafa Abedel Fattah]

14.2 69.864 {'Normal'}

should be "Reverse" instead:"Normal"

[Moustafa Abedel Fattah]

Write the code completely please

[Walter Roberson]

should be "Reverse" instead:"Normal"

Not according to your code, it should not be. The atand is 69.864 which is <= 90 and you have defined <= 90 as being Normal.

From the documentation: atand

For real values of X, atand(X) returns values in the interval [-90, 90].

-90 to +90 is <= 90 so for real angles, 'Normal' will always be satisfied.

[Moustafa Abedel Fattah]

You are right Sir