%==========================================================================
% Baffled_Circular_Piston_F_Theory.m
%
% Written by Prof. Steven Errede   April  3, 2012  11:40 hr
%
% Last Updated:                    April 22, 2016  09:25 hr {SME}
% 
%==========================================================================
%
% We calculate the following frequency domain complex relations vs. f:
% =====================================================================
% P(f) = ReP(f) + i*ImP(f), |P(f)| = sqrt(ReP(f)^2 + ImP(f)^2) {Pascals}
% U(f) = ReU(f) + i*ImU(f), |U(f)| = sqrt(ReU(f)^2 + ImU(f)^2) {m/sec}
% Z(f) = ReZ(f) + i*ImZ(f), |Z(f)| = sqrt(ReZ(f)^2 + ImZ(f)^2) {Ac. Ohms}
% I(f) = ReI(f) + i*ImI(f), |I(f)| = sqrt(ReI(f)^2 + ImI(f)^2) {Watts/m^2}
%
% and:
%
% The complex particle displacement:     D(f) =   U(f)/iw,   w = 2pi*f
% The complex particle acceleration:     A(f) = iwU(f)
%
% and:
%
% Complex sound intensity relations:     I(f) = P(f)*U^*(f)
% Complex version of Ohm's Law:          Z(f) = P(f)/U(f)   = rho0*V(f)
% Complex linear  momentum density:      G(f) = I(f)/c0^2
% Complex acoustic energy flow velocity: V(f) = Z(f)/rho0   =   c0*N(f)
% Complex "index of refraction":         N(f) = V(f)/c0     =   Z(f)/Z0
% Complex wavelength:                    L(f) = V(f)/f
% Complex wavenumber:                    K(f) =  2pi/L(f)
% Purely real acoustic energy density:   W(f) = (rho0/2)|P(f)|^2/z0^2 + (rho0/2)|U(f)|^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;   % Referemce Sound Pressure,  Po (RMS Pascals)
Uref = 4.80e-5;   % Referemce Sound PartVeloc, Uo (RMS mm/sec)

%======================================================================
% Dimensions of Circular Piston & other input parameters:
%======================================================================
 Zposn =  0.010;  % Z-position (m) of p/u mics at +Z in front of BB
%Zposn =  1.000;  % Z-position (m) of p/u mics at +Z in front of BB

%Apis  = (0.0254/2.0); % Radius of circular piston (m) (1" dia PZT disk)
%Apis  =  0.1500; % Radius of circular piston (m) (1' dia spkr)
 Apis  =  0.3000; % Radius of circular piston (m) (2' dia spkr)

 Po    =  1.0000; % RMS pressure amplitude (Pascals)
 Uo    =   Po/Z0; % Corresponding RMS particle velocity amplitude

 npts  =   20000; % 29.5 < f < 2030.5 Hz in 0.1 Hz steps

frq_max = 2000.0; % Controls max freq of plots (only)!

% ========================
% preallocate all arrays:
% ========================
 frq   = zeros(1,npts);
 Re_k  = zeros(1,npts);

 arg   = zeros(1,npts);
 ftr   = zeros(1,npts);
 brg   = 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; 
   frq(j) =  0.1*(j - 1) + 29.5;    % 29.5 Hz < f < 2030.5 Hz, in 0.1 Hz steps
   
  Re_k(j) = (2.0*pi*frq(j))/c0; % for RH ~ 50% air @ NTP... (m^-1)
end

%====================================================================
% Calculate complex pressure & particle velocity @ +z at time t = 0:
%====================================================================

for j=1:npts;
    arg(j) =     (Re_k(j)*Zposn);
    ftr(j) =     sqrt(1.0 + ((Apis/Zposn)^2));
    brg(j) = 0.5*(Re_k(j)*Zposn)*(ftr(j) - 1.0);
    
    ReP(j) =  Po*(cos(arg(j)) - cos(arg(j)*ftr(j)));
    ImP(j) = -Po*(sin(arg(j)) - sin(arg(j)*ftr(j)));

 % Compute MgP = |P|, PhP = PhaseP = ArcTan(ImP/ReP) and Cos(PhP):
    MgP(j) =     sqrt((ReP(j)*ReP(j)) + (ImP(j)*ImP(j))); % correct
 %  MgP(j) =   2.0*Po*abs(sin(brg(j)));                   % also correct!
    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,frq(j),ReP(j),ImP(j),MgP(j),PhP(j),IFP(j));

%===================================================
% Euler eqn prediction for complex U (mm/sec):
%===================================================
    ReU(j) =  1000.0*Uo*(cos(arg(j)) - (1.0/ftr(j))*cos(arg(j)*ftr(j)));
    ImU(j) = -1000.0*Uo*(sin(arg(j)) - (1.0/ftr(j))*sin(arg(j)*ftr(j)));
  
 % Compute MgUL = |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,frq(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*frq(j)); % Units: RMS mm
    ImD(j) =    -ReU(j)/(2*pi*frq(j)); % 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*frq(j)); % Units: RMS mm/s^2
    ImA(j) =     ReU(j)*(2*pi*frq(j)); % 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/sec - 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,frq(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/sec - 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(ReI(j)/ImI(j));

  %fprintf('%i %f %f %f %f %f %f \n',j,frq(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) =       SILr(j)/SILi(j);          % Re/Im 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 in/out complex wave velocity V, index of refraction N, 
% propagation constant C, 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/frq(j))*ReV(j); % Units: m
      ImL(j) =        (1.0/frq(j))*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. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReP,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re P (RMS Pa)')
  title('Complex P');

subplot(4,2,2);
 plot(frq,PhP,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi P (degrees)')

subplot(4,2,3);
 plot(frq,ImP,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Im P (RMS Pa)')
 
subplot(4,2,4);
 semilogx(frq,PhP,'b',10.0,0.0,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi P (degrees)')

subplot(4,2,5);
 plot(frq,MgP,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|P| (RMS Pa)')

subplot(4,2,6);
 plot(frq,IFP,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi P)')

subplot(4,2,7);
%semilogy(frq,MgP,'b',xg,yg,'k');
 semilogx(frq,MgP,'b');
%loglog(frq,MgP,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|P| (RMS Pa)')

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

%==========================================================================
% Plot Complex U Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReU,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re U (RMS mm/sec)')
  title('Complex U');

subplot(4,2,2);
 plot(frq,PhU,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi U (degrees)')

subplot(4,2,3);
 plot(frq,ImU,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Im U (RMS mm/sec)')
 
subplot(4,2,4);
 semilogx(frq,PhU,'b',10.0,0.0,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi U (degrees)')

subplot(4,2,5);
 plot(frq,MgU,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|U| (RMS mm/sec)')

subplot(4,2,6);
 plot(frq,IFU,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi U)')

subplot(4,2,7);
%semilogy(frq,MgU,'b',xg,yg,'k');
 semilogx(frq,MgU,'b');
%loglog(frq,MgU,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|U| (RMS mm/sec)')

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

%==========================================================================
% Plot Complex Z Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReZ,'b',xg,yg,'k');
 axis tight 
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re Z (Ohms)')
  title('Complex Z');

subplot(4,2,2);
 plot(frq,PhZ,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi Z (degrees)')

subplot(4,2,3);
 plot(frq,ImZ,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Im Z (Ohms)')
 
subplot(4,2,4);
 semilogx(frq,PhZ,'b',10.0,0.0,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi Z (degrees)')

subplot(4,2,5);
 plot(frq,MgZ,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|Z| (Ohms)')

subplot(4,2,6);
 plot(frq,IFZ,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi Z)')

subplot(4,2,7);
%semilogy(frq,MgZ,'b',xg,yg,'k');
 semilogx(frq,MgZ,'b');
%loglog(frq,MgZ,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|Z| (Ohms)')

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

%==========================================================================
% Plot complex I Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReI,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re I (RMS nW/m^2)')
  title('Complex I');

subplot(4,2,2);
 plot(frq,PhI,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi I (degrees)')

subplot(4,2,3);
 plot(frq,ImI,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Im I (RMS nW/m^2)')
 
subplot(4,2,4);
 semilogx(frq,PhI,'b',10.0,0.0,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi I (degrees)')

subplot(4,2,5);
 plot(frq,MgI,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|I| (RMS nW/m^2)')

subplot(4,2,6);
 plot(frq,IFI,'b',xg,yg,'k');
 axis tight 
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi I)')

subplot(4,2,7);
%semilogy(frq,MgI,'b',xg,yg,'k');
 semilogx(frq,MgI,'b');
%loglog(frq,MgI,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|I| (RMS nW/m^2)')

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

%==========================================================================
% Plot SIL, SPL, SUL & their differences:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,SIL,'b',xg,yg,'k');
axis tight
grid on
xlabel('Frequency (Hz)')
ylabel('SIL (dB)')
 title('SIL, SPL, SUL');

subplot(4,2,3);
plot(frq,SPL,'g',xg,yg,'k');
axis tight
grid on
xlabel('Frequency (Hz)')
ylabel('SPL (dB)')

subplot(4,2,5);
plot(frq,SUL,'r',xg,yg,'k');
axis tight
grid on
xlabel('Frequency (Hz)')
ylabel('SUL (dB)')

subplot(4,2,7);
plot(frq,SIL,'b',frq,SPL,'g',frq,SUL,'r',xg,yg,'k');
axis tight
grid on
xlabel('Frequency (Hz)')
ylabel('SIL, SPL, SUL (dB)')

subplot(4,2,2);
plot(frq,dSLip,'b',xg,yg,'k');
axis tight
grid on
xlabel('Frequency (Hz)')
ylabel('dSPip = SIL-SPL (dB)')

subplot(4,2,4);
plot(frq,dSLiu,'g',xg,yg,'k');
axis tight
grid on
xlabel('Frequency (Hz)')
ylabel('dSPiu = SIL-SUL (dB)')

subplot(4,2,6);
plot(frq,dSLpu,'r',xg,yg,'k');
axis tight
grid on
xlabel('Frequency (Hz)')
ylabel('dSPpu = SPL-SUL (dB)')

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

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

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

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

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

figure(39);
loglog(frq,    MgI, 'b','linewidth',2);
hold    on;
loglog(frq,abs(ReI),'r','linewidth',2);
loglog(frq,abs(ImI),'g','linewidth',2);
hold   off;
axis tight;
grid    on;
xlabel('Frequency, f (Hz)')
ylabel('|Ia(f)|/|Ia_{r}(f)|/|Ia_{i}(f)|');
 title('|Ia(f)|, |Ia_{r}(f)|, |Ia_{i}(f)| vs. f');
legend('|Ia|','|Ia_{r}|','|Ia_{i}|');
 
xg(1) = min(frq);
xg(2) = max(frq);
yg(1) = 1.0;
yg(2) = 1.0;
 
figure(40);
loglog(frq,rI,'b','linewidth',2);
hold    on;
loglog( xg,yg,'b','linewidth',2);
hold   off;
axis tight;
grid    on;
xlabel('Frequency, f (Hz)')
ylabel('rIa(f) = |Ia_{r}(f)/Ia_{i}(f)|');
 title('Ratio |Ia_{r}(f)/Ia_{i}(f)| vs. f');
 
%==========================================================================
% Plot Complex D Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReD,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re D (RMS mm)')
  title('Complex D');

subplot(4,2,2);
 plot(frq,PhD,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi D (degrees)')

subplot(4,2,3);
 plot(frq,ImD,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Im D (RMS mm)')
 
subplot(4,2,4);
 semilogx(frq,PhD,'b',10.0,0.0,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi D (degrees)')

subplot(4,2,5);
 plot(frq,MgD,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|D| (RMS mm)')

subplot(4,2,6);
 plot(frq,IFD,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi D)')

subplot(4,2,7);
%semilogy(frq,MgD,'b',xg,yg,'k');
 semilogx(frq,MgD,'b');
%loglog(frq,MgD,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|D| (RMS mm)')

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

%==========================================================================
% Plot Complex A Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReA,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re A (RMS mm/sec2)')
  title('Complex A');

subplot(4,2,2);
 plot(frq,PhA,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi A (degrees)')

subplot(4,2,3);
 plot(frq,ImA,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Im A (RMS mm/sec2)')
 
subplot(4,2,4);
 semilogx(frq,PhA,'b',10.0,0.0,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi A (degrees)')

subplot(4,2,5);
 plot(frq,MgA,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|A| (RMS mm/sec2)')

subplot(4,2,6);
 plot(frq,IFA,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi A)')

subplot(4,2,7);
%semilogy(frq,MgA,'b',xg,yg,'k');
 semilogx(frq,MgA,'b');
%loglog(frq,MgA,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|A| (RMS mm/sec2)')

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

%==========================================================================
% Loglog plots of |P(f)|, |U(f)|, |Z(f)|, |I(f)|, |D(f)| & |A(f)| vs. f
%==========================================================================
figure(51);
subplot(3,2,1);
loglog(frq,MgP,'b')
axis tight;
grid    on;
xlabel('Frequency, f (Hz)')
ylabel('|P(f)| (RMS Pa)');
 title('|P(f)|, |U(f)|, |Z(f)|, |I(f)|, |D(f)| and |A(f)| vs. f');

subplot(3,2,2);
loglog(frq,MgU,'b')
axis tight;
grid    on;
xlabel('Frequency, f (Hz)')
ylabel('|U(f)| (RMS mm/s)');

subplot(3,2,3);
loglog(frq,MgZ,'b')
axis tight;
grid    on;
xlabel('Frequency, f (Hz)')
ylabel('|Z(f)| (Ohms)');

subplot(3,2,4);
loglog(frq,MgI,'b')
axis tight;
grid    on;
xlabel('Frequency, f (Hz)')
ylabel('|I(f)| (RMS nW/m^{2})');

subplot(3,2,5);
loglog(frq,MgD,'b')
axis tight;
grid    on;
xlabel('Frequency, f (Hz)')
ylabel('|D(f)| (RMS mm)');

subplot(3,2,6);
loglog(frq,MgA,'b')
axis tight;
grid    on;
xlabel('Frequency, f (Hz)')
ylabel('|A(f)| (RMS mm/s^{2})');

%==========================================================================
%  Draw 3D complex plane for P
%==========================================================================
figure(61);
view(3);
plot3(ReP,ImP,frq,'b');
xlabel('Re(P) [RMS Pa]');
ylabel('Im(P) [RMS Pa]');
zlabel('Frequency (Hz)');
 title('P in the Complex Plane');
axis tight;
grid on;

%==========================================================================
%  Draw 3D complex plane for U
%==========================================================================
figure(62);
view(3);
plot3(ReU,ImU,frq,'b');
xlabel('Re(U) [RMS mm/s]');
ylabel('Im(U) [RMS mm/s]');
zlabel('Frequency (Hz)');
 title('U in the Complex Plane');
axis tight;
grid on;

%==========================================================================
%  Draw 3D complex plane for D
%==========================================================================
figure(63);
view(3);
plot3(ReD,ImD,frq,'b');
xlabel('Re(D) [RMS mm]');
ylabel('Im(D) [RMS mm]');
zlabel('Frequency (Hz)');
 title('D in the Complex Plane');
axis tight;
grid on;

%==========================================================================
%  Draw 3D complex plane for A
%==========================================================================
figure(64);
view(3);
plot3(ReA,ImA,frq,'b');
xlabel('Re(A) [RMS mm/s2]');
ylabel('Im(A) [RMS mm/s2]');
zlabel('Frequency (Hz)');
 title('A in the Complex Plane');
axis tight;
grid on;

%==========================================================================
%  Draw 3D complex plane for Z
%==========================================================================
figure(65);
view(3);
plot3(ReZ,ImZ,frq,'b');
xlabel('Re(Z) [Ohms]');
ylabel('Im(Z) [Ohms]');
zlabel('Frequency (Hz)');
 title('Z in the Complex Plane');
axis tight;
grid on;

%==========================================================================
%  Draw 3D complex plane for I
%==========================================================================
figure(66);
view(3);
plot3(ReI,ImI,frq,'b');
xlabel('Re(I) [RMS nW/m2]');
ylabel('Im(I) [RMS nW/m2]');
zlabel('Frequency (Hz)');
 title('I in the Complex Plane');
axis tight;
grid on;

%==========================================================================
% Plot Complex V Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReV,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re V (m/s)')
  title('Complex V');

subplot(4,2,2);
 plot(frq,PhV,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi V (degrees)')

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

subplot(4,2,5);
 plot(frq,MgV,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|V| (m/s)')

subplot(4,2,6);
 plot(frq,IFV,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi V)')

subplot(4,2,7);
%semilogy(frq,MgV,'b',xg,yg,'k');
 semilogx(frq,MgV,'b');
%loglog(frq,MgV,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|V| (m/s)')

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


%==========================================================================
% Plot Complex N Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReN,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re N')
  title('Complex N');

subplot(4,2,2);
 plot(frq,PhN,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi N (degrees)')

subplot(4,2,3);
 plot(frq,ImN,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Im N')
 
subplot(4,2,4);
 semilogx(frq,PhN,'b',10.0,0.0,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi N (degrees)')

subplot(4,2,5);
 plot(frq,MgN,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|N|')

subplot(4,2,6);
 plot(frq,IFN,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi NL)')

subplot(4,2,7);
%semilogy(frq,MgN,'b',xg,yg,'k');
 semilogx(frq,MgN,'b');
%loglog(frq,MgN,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|N|')

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


%==========================================================================
% Plot Complex L Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReL,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re L (m)')
  title('Complex L');

subplot(4,2,2);
 plot(frq,PhL,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi L (degrees)')

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

subplot(4,2,5);
 plot(frq,MgL,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|L| (m)')

subplot(4,2,6);
 plot(frq,IFL,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi L)')

subplot(4,2,7);
%semilogy(frq,MgL,'b',xg,yg,'k');
 semilogx(frq,MgL,'b');
%loglog(frq,MgL,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|L| (m)')

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


%==========================================================================
% Plot Complex K Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReK,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re K (m-1)')
  title('Complex K');

subplot(4,2,2);
 plot(frq,PhK,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi K (degrees)')

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

subplot(4,2,5);
 plot(frq,MgK,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|K| (m-1)')

subplot(4,2,6);
 plot(frq,IFK,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi K)')

subplot(4,2,7);
%semilogy(frq,MgK,'b',xg,yg,'k');
 semilogx(frq,MgK,'b');
%loglog(frq,MgK,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|K| (m-1)')

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


%==========================================================================
% Plot Complex G Data vs. Frequency - one graphics window, 4x2 format:
%==========================================================================
xg(1) =   0.0;
xg(2) = frq_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(frq,ReG,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Re G (RMS nkg/s-m2)')
  title('Complex G');

subplot(4,2,2);
 plot(frq,PhG,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Phi G (degrees)')

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

subplot(4,2,5);
 plot(frq,MgG,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|G| (RMS nkg/s-m2)')

subplot(4,2,6);
 plot(frq,IFG,'b',xg,yg,'k');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('Cos(Phi G)')

subplot(4,2,7);
%semilogy(frq,MgG,'b',xg,yg,'k');
 semilogx(frq,MgG,'b');
%loglog(frq,MgG,'b');
 axis tight
 grid on
 xlabel('Frequency (Hz)')
 ylabel('|G| (RMS nkg/s-m2)')

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

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

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

figure(111);

subplot(2,2,1);
 semilogy(frq,Wpot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Wpot(f) (RMS nJ/m3)');
  title('Wpot(f), Wkin(f), Wtot(f) vs. r');

subplot(2,2,2);
 semilogy(frq,Wkin,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Wkin(f) (RMS nJ/m3)');

subplot(2,2,3);
 semilogy(frq,Wtot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Wtot(f) (RMS nJ/m3)');

subplot(2,2,4);
 semilogy(frq,Rwpk,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Ratio Wpot(f)/Wkin(f)');

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

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

figure(112);

subplot(2,2,1);
 semilogy(frq,Fpot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Fpot');
  title('Fpot(f), Fkin(f), Ftot(f) vs. r');

subplot(2,2,2);
 semilogy(frq,Fkin,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Fkin(f)');

subplot(2,2,3);
 semilogy(frq,Ftot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Ftot(f)');

subplot(2,2,4);
 semilogy(frq,Rwpk,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Ratio Fpot(f)/Fkin(f)');

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

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

figure(113);

subplot(2,2,1);
 semilogy(frq,Wrad,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Wrad(f) (RMS nJ/m3)');
  title('Wrad(f), Wvrt(f), Wtot(f) vs. r');

subplot(2,2,2);
 semilogy(frq,Wvrt,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Wvrt(f) (RMS nJ/m3)');

subplot(2,2,3);
 semilogy(frq,Wtot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Wtot(f) (RMS nJ/m3)');

subplot(2,2,4);
 semilogy(frq,Rwrv,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Ratio Wrad(f)/Wvrt(f)');

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

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

figure(114);

subplot(2,2,1);
 semilogy(frq,Frad,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Frad(f)');
  title('Frad(f), Fvrt(f), Ftot(f) vs. r');

subplot(2,2,2);
 semilogy(frq,Fvrt,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Fvrt(f)');

subplot(2,2,3);
 semilogy(frq,Ftot,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Ftot(f)');

subplot(2,2,4);
 semilogy(frq,Rwrv,'b',xg,yg,'k');
 axis tight;
 grid    on;
 xlabel('Frequency, f (Hz)');
 ylabel('Ratio Frad(f)/Fvrt(f)');
 
 %==========================================================================
 beep();
 
 fprintf('\n Baffled_Circular_Piston_F_Theory Completed !!! \n')
%==========================================================================