Half-Basin of Sedimentary Deposition

I wrote the following MATLAB code (Version R2014b) to represnt a half-basin of sedimentray deposition of seven- rock formations overlying basment rocks

%=========================================%

clc

close all;

clear;

workspace;

% Constants

G = 0.0067; % Gravitational constant

pi_value = 22/7; % Value of pi

% Define densities of rock formations and basement

rho = [2.460, 2.290, 2.500, 2.300, 2.400, 2.500, 2.570];

rho_b= 2.680; % basement rock

% Number of columns

M = length(rho);

% Depth of each rock formation in each column

z = 1:M;

depth_baement=3.700;

z = [0.401, 0.654, 0.896, 1.274, 1.371, 2.503, 2.704];

%=========================================================================%

% Calculate average density contrasts for each column

rho_av = zeros(1, M);

for i = 1:M

rho_av(i) = (sum(rho(1:i)) / i) - rho_b;

end

% Calculate gravity effect for each column along the profile

g_z = pi_value * G * rho_av(i) .* z;

% Plot gravity effect curve

figure;

subplot(2,1,1);

plot(z, g_z, 'b', 'LineWidth', 2);

xlabel('Depth of rock formation');

ylabel('Gravity effect (g_z)');

title('Gravity Effect along the Profile');

grid on

% Generate colors for rock formations

colors = parula(M);

%=========================================================================%

% Plot juxtaposing columns

subplot(2,1,2);

hold on;

for i = 1:M

% Stacking rock formations within each column

y_bottom = sum(z(1:i-1)); % Calculate the bottom y-coordinate of the current column

y_top = y_bottom + z(i); % Calculate the top y-coordinate of the current column

rectangle('Position', [i-0.5, y_bottom, 1, z(i)], 'FaceColor', colors(i,:), 'EdgeColor', 'none');

% Add text labels indicating rock formation numbers enclosed within each column

text(i, y_bottom + z(i)/2, num2str(i), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Color', 'k');

end

xlabel('Column');

ylabel('Depth');

title('Juxtaposing Columns of the Basin');

xlim([0.5, M + 0.5]);

ylim([0, max(sum(z)) + depth_baement]); % Adjust ylim to accommodate basement rock

% Create custom color bar

colormap(colors);

cb = colorbar;

cb.Location = 'eastoutside';

cb.Ticks = linspace(0,1,M);

cb.TickLabels = 1:M;

cb.Label.String = 'Rock Formation';

cb.Label.FontSize = 12;

hold off;

% Reverse the vertical scale

set(gca, 'YDir', 'reverse');

% Super title

suptitle('Gravity Effect and Juxtaposing Columns');

%==========================================%

I need help of you to modefied the above MATLAB code to give me the results as in the annexed file (as possible as)

Answer 1

[Abhinaya Kennedy]

1. Adjust the rho and z arrays to match the values from your annexed file.
2. Correct the calculation of g_z to ensure it reflects the gravity effect accurately.
3. Update the plotting commands to match the visual style and labels shown in your annexed file.

clc;
close all;
clear;
workspace;
% Constants
G = 0.0067; % Gravitational constant
pi_value = 22/7; % Value of pi
% Define densities of rock formations and basement
rho = [2.460, 2.290, 2.500, 2.300, 2.400, 2.500, 2.570];
rho_b = 2.680; % basement rock
% Number of columns
M = length(rho);
% Depth of each rock formation in each column
z = [0.401, 0.654, 0.896, 1.274, 1.371, 2.503, 2.704];
depth_basement = 3.700;

% Calculate average density contrasts for each column
rho_av = zeros(1, M);
for i = 1:M
    rho_av(i) = (sum(rho(1:i)) / i) - rho_b;
end
% Calculate gravity effect for each column along the profile
g_z = pi_value * G * rho_av .* z;
% Plot gravity effect curve
figure;
subplot(2,1,1);
plot(z, g_z, 'b', 'LineWidth', 2);
xlabel('Depth of rock formation');
ylabel('Gravity effect (g_z)');
title('Gravity Effect along the Profile');
grid on;
% Generate colors for rock formations
colors = parula(M);

% Plot juxtaposing columns
subplot(2,1,2);
hold on;
for i = 1:M
    % Stacking rock formations within each column
    y_bottom = sum(z(1:i-1)); % Calculate the bottom y-coordinate of the current column
    y_top = y_bottom + z(i);  % Calculate the top y-coordinate of the current column
    rectangle('Position', [i-0.5, y_bottom, 1, z(i)], 'FaceColor', colors(i,:), 'EdgeColor', 'none');
    
    % Add text labels indicating rock formation numbers enclosed within each column
    text(i, y_bottom + z(i)/2, num2str(i), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Color', 'k');
end
xlabel('Column');
ylabel('Depth');
title('Juxtaposing Columns of the Basin');
xlim([0.5, M + 0.5]);
ylim([0, max(sum(z)) + depth_basement]);  % Adjust ylim to accommodate basement rock
% Create custom color bar
colormap(colors);
cb = colorbar;
cb.Location = 'eastoutside';
cb.Ticks = linspace(0,1,M);
cb.TickLabels = 1:M;
cb.Label.String = 'Rock Formation';
cb.Label.FontSize = 12;
hold off;
% Reverse the vertical scale
set(gca, 'YDir', 'reverse');
% Super title
subtitle('Gravity Effect and Juxtaposing Columns');