Need to put symbols on curve

I write the following Matlab code (R2014b) for detecting the geological fault (break down) and the iclination angle and type of fault and I need to put the symboles of fault as in the secobd figuer of the annexed fike - the RasGhari BouguetProfile_AB.xlsx is anned

Thanks in advance

%===========================================================%
clc
close all;
clear;
workspace;
%===========================================================%
% Constants
G = 0.0067;
delt_rohj = 0.050;

% Load data from Excel file
% data = xlsread('RasGhari BouguetProfile_AB.xlsx');
data = xlsread('RasGhariSepBougtProfiles.xlsx');
xc  = data(:, 1); % x-coordinate
gBt = data(:, 2); % Bouguer gravity anomaly

% Calculate gradient of the gravity anomaly profile
gradient_gBt = abs(gradient(gBt));

% Estimate the gradient threshold automatically
gradient_threshold = median(abs(gradient_gBt)) ; % Can adjust by multiplier as needed

% Detect fault locations
fault_locations = find(abs(gradient_gBt) > gradient_threshold);

% Initialize arrays to store fault angles and types
fault_angles = zeros(size(fault_locations));
fault_types = cell(size(fault_locations));

% Parameters
G = 6.67430e-3; % Gravitational constant
delt_rohj = 0.050; % Density contrast (example value, adjust as needed)

for i = 1:numel(fault_locations)
    % Get indices for calculating fault properties
    j = fault_locations(i);
    if j == 1 || j == numel(xc)
        continue;       % Skip if fault is at the boundaries
    end
    
    % Calculate fault angle
    C = abs(atan((gBt(j+1) - gBt(j-1)) / (xc(j+1) - xc(j-1)))); % Angle in Radians 
    B = rad2deg(C);      % Angle in Degrees 
    
    % Calculate h (difference in gBt values)
    h = gBt(i+1) - gBt(i);
    
    g_fault_i = 2 * pi * G * delt_rohj * h *...
        (1 + (1/pi) * atan(xc(i)./gBt(i) + cotd(B)) - (1/pi)...
        * atan(xc(i)./gBt(i+1) + cotd(B))); 
    
    % Determine fault type
    if B <= 90;                  
        fault_type = 'Normal';     
    else
        fault_type = 'Reverse';   
    end
    
    % Store fault angle and type
    fault_angles(i) = 90-B;
    fault_types{i} = fault_type;
end

% Display results
disp('Detected Faults:');

Detected Faults:

disp('Location   Angle (degrees)   Type');

Location   Angle (degrees)   Type

for i = 1:numel(fault_locations)
    disp([num2str(xc(fault_locations(i))), '   ', num2str(fault_angles(i)), '   ', fault_types{i}]);
end

32.5335   87.5575   Normal
34.1602   87.4417   Normal
35.7869   87.3003   Normal
37.4135   86.9421   Normal
39.0402   86.7289   Normal
40.6669   86.5637   Normal
42.2936   85.9853   Normal
43.9202   85.5759   Normal
45.5469   85.5344   Normal
47.1736   85.3   Normal
48.8003   85.1095   Normal
50.4269   85.0995   Normal
52.0536   84.9447   Normal
53.6803   84.8022   Normal
55.307   84.8042   Normal
56.9336   84.6707   Normal
58.5603   84.5332   Normal
60.187   84.5309   Normal
61.8137   84.0933   Normal
63.4403   83.5875   Normal
65.067   83.5192   Normal
66.6937   83.2219   Normal
68.3204   82.8312   Normal
69.947   82.7367   Normal
71.5737   82.5562   Normal
73.2004   82.2813   Normal
74.8271   82.1869   Normal
76.4537   81.9931   Normal
78.0804   81.6392   Normal
79.7071   81.4811   Normal
81.3338   81.6211   Normal
82.9604   81.935   Normal
84.5871   82.1116   Normal
86.2138   82.3994   Normal
87.8405   83.2564   Normal
89.4671   83.8262   Normal
91.0938   83.8981   Normal
92.7205   84.2241   Normal
94.3472   84.4756   Normal
95.9738   84.4911   Normal
97.6005   84.634   Normal
99.2272   84.7591   Normal
100.8539   84.7591   Normal
102.4805   84.8852   Normal
104.1072   85.0095   Normal
105.7339   85.0079   Normal
107.3606   85.0745   Normal
108.9872   85.1483   Normal
110.6139   85.1542   Normal
112.2406   85.1973   Normal
113.8673   85.2511   Normal
115.4939   85.2608   Normal
117.1206   85.2625   Normal
118.7473   85.2651   Normal
120.374   85.2644   Normal
122.0006   85.2911   Normal
123.6273   85.3379   Normal
125.254   85.3566   Normal
126.8807   85.4039   Normal
128.5073   85.5032   Normal
130.134   85.5545   Normal
131.7607   85.5996   Normal
133.3874   85.72   Normal
135.014   85.7957   Normal
136.6407   85.8559   Normal
138.2674   86.0763   Normal
139.8941   86.2382   Normal
141.5207   86.4028   Normal
143.1474   87.3928   Normal
250.508   87.7356   Normal
252.1346   87.722   Normal
253.7613   87.5996   Normal
255.388   87.4512   Normal
257.0147   87.4251   Normal
258.6413   87.3466   Normal
260.268   87.2373   Normal
261.8947   87.2064   Normal
263.5214   87.1593   Normal
265.148   87.0823   Normal
266.7747   87.0521   Normal
268.4014   87.0475   Normal
270.0281   87.0386   Normal
271.6547   87.0335   Normal
273.2814   87.0614   Normal
274.9081   87.1328   Normal
276.5348   87.1752   Normal
278.1615   87.2286   Normal
279.7881   87.4138   Normal
281.4148   87.5449   Normal
283.0415   87.6241   Normal
310.695   86.8596   Normal
312.3216   86.5014   Normal
313.9483   86.029   Normal
315.575   85.1678   Normal
317.2017   84.7776   Normal
318.8283   84.9774   Normal
320.455   85.4346   Normal
322.0817   85.6909   Normal
323.7084   85.9865   Normal
325.335   0   

% Plotting
figure(1)
plot(xc, gBt, 'LineWidth', 1); 
hold on;
grid on;
title('Bouguer Anomaly'); 
xlabel('Profile Points-xc'); 
ylabel('Gravity Anomaly');
legend('Bouguer Anomaly');

figure(2)
plot(xc, gBt, 'LineWidth', 1); 
hold on;
grid on;
title('Bouguer Anomaly'); 
xlabel('Profile Points-xc'); 
ylabel('Gravity Anomaly');

% Add symbols for fault angles
for i = 1:numel(fault_locations)
    x_fault = xc(fault_locations(i));
    y_fault = gBt(fault_locations(i));
    angle = fault_angles(i);
    if strcmp(fault_types{i}, 'Normal')
        symbol = 'o'; % Use circle for normal faults
    else
        symbol = 'x'; % Use cross for reverse faults
    end
    scatter(x_fault, y_fault, 100, 'Marker', symbol, 'MarkerEdgeColor', 'r', 'LineWidth', 2);
    % text(x_fault, y_fault, sprintf('%0.1f°', angle), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
end

legend('Bouguer Anomaly', 'Faults', 'Location', 'best');

[Moustafa Abedel Fattah]

I for got the files

Answer 1

[Mathieu NOE]

Like that ?

%===========================================================%
clc
close all;
clear;
workspace;
%===========================================================%
% Constants
G = 0.0067;
delt_rohj = 0.050;
% Load data from Excel file
% data = xlsread('RasGhari BouguetProfile_AB.xlsx');
data = xlsread('RasGhariSepBougtProfiles.xlsx');
xc  = data(:, 1); % x-coordinate
gBt = data(:, 2); % Bouguer gravity anomaly
% Calculate gradient of the gravity anomaly profile
gradient_gBt = abs(gradient(gBt));
% Estimate the gradient threshold automatically
gradient_threshold = median(abs(gradient_gBt)) ; % Can adjust by multiplier as needed
% Detect fault locations
fault_locations = find(abs(gradient_gBt) > gradient_threshold);
% Initialize arrays to store fault angles and types
fault_angles = zeros(size(fault_locations));
fault_types = cell(size(fault_locations));
% Parameters
G = 6.67430e-3; % Gravitational constant
delt_rohj = 0.050; % Density contrast (example value, adjust as needed)
for i = 1:numel(fault_locations)
    % Get indices for calculating fault properties
    j = fault_locations(i);
    if j == 1 || j == numel(xc)
        continue;       % Skip if fault is at the boundaries
    end
    
    % Calculate fault angle
    C = abs(atan((gBt(j+1) - gBt(j-1)) / (xc(j+1) - xc(j-1)))); % Angle in Radians 
    B = rad2deg(C);      % Angle in Degrees 
    
    % Calculate h (difference in gBt values)
    h = gBt(i+1) - gBt(i);
    
    g_fault_i = 2 * pi * G * delt_rohj * h *...
        (1 + (1/pi) * atan(xc(i)./gBt(i) + cotd(B)) - (1/pi)...
        * atan(xc(i)./gBt(i+1) + cotd(B))); 
    
    % Determine fault type
    if B <= 90             
        fault_type = 'Normal';     
    else
        fault_type = 'Reverse';   
    end
    
    % Store fault angle and type
    fault_angles(i) = 90-B;
    fault_types{i} = fault_type;
end
% Display results
disp('Detected Faults:');
disp('Location   Angle (degrees)   Type');
for i = 1:numel(fault_locations)
    disp([num2str(xc(fault_locations(i))), '   ', num2str(fault_angles(i)), '   ', fault_types{i}]);
end
% Plotting
figure(1)
plot(xc, gBt, 'LineWidth', 1); 
hold on;
grid on;
title('Bouguer Anomaly'); 
xlabel('Profile Points-xc'); 
ylabel('Gravity Anomaly');
legend('Bouguer Anomaly');

figure(2)
plot(xc, gBt, 'LineWidth', 1); 
hold on;
grid on;
title('Bouguer Anomaly'); 
xlabel('Profile Points-xc'); 
ylabel('Gravity Anomaly');
% Add symbols for fault angles
for i = 1:numel(fault_locations)
    x_fault = xc(fault_locations(i));
    y_fault = gBt(fault_locations(i));
    angle = fault_angles(i);
    if strcmp(fault_types{i}, 'Normal')
        symbol = 'o'; % Use circle for normal faults
    else
        symbol = 'x'; % Use cross for reverse faults
    end
    scatter(x_fault, y_fault, 100, 'Marker', symbol, 'MarkerEdgeColor', 'r', 'LineWidth', 2);
    % text(x_fault, y_fault, sprintf('%0.1f°', angle), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
end

% start and end indexes of fault locations 
ind_end = find(diff(fault_locations)>1);
ind_start = [1; ind_end+1];
ind_all = [ind_start;ind_end];
for k = 1:numel(ind_all)
    xx = xc(fault_locations(ind_all(k)));
    yy = gBt(fault_locations(ind_all(k)));
    plot([xx xx],[yy-2 yy+2],'-k','linewidth',2)
end

legend('Bouguer Anomaly', 'Faults', 'Location', 'best');

[Moustafa Abedel Fattah]

Thank you for quick response

I need an automatic way to detect the five-mark positions automaticaly depend gBt data - can you help me to do this ?

[Mathieu NOE]

Sure , but isn't it already doing it now? What else do you need?

[Moustafa Abedel Fattah]

Thank yo I mean fha we needint the red symboles its enough the five-vertical line as indicated for the location of fault. but I need to located it automatical if I load another data of (xc, gBt)

Thanks for your accepted answer

[Mathieu NOE]

hello again

ok, so what is the problem with another data set ? the code would need to be tweaked or ?

What we are looking at are segments with a slope above a treshold value , if I understand correctly

do you fear that this approach would not be robust enough with another data set ?