% 2-D Acoustic FDTD Simulation Demo in JavaScript (rev. Nov. 9th, 2014)
%    by Yoshiki NAGATANI
%      http://ultrasonics.jp/nagatani/e/
%      https://twitter.com/nagataniyoshiki
% 
%   This is a truly physical simulation program of the sound wave propagation
%   in a two-dimensional field filled with fluid media surrounded by total
%   reflecting walls.
% 
%   This program solves the equations of the "2D Acoustic FDTD (finite-difference
%   time-domain) method".
%   The field consists of air (332 m/s). The grid resolution is 10 mm/pixel
%   and the time step is 20 us/step (us = micro second = 1/1000000 second). 
%   A single pulse of sinusoidal sound wave at 1 kHz with Hann window is transmitted.
% 
%   For more detailed information about FDTD method (including 3-D elastic simulation),
%   please refer our papers on simulation. ==> http://ultrasonics.jp/nagatani/fdtd/
%   Thank you.

dx = 0.01;						% Spatical Resolution [m]
dt = 20.0e-6;					% Temporal Resolution [s]

rho = 1.3;						% Density [kg/m^3]
kappa = 142.0e3;				% Volume Elasticily (Bulk modulus) [Pa] */

freq = 1.0e3;					% Frequency of Initial Waveform [Hz]

NX = 300;						% Spatial Size for X direction [pixels]
NY = 400;						% Spatial Size for Y direction [pixels]

Vx = zeros(NX+1,NY  );			% Particle Velocity for X direction [m/s]
Vy = zeros(NX,  NY+1);			% Particle Velocity for Y direction [m/s]
P  = zeros(NX,  NY  );			% Sound Pressure [Pa]

% Window Preparation
h = figure();
image_intensity = 5000;			% Brightness of Acoustic Field (>0)

n = 0;

% Main Loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while true,

	% Update Particle Velocities
	Vx(2:NX,:) = Vx(2:NX,:) - dt / (rho * dx) * ( P(2:NX,:) - P(1:NX-1,:) );
	Vy(:,2:NY) = Vy(:,2:NY) - dt / (rho * dx) * ( P(:,2:NY) - P(:,1:NY-1) );

	% Update Sound Pressure
	P(1:NX,1:NY) = P(1:NX,1:NY) - ( kappa * dt / dx ) * ( ( Vx(2:NX+1,:) - Vx(1:NX,:) ) + ( Vy(:,2:NY+1) - Vy(:,1:NY) ) );

	% Initial Waveform from a Point Source (1 pulse of sinusoidal wave with Hann window)
	if n < (1.0/freq)/dt,
		P(floor(NX/4+1),floor(NY/3+1)) = P(floor(NX/4+1),floor(NY/3+1)) + (1.0-cos(2.0*pi*freq*n*dt))/2.0 * sin(2.0*pi*freq*n*dt);
	end
	n = n + 1;

	% Draw the Acoustic Field
	image( P*image_intensity+32 );
	axis equal; axis tight;
	title(['step: ' sprintf('%5d',n) '   (' sprintf('%6.2f',n*dt*1000) ' ms)']);
	drawnow;

end
