%==========================================================================
% Simple_Standing_Wave.m
%
% Study of the acoustical physics associated with the linear superposition
% of two counter-propagating 1-D monochromatic plane/traveling waves,
% = a standing wave.
%
%==========================================================================
% Author: Prof. Steven Errede  4/5/2014 07:50 hr
%
% Last Update:                 4/7/2014 10:20 hr {SME}
%==========================================================================
clear all;
close all;

npts = 40000;

thetar = zeros(1,npts);

MgPtot = zeros(1,npts);
MgUtot = zeros(1,npts);
MgZtot = zeros(1,npts);
MgItot = zeros(1,npts);

PhPtot = zeros(1,npts);
PhUtot = zeros(1,npts);
PhZtot = zeros(1,npts);
PhItot = zeros(1,npts);

Wapot  = zeros(1,npts);
Wakin  = zeros(1,npts);
Watot  = zeros(1,npts);

rho0 =  1.204; % Density of air @ NTP (kg/m^3)
   c =  344.0; % Longitudinal speed of sound @ NTP (m/s)
  z0 = rho0*c; % Longitudinal specific acoustic impedance of free air (Rayls)
  
  p0 =  1.000; % Over-pressure amplitude (RMS Pascals)
  
% |~R| = |~B|/|~A|
%MgR =  1.000 - 1.0e-7; 
 MgR =  0.100;
%MgR =  0.000 + 1.0e-7; 
 
  delphi0_ba  = 0.0; % = phi0b - phi0a   (degrees)
  delphi0_bar = (pi/180.0)*delphi0_ba; % (radians)
  
  thetar_lo = -2.0*pi + 1.0e-7; % = kx_lo
  thetar_hi =  2.0*pi;          % = kx_hi
  
  dthetar = (thetar_hi - thetar_lo)/npts;
  
  for j = 1:npts;
      thetar(j) = thetar_lo + (j-1)*dthetar; % = kx
      
      MgPtot(j) =  p0          *sqrt(1.0 + (2.0*MgR*cos(2.0*thetar(j) + delphi0_bar)) + (MgR*MgR));
      MgUtot(j) = (p0/(rho0*c))*sqrt(1.0 - (2.0*MgR*cos(2.0*thetar(j) + delphi0_bar)) + (MgR*MgR));
      MgZtot(j) = MgPtot(j)/MgUtot(j);
      MgItot(j) = MgPtot(j)*MgUtot(j)/2.0;
      
      Pnum      = sin(thetar(j))*(1.0 + (MgR*cos(2.0*thetar(j) + delphi0_bar))) + cos(thetar(j))*(MgR*sin(2.0*thetar(j) + delphi0_bar));
      Pden      = cos(thetar(j))*(1.0 + (MgR*cos(2.0*thetar(j) + delphi0_bar))) - sin(thetar(j))*(MgR*sin(2.0*thetar(j) + delphi0_bar));
      PhPtot(j) = (180.0/pi)*atan2(Pnum,Pden);
      
      Unum      = sin(thetar(j))*(1.0 - (MgR*cos(2.0*thetar(j) + delphi0_bar))) - cos(thetar(j))*(MgR*sin(2.0*thetar(j) + delphi0_bar));
      Uden      = cos(thetar(j))*(1.0 - (MgR*cos(2.0*thetar(j) + delphi0_bar))) + sin(thetar(j))*(MgR*sin(2.0*thetar(j) + delphi0_bar));
      PhUtot(j) = (180.0/pi)*atan2(Unum,Uden);
      
      Znum      = (2.0*MgR*sin(2.0*thetar(j) + delphi0_bar));
      Zden      = (1.0 - (MgR*MgR));
      PhZtot(j) = (180.0/pi)*atan2(Znum,Zden);
      
      Inum      = (2.0*MgR*sin(2.0*thetar(j) + delphi0_bar));
      Iden      = (1.0 - (MgR*MgR));
      PhItot(j) = (180.0/pi)*atan2(Inum,Iden);
      
       Wapot(j) = 0.5     *(MgPtot(j)*MgPtot(j))/(z0*c);
       Wakin(j) = 0.5*rho0*(MgUtot(j)*MgUtot(j));
       Watot(j) = Wapot(j) + Wakin(j);
  end
  
 figure (01);
  subplot(2,2,1);
  plot (thetar,MgPtot,'m',thetar,MgPtot,'m.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('|p_{tot}({\theta})| (RMS Pascals)');
   title ('|p_{tot}({\theta})| vs. {\theta}');
   
  subplot(2,2,2);
  plot (thetar,MgUtot,'g.',thetar,MgUtot,'g.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('|u_{tot}({\theta})| (RMS m/s)');
   title ('|u_{tot}({\theta})| vs. {\theta}');
   
  subplot(2,2,3);
  %plot   (thetar,MgZtot,'b',thetar,MgZtot,'b.');
  semilogy(thetar,MgZtot,'b',thetar,MgZtot,'b.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('|Z_{tot}({\theta})| (Rayls)');
   title ('|Z_{tot}({\theta})| vs. {\theta}');
   
  subplot(2,2,4);
  plot (thetar,MgItot,'k',thetar,MgItot,'k.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('|I_{tot}({\theta})| (RMS Watts)');
   title ('|I_{tot}({\theta})| vs. {\theta}');

   
 figure (02);
  subplot(2,2,1);
  plot (thetar,PhPtot,'m',thetar,PhPtot,'m.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('{\phi}_{ptot}({\theta}) (degrees)');
   title ('{\phi}_{ptot}({\theta}) vs. {\theta}');
   
  subplot(2,2,2);
  plot (thetar,PhUtot,'g',thetar,PhUtot,'g.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('{\phi}_{utot}({\theta}) (degrees)');
   title ('{\phi}_{utot}({\theta}) vs. {\theta}');
   
  subplot(2,2,3);
  plot    (thetar,PhZtot,'b',thetar,PhZtot,'b.');
 %semilogy(thetar,PhZtot,'b',thetar,PhZtot,'b.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('{\phi}_{Ztot}({\theta}) (degrees)');
   title ('{\phi}_{Ztot}({\theta}) vs. {\theta}');
   
  subplot(2,2,4);
  plot (thetar,PhItot,'k',thetar,PhItot,'k.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('{\phi}_{Itot}({\theta}) (degrees)');
   title ('{\phi}_{Itot}({\theta}) vs. {\theta}');

   
 figure (03);
  plot (thetar,Wapot,'m.',thetar,Wakin,'g.',thetar,Watot,'b.');
  axis tight;
  grid on;
  xlabel ('{\theta} = kx (radians)');
  ylabel ('w_{a}({\theta}) (RMS Joules/m^{3})');
   title ('w_{a}({\theta}) vs. {\theta}');
  legend ('w_{a}^{potl}','W_{a}^{kin}','W_{a}^{tot}');