% 2D Acoustic Wave Equation with Circular Obstacle
% Simulates wave propagation and produces animation

clear all; close all; clc;

% Simulation parameters
Lx = 10;              % Domain size in x (meters)
Ly = 10;              % Domain size in y (meters)
Nx = 300;             % Number of grid points in x
Ny = 300;             % Number of grid points in y
dx = Lx / (Nx - 1);   % Spatial step in x
dy = Ly / (Ny - 1);   % Spatial step in y
c = 1.0;              % Wave speed (m/s)
dt = 0.4 * min(dx, dy) / c;  % Time step (CFL condition)
T = 15;               % Total simulation time (seconds)
Nt = round(T / dt);   % Number of time steps

% Create spatial grid
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);

% Initialize wave field (pressure)
p = zeros(Ny, Nx);        % Current time step
p_old = zeros(Ny, Nx);    % Previous time step
p_new = zeros(Ny, Nx);    % Next time step

% Create circular obstacle
obstacle_center_x = Lx / 2;
obstacle_center_y = Ly / 2;
obstacle_radius = 1.5;
obstacle = sqrt((X - obstacle_center_x).^2 + (Y - obstacle_center_y).^2) <= obstacle_radius;

% Initial condition: Gaussian pulse
source_x = Lx / 4;
source_y = Ly / 2;
sigma = 0.5;
p = exp(-((X - source_x).^2 + (Y - source_y).^2) / (2 * sigma^2));
p(obstacle) = 0;  % Apply obstacle
p_old = p;

% Set up figure for animation
figure('Position', [100, 100, 800, 600]);
colormap jet;

% Stability factor
r = (c * dt / dx)^2;

% Time-stepping loop
for n = 1:Nt
    % Update interior points using finite difference scheme
    % 2D wave equation: ∂²p/∂t² = c²(∂²p/∂x² + ∂²p/∂y²)
    for i = 2:Ny-1
        for j = 2:Nx-1
            if ~obstacle(i, j)  % Only update non-obstacle points
                % Second-order finite difference in space
                laplacian = (p(i+1, j) - 2*p(i, j) + p(i-1, j)) / dy^2 + ...
                           (p(i, j+1) - 2*p(i, j) + p(i, j-1)) / dx^2;

                % Time integration
                p_new(i, j) = 2*p(i, j) - p_old(i, j) + (c*dt)^2 * laplacian;
            end
        end
    end

    % Apply boundary conditions (absorbing/Dirichlet)
    p_new(1, :) = 0;
    p_new(Ny, :) = 0;
    p_new(:, 1) = 0;
    p_new(:, Nx) = 0;

    % Apply obstacle boundary (rigid wall)
    p_new(obstacle) = 0;

    % Update for next time step
    p_old = p;
    p = p_new;

    % Plot every 10 time steps
    if mod(n, 10) == 0
        surf(X, Y, p, 'EdgeColor', 'none');
        view(2);
        shading interp;
        axis equal;
        xlim([0 Lx]);
        ylim([0 Ly]);
        caxis([-0.5 0.5]);

        % Draw obstacle boundary
        hold on;
        theta = linspace(0, 2*pi, 100);
        plot3(obstacle_center_x + obstacle_radius*cos(theta), ...
              obstacle_center_y + obstacle_radius*sin(theta), ...
              ones(size(theta))*0.6, 'k-', 'LineWidth', 2);
        hold off;

        title(sprintf('2D Acoustic Wave Propagation (t = %.2f s)', n*dt));
        xlabel('x (m)');
        ylabel('y (m)');
        colorbar;

        drawnow;
    end
end

disp('Simulation complete!');
