%% Modeling Acoustic Wave Behavior in Shale-Sandstone Systems: A 2D Discrete Approach with Real-World Variability

%% Preprogram Cleanup:

clc

close all

clear all

%% Initial Parameters:

Gx = 1000; %X-axis grid size (ft)

Gy = 1000; %Y-axis grid size (ft)

Dx = 1; %Spatial step (ft)

Dt = 1; %Time step (microseconds)

Nt = 900; %Number of time steps

%% Acoustic Properties:

SSD = 2.65; %Mean density of sandstones (g/cm^3)

SD = 2.4; %Mean density of shales (g/cm^3)

VSS = 4; %Mean velocity in sandstone (km/s)

VS = 2.75; %Mean velocity in shale (km/s)

%% Real-World Variability:

SSV = 0.25; %Sandstone density + velocity variation of 25%

SV = 0.25; %Shale density + velocity variation of 25%

%% Creating Density and Velocity Grids:

D_grid = zeros(Gx, Gy); %Density grid (empty)

V_grid = zeros(Gx, Gy); %Velocity grid (empty)

S_depth = round(Gy / 4); %Top 250ft is shale

%% Creating the Loop to Assign Random Variation to Density and Velocity Values in Layers:

for j = 1:Gy

    if j <= S_depth

        D_grid(:, j) = SD * (1 + SSV * (2*rand(Gx,1) - 1)); %Randomly vary shale density

        V_grid(:, j) = VS * (1 + SV * (2*rand(Gx,1) - 1)); %Randomly vary shale velocity

    else

        D_grid(:, j) = SSD * (1 + SSV * (2*rand(Gx,1) - 1)); %Randomly vary sandstone density

        V_grid(:, j) = VSS * (1 + SV * (2*rand(Gx,1) - 1)); %Randomly vary shale density

    end

end

%% Computing the Impedance Grid:

I_grid = D_grid .* V_grid; %Density * Velocity = Impedence

%% Initializing Displacement Matrices:

d_current = zeros(Gx, Gy); %Current displacement

d_previous = zeros(Gx, Gy); %Previous displacement

d_next = zeros(Gx, Gy); %Next displacement

%% Initializing the Initial Excitation:

%Top center, focused downward

source_x = round(Gx / 2); %Center of X-axis

source_y = 1; %Top of the Y-axis

spread = 20; %Width of the pulse

%% Creating the Initial 2D Gaussian Pulse to Simulate Wave Generation:

for dx_offset = -spread:spread %Defines how far left and right the pulse spreads

    for dy_offset = -spread:spread %Defines how far up and down the pulse will extend

        x = source_x + dx_offset; %Calculates the horizontal position in the grid

        y = source_y + dy_offset; %Calculates the vertical position in the grid

        if x > 0 && x <= Gx && y > 0 && y <= Gy %Makes sure the location is inside the boundaries

            distance_squared = dx_offset^2 + dy_offset^2;
            %Calculates the squared radial distance from the center of the pulse
            %r² = dx² + dy².
            %The pulse will decrease in strength with the radial distance
            %from the center

            d_current(x, y) = exp(-distance_squared / 10); 
            %Applies a 2D Gaussian function at the position (x, y)
            %The exponent creates a profile centered at (source_x, source_y)
            %The denominator controls the width of the pulse.

        end

    end
    
end

%% Sensor Creation to Model Entire Depth Range:

sensor_x = round(Gx / 2); %Horizontally centered

sensor_y_start = 1; %Starts at depth of 0

sensor_y_end = Gy; %Stops at depth of 1000

DTC = zeros(1, Gy); %Initializes an array to hold the sonic transit time for each depth (DTC)

dtc_detected = false(1, Gy); %Boolean array tracks if DTC has been detected for a depth

threshold = 0.005; %Detection threshold for the sensor (displacement magnitude)

%% Creating the Time Loop for the Wave Propogation Simulation:

figure;

for t = 1:Nt

%Update wave propogation field using the finite difference method
    
    for i = 2:Gx-1

        for j = 2:Gy-1

%Compute the average impedance over neighboring points

            neighbors = [I_grid(i,j), I_grid(i+1,j), I_grid(i-1,j), I_grid(i,j+1), I_grid(i,j-1), I_grid(i+1,j+1), I_grid(i-1,j-1), I_grid(i+1,j-1), I_grid(i-1,j+1)];

            z_avg = sum(neighbors) / numel(neighbors);
            
%Apply the finite difference method to compute next displacement values for
%horizontal, vertical, and diagonal second differences

            laplacian = (d_current(i+1,j) - 2*d_current(i,j) + d_current(i-1,j)) + (d_current(i,j+1) - 2*d_current(i,j) + d_current(i,j-1)) + 0.25 * ((d_current(i+1,j+1) - d_current(i,j)) + (d_current(i-1,j-1) - d_current(i,j)) + (d_current(i+1,j-1) - d_current(i,j)) + (d_current(i-1,j+1) - d_current(i,j)));

%Update the displacement at the next time step

            d_next(i,j) = 2*d_current(i,j) - d_previous(i,j) + (Dt^2 / Dx^2) * laplacian / z_avg;

        end

    end

%% Time Update to Move to the Next Time Step:

    d_previous = d_current;

    d_current = d_next;

%% Sensor detection to Monitor Displacement at Sensor Location:

    for y = sensor_y_start:sensor_y_end

        val = abs(d_current(sensor_x, y)); %Get the absolute displacement at the sensor location

        if ~dtc_detected(y) && val > threshold

            DTC(y) = t * Dt; %Record DTC when the threshold is first exceeded

            dtc_detected(y) = true; %Mark that DTC has been detected for this depth

        end

    end

%% Wave Propagation Visualization Every 25 Time Steps:

if mod(t, 25) == 0 %Executes this block every 25 time steps

    cla; %Clears the current axes to remove the previous visual

    hold on; %Retains the current plot

    rectangle('Position', [1, 1, Gx, Gy/4], 'FaceColor', [0 0 1],'EdgeColor', 'none','FaceAlpha', 0.3);                  

    rectangle('Position', [1, Gy/4 + 1, Gx, Gy*(3/4)],'FaceColor', [1 1 0],'EdgeColor', 'none','FaceAlpha', 0.3);

    h_img = imagesc(d_current'); %Displays the current wave displacement field

    colormap(gca, hot); %Applies a red-yellow colormap for better visibility

    set(h_img, 'AlphaData', 0.9); %Makes the wavefield semi-transparent so the rectangles remain visible underneath

    title(['Wave Propagation at t = ', num2str(t*Dt, '%.1f'), ' \mus']); %Displays current simulation time in microseconds

    xlabel('Horizontal Position'); %X-axis label

    ylabel('Depth (ft)'); %Y-axis label

    axis equal tight; %Equal aspect ratio

    set(gca, 'YDir', 'reverse'); %Flips the Y-axis so that depth increases downward

    ylim([0 Gy]); %Set Y-axis limits

    colorbar; %Indicates wave amplitude/magnitude

    plot([sensor_x, sensor_x], [0, Gy], 'r', 'LineWidth', 2); %Plots the sensor line in red

    drawnow; %Creates the plot on the screen immediately

    hold off; %Releases the hold on the plot

end

end

%% Plot DTC vs. Depth

depth = 1:Gy; %Creates a depth vector from 1 to the number of vertical points

figure; %Creates the figure

scatter(DTC, depth, 10, DTC, 'filled'); hold on; %Plots DTC vs. depth with color coding based on DTC values

colormap(parula); %Parula colormap (I like it)

colorbar; %Displays a colorbar showing the range of DTC values

set(gca, 'YDir', 'reverse'); %Reverses Y-axis for depth

xlabel('Sonic Transit Time (\mus/ft)'); %Label for X-axis

ylabel('Depth (ft)'); %Label for Y-axis

title('DTC vs. Depth'); %Title for the plot

grid on; %Displays grid lines on the plot

%% Plot Impedance vs. Depth:

figure; %Creates the figure

impedance = mean(I_grid, 1);  % Computes the average impedance horizontally at each depth

plot(impedance, depth, 'LineWidth', 2); %Plots the average impedance as a function of depth

set(gca, 'YDir', 'reverse'); %Reverses the Y-axis for depth

xlabel('Impedance (kg/m^2s)'); %Label for the X-axis

ylabel('Depth (ft)'); %Label for the Y-axis

title('Impedance vs. Depth'); %Title for the plot

grid on; %Displays grid lines on the plot

%% Plot Velocity vs. Depth:

figure; %Creates the figure

velocity = mean(V_grid, 1);  %Computes the average velocity horizontally at each depth

plot(velocity, depth, 'LineWidth', 2); %Plot the average velocity as a function of depth

set(gca, 'YDir', 'reverse'); %Reverses the Y-axis for depth

xlabel('Velocity (m/s)'); %Label for the X-axis

ylabel('Depth (ft)'); %Label for the Y-axis

title('Velocity vs. Depth'); %Title for the plot

grid on; %Displays the grid lines on the plot

%% Plot Density vs. Depth:

figure; %Creates the figure

density = mean(D_grid, 1); %Computes the average density horizontally at each depth

plot(density, depth, 'LineWidth', 2); %Plots the average density as a function of depth

set(gca, 'YDir', 'reverse'); %Reverses the Y-axis for depth

xlabel('Density (kg/m^3)'); %Label for the X-axis

ylabel('Depth (ft)'); %Label for the Y-axis

title('Density vs. Depth'); %Title for the plot

grid on; %Displays grid lines for the plot

%% End of Program