Drawing planes in geology
Dear all
I need help to add an algoritm or modified the exsit one to get the figures annexed in file:the following is my code:
clc; clear; close all;
%% —————- 1. Load Data ————————————
filename = ‘ThreeFaultModel_Modified.xlsx’;
data = xlsread(filename);
x = data(:,1);
g_Bouguer = data(:,2);
x = x(:); g_Bouguer = g_Bouguer(:);
%% —————- 2. User Inputs ———————————
N_layers = input(‘Enter number of layers: ‘);
rho_profiles = zeros(1, N_layers);
for i = 1:N_layers
rho_profiles(i) = input([‘Density of layer ‘ num2str(i) ‘: ‘]);
end
rho_basement = input(‘Basement density: ‘);
rho_contrasts = rho_profiles – rho_basement;
%% —————- 3. Inversion ————————
G = 6.67e-3;
z_inv = bsxfun(@rdivide, g_Bouguer, (2*pi*G*rho_contrasts));
%% —————- 4. Detect ONE Main Fault ————————
VG = gradient(z_inv(:,1), x);
% اختيار أقوى تغير (fault)
[~, idx_fault] = max(abs(VG));
MainFault.x = x(idx_fault);
MainFault.dip = atan2d(VG(idx_fault),1); % تقدير تقريبي
MainFault.throw = max(z_inv(:,1)) * 0.3; % تقدير بسيط
MainFault.type = ‘Normal’;
if MainFault.dip > 0
MainFault.type = ‘Reverse’;
end
if abs(MainFault.dip) > 85
MainFault.type = ‘Vertical’;
end
%% —————- 5. Apply ONE Fault to All Layers —————-
z_layers = z_inv;
for L = 1:N_layers
if strcmp(MainFault.type,’Reverse’)
z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) + MainFault.throw;
elseif strcmp(MainFault.type,’Normal’)
z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) – MainFault.throw;
else
% Vertical → no vertical shift
end
end
%% —————- 6. Plot Depth ————————
figure; hold on;
for i = 1:N_layers
plot(x, z_layers(:,i),’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
grid on;
title(‘Depth with Single Fault’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
legend(arrayfun(@(k) sprintf(‘Layer %d’,k),1:N_layers,’UniformOutput’,false));
%% —————- 7. Plot Single Fault Plane ————————
z_max = max(z_layers(:));
if strcmp(MainFault.type,’Vertical’)
plot([MainFault.x MainFault.x],[0 z_max],’r–‘,’LineWidth’,2);
else
dx = z_max / tand(MainFault.dip);
if MainFault.dip > 0
x2 = MainFault.x – dx;
else
x2 = MainFault.x + dx;
end
plot([MainFault.x x2],[0 z_max],’r–‘,’LineWidth’,2);
end
%% —————- 8. Dip Line from Surface ————————
figure; hold on;
for i = 1:N_layers
plot(x, z_layers(:,i),’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
grid on;
title(‘Single Fault Dip Line’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
z1 = max(z_layers(:,1));
dip = MainFault.dip;
x0 = MainFault.x;
if abs(dip) > 89
x1 = x0;
else
x1 = x0 + z1 / tand(dip);
end
plot([x0 x1],[0 z1],’r–‘,’LineWidth’,2);
text((x0+x1)/2,0.3,’Main Fault’,’Color’,’r’);
%% —————- 9. Geological Section ————————
figure; hold on;
colors = lines(N_layers);
for i = N_layers:-1:1
if i == 1
top = zeros(size(x));
else
top = z_layers(:,i-1);
end
bot = z_layers(:,i);
fill([x;flipud(x)], [top;flipud(bot)], colors(i,:), …
‘EdgeColor’,’k’,’FaceAlpha’,0.6);
end
% رسم الصدع
if strcmp(MainFault.type,’Vertical’)
plot([MainFault.x MainFault.x],[0 z_max],’k–‘,’LineWidth’,2);
else
dx = z_max / tand(MainFault.dip);
if MainFault.dip > 0
x2 = MainFault.x – dx;
else
x2 = MainFault.x + dx;
end
plot([MainFault.x x2],[0 z_max],’k–‘,’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
title(‘Geological Cross Section with Single Fault’);
legend(arrayfun(@(i) sprintf(‘Layer %d’,i),1:N_layers,’UniformOutput’,false));
we use the following paramters:
Enter number of layers: 3
Density of layer 1: 2.20
Density of layer 2: 2.25
Density of layer 3: 2.30
Basement density: 2.67Dear all
I need help to add an algoritm or modified the exsit one to get the figures annexed in file:the following is my code:
clc; clear; close all;
%% —————- 1. Load Data ————————————
filename = ‘ThreeFaultModel_Modified.xlsx’;
data = xlsread(filename);
x = data(:,1);
g_Bouguer = data(:,2);
x = x(:); g_Bouguer = g_Bouguer(:);
%% —————- 2. User Inputs ———————————
N_layers = input(‘Enter number of layers: ‘);
rho_profiles = zeros(1, N_layers);
for i = 1:N_layers
rho_profiles(i) = input([‘Density of layer ‘ num2str(i) ‘: ‘]);
end
rho_basement = input(‘Basement density: ‘);
rho_contrasts = rho_profiles – rho_basement;
%% —————- 3. Inversion ————————
G = 6.67e-3;
z_inv = bsxfun(@rdivide, g_Bouguer, (2*pi*G*rho_contrasts));
%% —————- 4. Detect ONE Main Fault ————————
VG = gradient(z_inv(:,1), x);
% اختيار أقوى تغير (fault)
[~, idx_fault] = max(abs(VG));
MainFault.x = x(idx_fault);
MainFault.dip = atan2d(VG(idx_fault),1); % تقدير تقريبي
MainFault.throw = max(z_inv(:,1)) * 0.3; % تقدير بسيط
MainFault.type = ‘Normal’;
if MainFault.dip > 0
MainFault.type = ‘Reverse’;
end
if abs(MainFault.dip) > 85
MainFault.type = ‘Vertical’;
end
%% —————- 5. Apply ONE Fault to All Layers —————-
z_layers = z_inv;
for L = 1:N_layers
if strcmp(MainFault.type,’Reverse’)
z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) + MainFault.throw;
elseif strcmp(MainFault.type,’Normal’)
z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) – MainFault.throw;
else
% Vertical → no vertical shift
end
end
%% —————- 6. Plot Depth ————————
figure; hold on;
for i = 1:N_layers
plot(x, z_layers(:,i),’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
grid on;
title(‘Depth with Single Fault’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
legend(arrayfun(@(k) sprintf(‘Layer %d’,k),1:N_layers,’UniformOutput’,false));
%% —————- 7. Plot Single Fault Plane ————————
z_max = max(z_layers(:));
if strcmp(MainFault.type,’Vertical’)
plot([MainFault.x MainFault.x],[0 z_max],’r–‘,’LineWidth’,2);
else
dx = z_max / tand(MainFault.dip);
if MainFault.dip > 0
x2 = MainFault.x – dx;
else
x2 = MainFault.x + dx;
end
plot([MainFault.x x2],[0 z_max],’r–‘,’LineWidth’,2);
end
%% —————- 8. Dip Line from Surface ————————
figure; hold on;
for i = 1:N_layers
plot(x, z_layers(:,i),’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
grid on;
title(‘Single Fault Dip Line’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
z1 = max(z_layers(:,1));
dip = MainFault.dip;
x0 = MainFault.x;
if abs(dip) > 89
x1 = x0;
else
x1 = x0 + z1 / tand(dip);
end
plot([x0 x1],[0 z1],’r–‘,’LineWidth’,2);
text((x0+x1)/2,0.3,’Main Fault’,’Color’,’r’);
%% —————- 9. Geological Section ————————
figure; hold on;
colors = lines(N_layers);
for i = N_layers:-1:1
if i == 1
top = zeros(size(x));
else
top = z_layers(:,i-1);
end
bot = z_layers(:,i);
fill([x;flipud(x)], [top;flipud(bot)], colors(i,:), …
‘EdgeColor’,’k’,’FaceAlpha’,0.6);
end
% رسم الصدع
if strcmp(MainFault.type,’Vertical’)
plot([MainFault.x MainFault.x],[0 z_max],’k–‘,’LineWidth’,2);
else
dx = z_max / tand(MainFault.dip);
if MainFault.dip > 0
x2 = MainFault.x – dx;
else
x2 = MainFault.x + dx;
end
plot([MainFault.x x2],[0 z_max],’k–‘,’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
title(‘Geological Cross Section with Single Fault’);
legend(arrayfun(@(i) sprintf(‘Layer %d’,i),1:N_layers,’UniformOutput’,false));
we use the following paramters:
Enter number of layers: 3
Density of layer 1: 2.20
Density of layer 2: 2.25
Density of layer 3: 2.30
Basement density: 2.67 Dear all
I need help to add an algoritm or modified the exsit one to get the figures annexed in file:the following is my code:
clc; clear; close all;
%% —————- 1. Load Data ————————————
filename = ‘ThreeFaultModel_Modified.xlsx’;
data = xlsread(filename);
x = data(:,1);
g_Bouguer = data(:,2);
x = x(:); g_Bouguer = g_Bouguer(:);
%% —————- 2. User Inputs ———————————
N_layers = input(‘Enter number of layers: ‘);
rho_profiles = zeros(1, N_layers);
for i = 1:N_layers
rho_profiles(i) = input([‘Density of layer ‘ num2str(i) ‘: ‘]);
end
rho_basement = input(‘Basement density: ‘);
rho_contrasts = rho_profiles – rho_basement;
%% —————- 3. Inversion ————————
G = 6.67e-3;
z_inv = bsxfun(@rdivide, g_Bouguer, (2*pi*G*rho_contrasts));
%% —————- 4. Detect ONE Main Fault ————————
VG = gradient(z_inv(:,1), x);
% اختيار أقوى تغير (fault)
[~, idx_fault] = max(abs(VG));
MainFault.x = x(idx_fault);
MainFault.dip = atan2d(VG(idx_fault),1); % تقدير تقريبي
MainFault.throw = max(z_inv(:,1)) * 0.3; % تقدير بسيط
MainFault.type = ‘Normal’;
if MainFault.dip > 0
MainFault.type = ‘Reverse’;
end
if abs(MainFault.dip) > 85
MainFault.type = ‘Vertical’;
end
%% —————- 5. Apply ONE Fault to All Layers —————-
z_layers = z_inv;
for L = 1:N_layers
if strcmp(MainFault.type,’Reverse’)
z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) + MainFault.throw;
elseif strcmp(MainFault.type,’Normal’)
z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) – MainFault.throw;
else
% Vertical → no vertical shift
end
end
%% —————- 6. Plot Depth ————————
figure; hold on;
for i = 1:N_layers
plot(x, z_layers(:,i),’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
grid on;
title(‘Depth with Single Fault’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
legend(arrayfun(@(k) sprintf(‘Layer %d’,k),1:N_layers,’UniformOutput’,false));
%% —————- 7. Plot Single Fault Plane ————————
z_max = max(z_layers(:));
if strcmp(MainFault.type,’Vertical’)
plot([MainFault.x MainFault.x],[0 z_max],’r–‘,’LineWidth’,2);
else
dx = z_max / tand(MainFault.dip);
if MainFault.dip > 0
x2 = MainFault.x – dx;
else
x2 = MainFault.x + dx;
end
plot([MainFault.x x2],[0 z_max],’r–‘,’LineWidth’,2);
end
%% —————- 8. Dip Line from Surface ————————
figure; hold on;
for i = 1:N_layers
plot(x, z_layers(:,i),’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
grid on;
title(‘Single Fault Dip Line’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
z1 = max(z_layers(:,1));
dip = MainFault.dip;
x0 = MainFault.x;
if abs(dip) > 89
x1 = x0;
else
x1 = x0 + z1 / tand(dip);
end
plot([x0 x1],[0 z1],’r–‘,’LineWidth’,2);
text((x0+x1)/2,0.3,’Main Fault’,’Color’,’r’);
%% —————- 9. Geological Section ————————
figure; hold on;
colors = lines(N_layers);
for i = N_layers:-1:1
if i == 1
top = zeros(size(x));
else
top = z_layers(:,i-1);
end
bot = z_layers(:,i);
fill([x;flipud(x)], [top;flipud(bot)], colors(i,:), …
‘EdgeColor’,’k’,’FaceAlpha’,0.6);
end
% رسم الصدع
if strcmp(MainFault.type,’Vertical’)
plot([MainFault.x MainFault.x],[0 z_max],’k–‘,’LineWidth’,2);
else
dx = z_max / tand(MainFault.dip);
if MainFault.dip > 0
x2 = MainFault.x – dx;
else
x2 = MainFault.x + dx;
end
plot([MainFault.x x2],[0 z_max],’k–‘,’LineWidth’,2);
end
set(gca,’YDir’,’reverse’);
xlabel(‘Distance (km)’);
ylabel(‘Depth (km)’);
title(‘Geological Cross Section with Single Fault’);
legend(arrayfun(@(i) sprintf(‘Layer %d’,i),1:N_layers,’UniformOutput’,false));
we use the following paramters:
Enter number of layers: 3
Density of layer 1: 2.20
Density of layer 2: 2.25
Density of layer 3: 2.30
Basement density: 2.67 matlab, drawing, plane MATLAB Answers — New Questions
