%==========================================================================
% Acoustic_Monopole_F_Theory.m
%
% Written by Prof. Steven Errede   April  3, 2012  11:40 hr
%
% Last Updated:                    April 22, 2016  09:20 hr {SME}
% 
%==========================================================================
% We calculate the following frequency domain complex relations vs. r:
% =====================================================================
% P(r) = ReP(r) + i*ImP(r), |P(r)| = sqrt(ReP(r)^2 + ImP(r)^2) {Pascals}
% U(r) = ReU(r) + i*ImU(r), |U(r)| = sqrt(ReU(r)^2 + ImU(r)^2) {m/sec}
% Z(r) = ReZ(r) + i*ImZ(r), |Z(r)| = sqrt(ReZ(r)^2 + ImZ(r)^2) {Ac. Ohms}
% I(r) = ReI(r) + i*ImI(r), |I(r)| = sqrt(ReI(r)^2 + ImI(r)^2) {Watts/m^2}
%
% and:
%
% Complex particle displacement:         D(r) =   U(r)/iw,   w = 2pi*f
% Complex particle acceleration:         A(r) = iwU(r)
%
% and:
%
% Complex sound intensity relations:     I(r) = P(r)*U^*(r)
% Complex version of Ohm's Law:          Z(r) = P(r)/U(r)   = rho0*V(r)
% Complex linear  momentum density:      G(r) = I(r)/c0^2
% Complex acoustic energy flow velocity: V(r) = Z(r)/rho0   =   c0*N(r)
% Complex "index of refraction":         N(r) = V(r)/c0     =   Z(r)/Z0
% Complex wavelength:                    L(r) = V(r)/f
% Complex wavenumber:                    K(r) =  2pi/L(r)
% Purely real acoustic energy density:   W(r) = (rho0/2)|P(r)|^2/z0^2 + (rho0/2)|U(r)|^2
%
%==========================================================================
close all; % close out all previously displayed plots.
clear all; % clear/zero-out/initialize all variables & arrays.

%===================================================================
% Free-Air Parameters:
%===================================================================
rho0 =  1.2040; % mass density of air (kg/m^3) @ NTP (20 deg. C, 1 Atm)
  c0 =  345.00; % nominal free-field/free-air speed of sound @ NTP in 6105ESB
  Z0 = rho0*c0; % specific long. acoust. impedance free-air (~ 415 Ohms @ NTP)
 mu0 = 1.85e-5; % dynamic viscosity of air @ NTP (Pascal-sec)
 
Iref = 1.00e-3; % Reference Sound Intensity, Io (RMS nano-Watts/m^2)
Pref = 2.00e-5; % Reference Sound Pressure,  Po (RMS Pascals)
Uref = 4.80e-5; % Reference Sound PartVeloc, Uo (RMS mm/s)

%======================================================================
% Other needed input parameters:
%======================================================================
freq  =  300.0; % Fixed frequency (Hz)

B0    = 1.0000; % "Amplitude" Coefficient (SI Units: RMS Pascal - m)

npts  =  10000; % 0.001 < r < 10.0 m in 0.001 m steps

r_max =   10.0; % controls max radial dist (m) of plots (only)!

% ========================
% preallocate all arrays:
% ========================
rps    = zeros(1,npts); % radial distance from source (m)
Re_k   = zeros(1,npts);

 ReP   = zeros(1,npts);
 ImP   = zeros(1,npts);
 MgP   = zeros(1,npts);
 PhP   = zeros(1,npts);
 IFP   = zeros(1,npts);
 
 ReU   = zeros(1,npts);
 ImU   = zeros(1,npts);
 MgU   = zeros(1,npts);
 PhU   = zeros(1,npts);
 IFU   = zeros(1,npts);
 
 ReD   = zeros(1,npts);
 ImD   = zeros(1,npts);
 MgD   = zeros(1,npts);
 PhD   = zeros(1,npts);
 IFD   = zeros(1,npts);
 
 ReA   = zeros(1,npts);
 ImA   = zeros(1,npts);
 MgA   = zeros(1,npts);
 PhA   = zeros(1,npts);
 IFA   = zeros(1,npts);
 
 ReZ   = zeros(1,npts);
 ImZ   = zeros(1,npts);
 MgZ   = zeros(1,npts);
 PhZ   = zeros(1,npts);
 IFZ   = zeros(1,npts);
 
 ReI   = zeros(1,npts);
 ImI   = zeros(1,npts);
 MgI   = zeros(1,npts);
 PhI   = zeros(1,npts);
 IFI   = zeros(1,npts);
  rI   = zeros(1,npts);
 
 ReG   = zeros(1,npts);
 ImG   = zeros(1,npts);
 MgG   = zeros(1,npts);
 PhG   = zeros(1,npts);
 IFG   = zeros(1,npts);
 
 ReV   = zeros(1,npts);
 ImV   = zeros(1,npts);
 MgV   = zeros(1,npts);
 PhV   = zeros(1,npts);
 IFV   = zeros(1,npts);
 
 ReN   = zeros(1,npts);
 ImN   = zeros(1,npts);
 MgN   = zeros(1,npts);
 PhN   = zeros(1,npts);
 IFN   = zeros(1,npts);
 
 ReK   = zeros(1,npts);
 ImK   = zeros(1,npts);
 MgK   = zeros(1,npts);
 PhK   = zeros(1,npts);
 IFK   = zeros(1,npts);
 
 ReL   = zeros(1,npts);
 ImL   = zeros(1,npts);
 MgL   = zeros(1,npts);
 PhL   = zeros(1,npts);
 IFL   = zeros(1,npts);
 
 SIL   = zeros(1,npts);
 SILr  = zeros(1,npts);
 SILi  = zeros(1,npts);
dSIL   = zeros(1,npts);
rSIL   = zeros(1,npts);

 SPL   = zeros(1,npts);
 SUL   = zeros(1,npts);
  
dSLip  = zeros(1,npts);
dSLiu  = zeros(1,npts);
dSLpu  = zeros(1,npts);

  Wpot = zeros(1,npts);
  Wkin = zeros(1,npts);
  Wtot = zeros(1,npts);
  Rwpk = zeros(1,npts);
  Fpot = zeros(1,npts);
  Fkin = zeros(1,npts);
  Ftot = zeros(1,npts);
  
  Wrad = zeros(1,npts);
  Wvrt = zeros(1,npts);
  Rwrv = zeros(1,npts);
  Frad = zeros(1,npts);
  Fvrt = zeros(1,npts);
  
for j=1:npts; 
   rps(j) = (0.001*(j-1)) + 0.010; % 0.010 < r < 10.0 m in 0.001 m steps
   
  Re_k(j) = (2.0*pi*freq)/c0; % for RH ~ 50% air @ NTP... (m^-1)
end

%====================================================================
% Calculate complex pressure & particle velocity @ +z at time t = 0:
%====================================================================
for j=1:npts;
    ReP(j) =  (B0/rps(j))*cos(-Re_k(j)*rps(j));
    ImP(j) =  (B0/rps(j))*sin(-Re_k(j)*rps(j));

 % Compute MgP = |P|, PhP = PhaseP = ArcTan(ImP/ReP) and Cos(PhP):
    MgP(j) =     sqrt((ReP(j)*ReP(j)) + (ImP(j)*ImP(j)));
    PhP(j) = (180.0/pi)*atan2(ImP(j),ReP(j));
    IFP(j) =        cos(atan2(ImP(j),ReP(j)));
  
 %fprintf('%i %f %f %f %f %f %f \n',j,rps(j),ReP(j),ImP(j),MgP(j),PhP(j),IFP(j));

 %===================================================
 % Euler eqn prediction for complex U (mm/s):
 %===================================================
    ReU(j) =  1000.0*(B0/(Z0*rps(j)))*(cos(-Re_k(j)*rps(j)) + (1.0/(Re_k(j)*rps(j)))*sin(-Re_k(j)*rps(j))); % n.b. mm/s !!!
    ImU(j) =  1000.0*(B0/(Z0*rps(j)))*(sin(-Re_k(j)*rps(j)) - (1.0/(Re_k(j)*rps(j)))*cos(-Re_k(j)*rps(j))); % n.b. mm/s !!!
  
 % Compute MgU = |U| and PhU = PhaseU = ArcTan(ImU/ReU):
    MgU(j) =     sqrt((ReU(j)*ReU(j)) + (ImU(j)*ImU(j)));
    PhU(j) = (180.0/pi)*atan2(ImU(j),ReU(j));
    IFU(j) =        cos(atan2(ImU(j),ReU(j)));
  
 %fprintf('%i %f %f %f %f %f %f \n',j,rps(j),ReU(j),ImU(j),MgU(j),PhU(j),IFU(j));
end

    PhP    = (180.0/pi)*unwrap((pi/180.0)*PhP); % Unwrap P-phase
    PhU    = (180.0/pi)*unwrap((pi/180.0)*PhU); % Unwrap U-phase
    
%==========================================================================
% Compute complex longitudinal particle displacement D and acceleration A:
%==========================================================================
for j=1:npts; 
    ReD(j) =     ImU(j)/(2*pi*freq); % Units: RMS mm
    ImD(j) =    -ReU(j)/(2*pi*freq); % Units: RMS mm
      
    MgD(j) =     sqrt((ReD(j)*ReD(j)) + (ImD(j)*ImD(j)));
    PhD(j) = (180.0/pi)*atan2(ImD(j),ReD(j));
    IFD(j) =        cos(atan2(ImD(j),ReD(j)));

    ReA(j) =    -ImU(j)*(2*pi*freq); % Units: RMS mm/s^2
    ImA(j) =     ReU(j)*(2*pi*freq); % Units: RMS mm/s^2
      
    MgA(j) =     sqrt((ReA(j)*ReA(j)) + (ImA(j)*ImA(j)));
    PhA(j) = (180.0/pi)*atan2(ImA(j),ReA(j));
    IFA(j) =        cos(atan2(ImA(j),ReA(j)));
end

    PhD    = (180.0/pi)*unwrap((pi/180.0)*PhD); % Unwrap D-phase
    PhA    = (180.0/pi)*unwrap((pi/180.0)*PhA); % Unwrap A-phase

%==========================================================================
% Calculate complex longitudinal specific acoustic impedance:
%==========================================================================
for j=1:npts; 
  % n.b. don't forget that U is in mm/s - convert back to m/sec!
    ReZ(j) = 1000.0*((ReP(j)*ReU(j)) + (ImP(j)*ImU(j)))/(MgU(j)*MgU(j));
    ImZ(j) = 1000.0*((ImP(j)*ReU(j)) - (ReP(j)*ImU(j)))/(MgU(j)*MgU(j));
  
  % Compute MgZ = |Z|, PhZ = PhaseZ = ArcTan(ImZ/ReZ) and IFZ = Cos(PhZ):
    MgZ(j) =     sqrt((ReZ(j)*ReZ(j)) + (ImZ(j)*ImZ(j)));
    PhZ(j) = (180.0/pi)*atan2(ImZ(j),ReZ(j));
    IFZ(j) =        cos(atan2(ImZ(j),ReZ(j)));
  
  %fprintf('%i %f %f %f %f %f \n',j,rps(j),ReZ(j),ImZ(j),MgZ(j),PhZ(j));
end

%==========================================================================
% Compute complex longitudinal acoustic intensity I = P*U^cc:
%==========================================================================
for j=1:npts; 
  % n.b. don't forget that U in mm/s - convert to m/sec, then express I in RMS nano-Watts/m^2!
    ReI(j) = 1.0E+6*((ReP(j)*ReU(j)) + (ImP(j)*ImU(j)));             % (RMS nano-Watts/m^2)
    ImI(j) = 1.0E+6*((ImP(j)*ReU(j)) - (ReP(j)*ImU(j)));             % (RMS nano-Watts/m^2)
  
  % Compute MgI = |I|, PhI = PhaseI = ArcTan(ImI/ReI) and IFI = Cos(PhI):
    MgI(j) =     sqrt((ReI(j)*ReI(j)) + (ImI(j)*ImI(j)));
    PhI(j) = (180.0/pi)*atan2(ImI(j),ReI(j));
    IFI(j) =        cos(atan2(ImI(j),ReI(j)));
    
     rI(j) =       abs(ImI(j)/ReI(j));

  %fprintf('%i %f %f %f %f %f %f \n',j,rps(j),ReI(j),ImI(j),MgI(j),PhI(j),IFI(j));

   SIL (j) = 10.0*log10(    MgI(j) /Iref);   %   |Sound Intensity Level| (dB)
   SILr(j) = 10.0*log10(abs(ReI(j))/Iref);   % Re{Sound Intensity Level} (dB)
   SILi(j) = 10.0*log10(abs(ImI(j))/Iref);   % Im{Sound Intensity Level} (dB)
  dSIL (j) =       SILr(j)-SILi(j);          % Re-Im SIL Difference      (dB)
  rSIL (j) =       SILi(j)/SILr(j);          % Im/Re SIL Ratio
   
   SPL (j) = 20.0*log10(    MgP(j) /Pref);   %   |Sound Pressure  Level| (dB)
   SUL (j) = 20.0*log10(    MgU(j) /Uref);   %   |Sound PartVeloc Level| (dB)
  
  dSLip(j) = SIL(j) - SPL(j);           % I - P^2 SL Difference (dB)
  dSLiu(j) = SIL(j) - SUL(j);           % I - U^2 SL Difference (dB)
  dSLpu(j) = SPL(j) - SUL(j);           % P^2-U^2 SL Difference (dB)
end

%==========================================================================
% Compute complex wave energy flow velocity V, index of refraction N, 
% wavelength L, wavenumber K, linear momentum density G, angular momentum
% density M & energy density W (purely real):
%==========================================================================
for j = 1:npts;
      ReV(j) =          (1.0/rho0)*ReZ(j); % Units: m/sec
      ImV(j) =          (1.0/rho0)*ImZ(j); % Units: m/sec
      
      MgV(j) =        sqrt((ReV(j)*ReV(j)) + (ImV(j)*ImV(j)));
      PhV(j) =    (180.0/pi)*atan2(ImV(j),ReV(j));
      IFV(j) =           cos(atan2(ImV(j),ReV(j)));

      ReN(j) =              ReV(j)/c0;     % dimensionless
      ImN(j) =              ImV(j)/c0;     % dimensionless
      
      MgN(j) =        sqrt((ReN(j)*ReN(j)) + (ImN(j)*ImN(j)));
      PhN(j) =    (180.0/pi)*atan2(ImN(j),ReN(j));
      IFN(j) =           cos(atan2(ImN(j),ReN(j)));

      ReL(j) =        (1.0/freq)*ReV(j); % Units: m
      ImL(j) =        (1.0/freq)*ImV(j); % Units: m

      MgL(j) =        sqrt((ReL(j)*ReL(j)) + (ImL(j)*ImL(j)));
      PhL(j) =    (180.0/pi)*atan2(ImL(j),ReL(j));
      IFL(j) =           cos(atan2(ImL(j),ReL(j)));

      ReK(j) =             (2.0*pi*ReL(j))/(MgL(j)*MgL(j)); % Units: m-1
      ImK(j) =            -(2.0*pi*ImL(j))/(MgL(j)*MgL(j)); % Units: m-1

      MgK(j) =        sqrt((ReK(j)*ReK(j)) + (ImK(j)*ImK(j)));
      PhK(j) =    (180.0/pi)*atan2(ImK(j),ReK(j));
      IFK(j) =           cos(atan2(ImK(j),ReK(j)));

      ReG(j) =       (1.0/(c0*c0))*ReI(j); % Units: nano-kg/s-m^2
      ImG(j) =       (1.0/(c0*c0))*ImI(j); % Units: nano-kg/s-m^2
      
      MgG(j) =        sqrt((ReG(j)*ReG(j)) + (ImG(j)*ImG(j)));
      PhG(j) =    (180.0/pi)*atan2(ImG(j),ReG(j));
      IFG(j) =           cos(atan2(ImG(j),ReG(j)));

     Wpot(j) =  1.0e9*rho0*(MgP(j)*MgP(j))/(Z0*Z0); % Units: nJ/m^3
     Wkin(j) =  1.0e3*rho0*(MgU(j)*MgU(j));         % Units: nJ/m^3 (U in mm/s)
     Wtot(j) =            Wpot(j)+Wkin(j);          % Units: nJ/m^3
     Rwpk(j) =            Wpot(j)/Wkin(j);
     
     Fpot(j) =            Wpot(j)/Wtot(j);
     Fkin(j) =            Wkin(j)/Wtot(j);
     Ftot(j) =            Fpot(j)+Fkin(j);

     Frad(j) =    (ReV(j)*ReV(j))/(MgV(j)*MgV(j));
     Fvrt(j) =    (ImV(j)*ImV(j))/(MgV(j)*MgV(j));
     
     Wrad(j) =            Frad(j)*Wtot(j);
     Wvrt(j) =            Fvrt(j)*Wtot(j);
     Rwrv(j) =            Wrad(j)/Wvrt(j);
end

%==========================================================================
% Plot Complex P Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xP(1) =   0.0;
xP(2) =   0.0;

yP(1) =   0.0;
yP(2) =   0.0;

figure(31);

subplot(4,2,1);
 plot(rps,ReP,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re P(r) (RMS Pa)');
  title('Complex P(r) vs. r');

subplot(4,2,2);
 plot(rps,PhP,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi P(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImP,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im P(r) (RMS Pa)');
 
subplot(4,2,4);
 semilogx(rps,PhP,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi P(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgP,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|P(r)| (RMS Pa)');

subplot(4,2,6);
 plot(rps,IFP,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi P(r))');

subplot(4,2,7);
 semilogy(rps,MgP,'b',xg,yg,'k');
 axis tight;
 grid on;
 xlabel('Radial Position, r (m)');
 ylabel('|P(r)| (RMS Pa)');

subplot(4,2,8);
 plot(ReP,ImP,'b',xP,yP,'k');
 axis tight;
 grid    on;
 xlabel('Re P(r) (RMS Pa)');
 ylabel('Im P(r) (RMS Pa)');

%==========================================================================
% Plot Complex U Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xU(1) =   0.0;
xU(2) =   0.0;

yU(1) =   0.0;
yU(2) =   0.0;

figure(32);

subplot(4,2,1);
 plot(rps,ReU,'b',xg,yg,'k');
 axis tight;
 grid on;
 xlabel('Radial Position, r (m)');
 ylabel('Re U(r) (RMS mm/s)');
  title('Complex U(r) vs. r');

subplot(4,2,2);
 plot(rps,PhU,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi U(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImU,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im U(r) (RMS mm/s)');
 
subplot(4,2,4);
 semilogx(rps,PhU,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi U(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgU,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|U(r)| (RMS mm/s)');

subplot(4,2,6);
 plot(rps,IFU,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi U(r))');

subplot(4,2,7);
 semilogy(rps,MgU,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|U(r)| (RMS mm/s)');

subplot(4,2,8);
 plot(ReU,ImU,'b',xU,yU,'k');
 axis tight;
 grid    on;
 xlabel('Re U(r) (RMS mm/s)');
 ylabel('Im U(r) (RMS mm/s)');

%==========================================================================
% Plot Complex Z Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xZ(1) =   0.0;
xZ(2) =   0.0;

yZ(1) =   0.0;
yZ(2) =   0.0;

figure(33);

subplot(4,2,1);
 plot(rps,ReZ,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re Z(r) (Ohms)');
  title('Complex Z(r) vs. r');

subplot(4,2,2);
 plot(rps,PhZ,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi Z(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImZ,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im Z(r) (Ohms)');
 
subplot(4,2,4);
 semilogx(rps,PhZ,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi Z(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgZ,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|Z(r)| (Ohms)');

subplot(4,2,6);
 plot(rps,IFZ,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi Z(r))');

subplot(4,2,7);
 semilogy(rps,MgZ,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|Z(r)| (Ohms)');

subplot(4,2,8);
 plot(ReZ,ImZ,'b',xZ,yZ,'k');
 axis tight;
 grid    on;
 xlabel('Re Z(r) (Ohms)');
 ylabel('Im Z(r) (Ohms)');

%==========================================================================
% Plot complex I Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xI(1) =   0.0;
xI(2) =   0.0;

yI(1) =   0.0;
yI(2) =   0.0;

figure(34);

subplot(4,2,1);
 plot(rps,ReI,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re I(r) (RMS nW/m^{2})');
  title('Complex I(r) vs. r');

subplot(4,2,2);
 plot(rps,PhI,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi I(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImI,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im I(r) (RMS nW/m^{2})');
 
subplot(4,2,4);
 semilogx(rps,PhI,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi I(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgI,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|I(r)| (RMS nW/m^{2})');

subplot(4,2,6);
 plot(rps,IFI,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi I(r))');

subplot(4,2,7);
 semilogy(rps,MgI,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position (m)');
 ylabel('|I(r)| (RMS nW/m^{2})');

subplot(4,2,8);
 plot(ReI,ImI,'b',xI,yI,'k');
 axis tight;
 grid    on;
 xlabel('Re I(r) (RMS nW/m^{2})');
 ylabel('Im I(r) (RMS nW/m^{2})');

%==========================================================================
% Plot SIL, SILr, SILi, SPL, SUL & their differences vs. Position:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xI(1) =  0.0;
xI(2) =  0.0;

yI(1) =  0.0;
yI(2) =  0.0;

figure (35);
subplot(4,2,1);
plot(rps,SIL,'b',xg,yg,'k');
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('SIL(r) (dB)');
 title('SIL(r), SPL(r), SUL(r) vs. r');

subplot(4,2,3);
plot(rps,SPL,'g',xg,yg,'k');
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('SPL(r) (dB)');

subplot(4,2,5);
plot(rps,SUL,'r',xg,yg,'k');
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('SUL(r) (dB)');

subplot(4,2,7);
plot(rps,SIL,'b',rps,SPL,'g',rps,SUL,'r',xg,yg,'k');
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('SIL(r), SPL(r), SUL(r) (dB)');

subplot(4,2,2);
plot(rps,dSLip,'b',xg,yg,'k');
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('dSPip(r) = SIL(r)-SPL(r) (dB)');

subplot(4,2,4);
plot(rps,dSLiu,'g',xg,yg,'k');
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('dSPiu(r) = SIL(r)-SUL(r) (dB)');

subplot(4,2,6);
plot(rps,dSLpu,'r',xg,yg,'k');
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('dSPpu(r) = SPL(r)-SUL(r) (dB)');

subplot(4,2,8);
plot(rps,dSLip,'b',rps,dSLiu,'g',rps,dSLpu,'r',xg,yg,'k');
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('dSPip(r), dSPiu(r), dSPpu(r) (dB)');

figure(36);
semilogx(rps,SIL, 'b','linewidth',2);
hold    on;
semilogx(rps,SILr,'r','linewidth',2);
semilogx(rps,SILi,'g','linewidth',2);
hold   off;
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('SIL(r)/SIL_{r}(r)/SIL_{i}(r) (dB)');
 title('SIL(r), SIL_{r}(r), SIL_{i}(r) vs. r');
legend('SIL','SIL_{r}','SIL_{i}');

xg(1) = min(rps);
xg(2) = max(rps);
yg(1) = 0.0;
yg(2) = 0.0;

figure(37);
semilogx(rps,dSIL,'b','linewidth',2);
hold    on;
semilogx( xg,  yg,'k','linewidth',2);
hold   off;
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('dSIL(r) = SIL_{r}(r) - SIL_{i}(r) (dB)');
 title('dSIL(r) = SIL_{r}(r) - SIL_{i}(r) vs. r');
 
xg(1) = min(rps);
xg(2) = max(rps);
yg(1) = 1.0;
yg(2) = 1.0;

figure(38);
semilogx(rps,rSIL,'b','linewidth',2);
hold    on;
semilogx( xg,  yg,'k','linewidth',2);
hold   off;
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('rSIL(r) = SIL_{r}(r)/SIL_{i}(r)');
 title('Ratio SIL_{r}(r)/SIL_{i}(r) vs. r');

figure(39);
loglog(rps,    MgI ,'b','linewidth',2);
hold    on;
loglog(rps,abs(ReI),'r','linewidth',2);
loglog(rps,abs(ImI),'g','linewidth',2);
hold   off;
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('|Ia(r)|/|Ia_{r}(r)|/|Ia_{i}(r)|');
 title('|Ia(r)|, |Ia_{r}(r)|, |Ia_{i}(r)| vs. r');
legend('|Ia|','|Ia_{r}|','|Ia_{i}|');

xg(1) = min(rps);
xg(2) = max(rps);
yg(1) = 1.0;
yg(2) = 1.0;
 
figure(40);
loglog(rps,rI,'b','linewidth',2);
hold    on;
loglog( xg,yg,'k','linewidth',2);
hold   off;
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('rIa(r) = |Ia_{r}(r)/Ia_{i}(r)|');
 title('Ratio |Ia_{r}(r)/Ia_{i}(r)| vs. r');
 
%==========================================================================
% Plot Complex D Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xD(1) =   0.0;
xD(2) =   0.0;

yD(1) =   0.0;
yD(2) =   0.0;

figure(41);

subplot(4,2,1);
 plot(rps,ReD,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re D(r) (RMS mm)');
  title('Complex D(r) vs. r');

subplot(4,2,2);
 plot(rps,PhD,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi D(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImD,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im D(r) (RMS mm)');
 
subplot(4,2,4);
 semilogx(rps,PhD,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi D(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgD,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|D(r)| (RMS mm)');

subplot(4,2,6);
 plot(rps,IFD,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi D(r))');

subplot(4,2,7);
 semilogy(rps,MgD,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|D(r)| (RMS mm)');

subplot(4,2,8);
 plot(ReD,ImD,'b',xD,yD,'k');
 axis tight;
 grid    on;
 xlabel('Re D(r) (RMS mm)');
 ylabel('Im D(r) (RMS mm)');

%==========================================================================
% Plot Complex A Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xA(1) =   0.0;
xA(2) =   0.0;

yA(1) =   0.0;
yA(2) =   0.0;

figure(42);

subplot(4,2,1);
 plot(rps,ReA,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re A(r) (RMS mm/s^{2})');
  title('Complex A(r) vs. r');

subplot(4,2,2);
 plot(rps,PhA,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi A(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImA,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im A(r) (RMS mm/s^{2})');
 
subplot(4,2,4);
 semilogx(rps,PhA,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi A(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgA,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|A(r)| (RMS mm/s^{2})');

subplot(4,2,6);
 plot(rps,IFA,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi A(r))');

subplot(4,2,7);
 semilogy(rps,MgA,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|A(r)| (RMS mm/s^{2})');

subplot(4,2,8);
 plot(ReA,ImA,'b',xA,yA,'k');
 axis tight;
 grid    on;
 xlabel('Re A(r) (RMS mm/s^{2})');
 ylabel('Im A(r) (RMS mm/s^{2})');



%==========================================================================
% Loglog plots of |P(r)|, |U(r)|, |Z(r)|, |I(r)|, |D(r)| & |A(r)| vs. r
%==========================================================================
figure(51);
subplot(3,2,1);
loglog(rps,MgP,'b')
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('|P(r)| (RMS Pa)');
 title('|P(r)|, |U(r)|, |Z(r)|, |I(r)|, |D(r)| & |A(r)| vs. r');

subplot(3,2,2);
loglog(rps,MgU,'b')
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('|U(r)| (RMS mm/s)');

subplot(3,2,3);
loglog(rps,MgZ,'b')
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('|Z(r)| (Ohms)');

subplot(3,2,4);
loglog(rps,MgI,'b')
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('|I(r)| (RMS nW/m^{2})');

subplot(3,2,5);
loglog(rps,MgD,'b')
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('|D(r)| (RMS mm)');

subplot(3,2,6);
loglog(rps,MgA,'b')
axis tight;
grid    on;
xlabel('Radial Position, r (m)');
ylabel('|A(r)| (RMS mm/s^{2})');

%==========================================================================
%  Draw 3D complex plane for P
%==========================================================================
figure(61);
view(3);
plot3(ReP,ImP,rps,'b');
xlabel('Re(P(r)) [RMS Pa]');
ylabel('Im(P(r)) [RMS Pa]');
zlabel('Radial Position, r (m)');
 title('Complex P(r) vs. r');
axis tight;
grid    on;

%==========================================================================
%  Draw 3D complex plane for U
%==========================================================================
figure(62);
view(3);
plot3(ReU,ImU,rps,'b');
xlabel('Re(U(r)) [RMS mm/s]');
ylabel('Im(U(r)) [RMS mm/s]');
zlabel('Radial Position, r (m)');
 title('Complex U(r) vs. r');
axis tight;
grid    on;

%==========================================================================
%  Draw 3D complex plane for D
%==========================================================================
figure(63);
view(3);
plot3(ReD,ImD,rps,'b');
xlabel('Re(D(r)) [RMS mm]');
ylabel('Im(D(r)) [RMS mm]');
zlabel('Radial Position, r (m)');
 title('Complex D(r) vs. r');
axis tight;
grid    on;

%==========================================================================
%  Draw 3D complex plane for A
%==========================================================================
figure(64);
view(3);
plot3(ReA,ImA,rps,'b');
xlabel('Re(A(r)) [RMS mm/s^{2}]');
ylabel('Im(A(r)) [RMS mm/s^{2}]');
zlabel('Radial Position, r (m)');
 title('Complex A(r) vs. r');
axis tight;
grid    on;

%==========================================================================
%  Draw 3D complex plane for Z
%==========================================================================
figure(65);
view(3);
plot3(ReZ,ImZ,rps,'b');
xlabel('Re(Z(r)) [Ohms]');
ylabel('Im(Z(r)) [Ohms]');
zlabel('Radial Position, r (m)');
 title('Complex Z(r) vs. r');
axis tight;
grid    on;

%==========================================================================
%  Draw 3D complex plane for I
%==========================================================================
figure(66);
view(3);
plot3(ReI,ImI,rps,'b');
xlabel('Re(I(r)) [RMS nW/m^{2}]');
ylabel('Im(I(r)) [RMS nW/m^{2}]');
zlabel('Radial Position, r (m)');
 title('Complex I(r) vs. r');
axis tight;
grid    on;

%==========================================================================
% Plot Complex V Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xV(1) =   0.0;
xV(2) =   0.0;

yV(1) =   0.0;
yV(2) =   0.0;

figure(91);

subplot(4,2,1);
 plot(rps,ReV,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re V(r) (m/s)');
  title('Complex V(r) vs. r');

subplot(4,2,2);
 plot(rps,PhV,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi V(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImV,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im V(r) (m/s)');
 
subplot(4,2,4);
 semilogx(rps,PhV,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi V(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgV,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|V(r)| (m/s)');

subplot(4,2,6);
 plot(rps,IFV,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi V(r))');

subplot(4,2,7);
 semilogy(rps,MgV,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|V(r)| (m/s)');

subplot(4,2,8);
 plot(ReV,ImV,'b',xV,yV,'k');
 axis tight;
 grid    on;
 xlabel('Re V(r) (m/s)');
 ylabel('Im V(r) (m/s)');


%==========================================================================
% Plot Complex N Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xN(1) =   0.0;
xN(2) =   0.0;

yN(1) =   0.0;
yN(2) =   0.0;

figure(92);

subplot(4,2,1);
 plot(rps,ReN,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re N(r)');
  title('Complex N(r) vs. r');

subplot(4,2,2);
 plot(rps,PhN,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi N(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImN,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im N(r)');
 
subplot(4,2,4);
 semilogx(rps,PhN,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi N(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgN,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|N(r)|');

subplot(4,2,6);
 plot(rps,IFN,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi N(r))');

subplot(4,2,7);
 semilogy(rps,MgN,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|N(r)|');

subplot(4,2,8);
 plot(ReN,ImN,'b',xN,yN,'k');
 axis tight;
 grid    on;
 xlabel('Re N(r)');
 ylabel('Im N(r)');

%==========================================================================
% Plot Complex L Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xL(1) =   0.0;
xL(2) =   0.0;

yL(1) =   0.0;
yL(2) =   0.0;

figure(94);

subplot(4,2,1);
 plot(rps,ReL,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re L(r) (m)');
  title('Complex L(r) vs. r');

subplot(4,2,2);
 plot(rps,PhL,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi L(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImL,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im L(r) (m)');
 
subplot(4,2,4);
 semilogx(rps,PhL,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi L(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgL,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|L(r)| (m)');

subplot(4,2,6);
 plot(rps,IFL,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi L(r))');

subplot(4,2,7);
 semilogy(rps,MgL,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|L(r)| (m)');

subplot(4,2,8);
 plot(ReL,ImL,'b',xL,yL,'k');
 axis tight;
 grid    on;
 xlabel('Re L(r) (m)');
 ylabel('Im L(r) (m)');

%==========================================================================
% Plot Complex K Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xK(1) =   0.0;
xK(2) =   0.0;

yK(1) =   0.0;
yK(2) =   0.0;

figure(95);

subplot(4,2,1);
 plot(rps,ReK,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re K(r) (m-1)');
  title('Complex K(r) vs. r');

subplot(4,2,2);
 plot(rps,PhK,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi K(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImK,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im K(r) (m-1)');
 
subplot(4,2,4);
 semilogx(rps,PhK,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi K(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgK,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|K(r)| (m-1)');

subplot(4,2,6);
 plot(rps,IFK,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)')
 ylabel('Cos(Phi K)')

subplot(4,2,7);
 semilogy(rps,MgK,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position (m)');
 ylabel('|K(r)| (m-1)');

subplot(4,2,8);
 plot(ReK,ImK,'b',xK,yK,'k');
 axis tight;
 grid    on;
 xlabel('Re K(r) (m-1)');
 ylabel('Im K(r) (m-1)');

%==========================================================================
% Plot Complex G Data vs. Position - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

xG(1) =   0.0;
xG(2) =   0.0;

yG(1) =   0.0;
yG(2) =   0.0;

figure(96);

subplot(4,2,1);
 plot(rps,ReG,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Re G(r) (RMS nkg/s-m2)');
  title('Complex G(r) vs. r');

subplot(4,2,2);
 plot(rps,PhG,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi G(r) (degrees)');

subplot(4,2,3);
 plot(rps,ImG,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Im G(r) (RMS nkg/s-m2)');
 
subplot(4,2,4);
 semilogx(rps,PhG,'b',10.0,0.0,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Phi G(r) (degrees)');

subplot(4,2,5);
 plot(rps,MgG,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|G(r)| (RMS nkg/s-m2)');

subplot(4,2,6);
 plot(rps,IFG,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Cos(Phi G(r))');

subplot(4,2,7);
 semilogy(rps,MgG,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('|G(r)| (RMS nkg/s-m2)');

subplot(4,2,8);
 plot(ReG,ImG,'b',xG,yG,'k');
 axis tight;
 grid    on;
 xlabel('Re G(r) (RMS nkg/s-m2)');
 ylabel('Im G(r) (RMS nkg/s-m2)');

%==========================================================================
% Plot Wpot, Wkin, Wtot vs. Position - one graphics window, 2x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

figure(111);

subplot(2,2,1);
 semilogy(rps,Wpot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Wpot(r) (RMS nJ/m3)');
  title('Wpot(r), Wkin(r), Wtot(r) vs. r');

subplot(2,2,2);
 semilogy(rps,Wkin,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Wkin(r) (RMS nJ/m3)');

subplot(2,2,3);
 semilogy(rps,Wtot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Wtot(r) (RMS nJ/m3)');

subplot(2,2,4);
 semilogy(rps,Rwpk,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Ratio Wpot(r)/Wkin(r)');

%==========================================================================
% Plot Fpot, Fkin, Ftot vs. Position - one graphics window, 2x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

figure(112);

subplot(2,2,1);
 semilogy(rps,Fpot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Fpot');
  title('Fpot(r), Fkin(r), Ftot(r) vs. r');

subplot(2,2,2);
 semilogy(rps,Fkin,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Fkin(r)');

subplot(2,2,3);
 semilogy(rps,Ftot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Ftot(r)');

subplot(2,2,4);
 semilogy(rps,Rwpk,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Ratio Fpot(r)/Fkin(r)');

%==========================================================================
% Plot Wrad, Wvrt, Wtot vs. Position - one graphics window, 2x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

figure(113);

subplot(2,2,1);
 semilogy(rps,Wrad,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Wrad(r) (RMS nJ/m3)');
  title('Wrad(r), Wvrt(r), Wtot(r) vs. r');

subplot(2,2,2);
 semilogy(rps,Wvrt,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Wvrt(r) (RMS nJ/m3)');

subplot(2,2,3);
 semilogy(rps,Wtot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Wtot(r) (RMS nJ/m3)');

subplot(2,2,4);
 semilogy(rps,Rwrv,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Ratio Wrad(r)/Wvrt(r)');

%==========================================================================
% Plot Frad, Fvrt, Ftot vs. Position - one graphics window, 2x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = r_max;

yg(1) =   0.0;
yg(2) =   0.0;

figure(114);

subplot(2,2,1);
 semilogy(rps,Frad,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Frad(r)');
  title('Frad(r), Fvrt(r), Ftot(r) vs. r');

subplot(2,2,2);
 semilogy(rps,Fvrt,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Fvrt(r)');

subplot(2,2,3);
 semilogy(rps,Ftot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Ftot(r)');

subplot(2,2,4);
 semilogy(rps,Rwrv,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Radial Position, r (m)');
 ylabel('Ratio Frad(r)/Fvrt(r)');

%==========================================================================
 beep;
 
 fprintf('\n Acoustic_Monopole_R_Theory Completed !!! \n')
%==========================================================================