Shifting the parallel lines

3 次查看（过去 30 天）
Moustafa
Moustafa 2025-3-14
评论： Voss 2025-3-14
Ran in:
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'});
  3 个评论
显示 1更早的评论
dpb
dpb 2025-3-14
@Moustafa, you forgot to attach the data file as well...
Moustafa
Moustafa 2025-3-14
移动：Torsten 2025-3-14
Dear Sir,
Sorry for forgetten attached data
please see the attached file

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Voss
Voss 2025-3-14
Ran in:
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'});
  2 个评论
Moustafa
Moustafa 2025-3-14
移动：Voss 2025-3-14
Dear Sir;
Thanks for quick response with accepted solution

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Display and Presentation 的更多信息

标签

产品


版本

R2014b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by