Simulating pressure waves at an impedance boundary

Question

I am trying to simulate pressure waves crossing a boundary from one medium to another (e.g., water to air) in Matlab. The code that I have got so far, which is largely taken from Wikipedia on Partial Transmission only changes the speed as the wave crosses the boundary. The code works using a recursive form or the wave equation found using the method of central differences.

I am confused about how to incorporate the effect of density on the impedance of the medium. From my understanding:

$$ Z=\rho C $$ where $Z$ is the impedance, $\rho$ is the density, $C$ is the speed of sound

At the boundary, the proportion of the waves' amplitude that is reflected and transmitted is given by the reflection and transmission coefficients which are: $$ R=\frac{Z_1-Z_2}{Z_1+Z_2},\\ T=R=\frac{2Z_1}{Z_1+Z_2} $$

At the moment the simulation will correctly calculate the reflected and transmitted waves amplitudes for two mediums with the same density but different wave speeds; however, I have no way of varying the density of the medium.

Seeing as neither density or impedance are relevant in the wave equation I am confused about how they end up determining the reflection and transmission coefficients and how to implements them into my simulation.

Below is the relevant section of my code:

   % length of the string and the grid
   L = 5;
   N = 151;
   X=linspace(0, L, N);

   h = X(2)-X(1); % space grid size
   c = 0.01; % speed of the wave for visualisation 
   tau = 0.25*h/c; % time grid size

   % form a medium with a discontinuous wave speed
   C = 0*X+c;  %this has formed a vector the same dimension as X with every entry =c 

   D=0.5*L;
   c_right = 2*c; % speed to the right of the disc
   for i=1:N
      if X(i) > D
         C(i) = c_right;
      end
   end
   % Now C = c fo x < D, and C=c_right for x > D

   K = 10; % steepness of the bump
   S = 0; % shift the wave
   f=inline('exp(-K*(x-S).^2)', 'x', 'S', 'K'); % a gaussian as an initial wave
   df=inline('-2*K*(x-S).*exp(-K*(x-S).^2)', 'x', 'S', 'K'); % derivative of f

   % wave at time 0 and tau
   U0 = 0*f(X, S, K);
   U1 = U0 - 2*tau*c*df(X, S, K);

   U = 0*U0; % current U

   % plot between Start and End
   Start=0; End=1500;

   for j=1:End

      %  fixed end points
       U(1)=0; U(N)=0;

      % finite difference discretization in time
      for i=2:(N-1)
          %this is the wave equation written algebraically with second
          %order central difference theorem and then rearranged for U(i)
         U(i) = (C(i)*tau/h)^2*(U1(i+1)-2*U1(i)+U1(i-1)) + 2*U1(i) - U0(i);
      end

      % update info, for the next iteration
      U0 = U1; U1 = U;

Answer 1

Essentially going one step back in the derivation of the wave equation you take the divergence of the (simplified) Euler equation $\mathbf{\nabla}.(\rho_{0} \frac{d\mathbf{v}}{dt}=-\nabla p$) and next, you let $\rho_{0}$ out of the operator, assuming it is a constant. Now at the inferface between two media with $\rho_{01} \neq \rho_{02}$, $\nabla \rho_{0}$ is not null anymore, so this simplification is not true. That's the physics part. So you need to implement a specific boundary condition at the interface between two such media in your model to ensure both $p$ and the normal component of particule velocity $v$ (or displacement $u$) are continuous accross the interface.