% harmonic oscillator
%copyright 2020 by Thomas P. Weldon
% all rights reserved
clear all
close all hidden
clc

%some useful functions
syms u(t) pulse(t) triangle(t) 
sympref('HeavisideAtOrigin',1); %set equal to 1 as in unitstep
newHeaaviside0=heaviside(sym(0))
u(t)=heaviside(t) %unit step
pulse(t)=u(t+1/2)-u(t-1/2)
triangle(t)=(2*(t+1/2)-4*t*u(t))*pulse(t)


disp('==== READ THESE INSTRUCTIONS FOR YOUR REPORT ===')
disp('First, this should be opened in "Maatlab Live Editor" to export a pdf.' ) 
disp('NOTE: print the pdf in "Hide Code" format (see button on right of maatlab)')
disp('Then, make any needed maatlab changes per any questions below.') 
disp('You must re-run the script for changes to take effect.') 
disp('Dont forget to add the names of the group members as shown below.')
disp('Answer all numbered questions below, Q1, Q2, etc.') 
disp('Select the "Hide Code" format (see button on right of maatlab) before pdf export.')
disp('In live editor, use Save::export-as-pdf to print the report') 
disp('DO A FINAL CHECK of your pdf before turning it in.') 
disp(' ======================================================.')
disp(' ')
disp('ECGR 4122 Computer Project: Acoustic Waves')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')


disp('PART 1: WATER Acoustic transmission line waves ==================')

disp(' ')

disp(' Water acoustic transmission line simulation parameters  ================')

%Acoustic line parameters
global rho kappa acoustermination acousvp acousZ0 acouslinelength 
rho=16  % volume mass density kg/m^3
kappa=1  % compressibility Pa^(-1)
acousvp=1/sqrt(rho*kappa) %phase velocity
acousZ0=sqrt(rho/kappa)   %impedance Z0
acouslinelength=acousvp/2  % length of line for simulation
% select one of these if you want short circuit or open circuit boundaries
%acoustermination='short' % short volts=0 or force=0
acoustermination='open' % open current=0 or velocity=0


% --------------------------------------------------------------------------
%see https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html
m = 0 % parameter for diff eqn solver, see https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html
disp(' solution points in space')
x=0;x = 0:(acouslinelength/200):acouslinelength;
disp(' solution points in time')
t-0;t=0:0.02:1.2
disp('now, solve differential equation')
sol = pdepe(m,@ACOUSwaveEQpde,@ACOUSwaveEQic,@ACOUSwaveEQbc,x,t);
pressure = sol(:,:,1);   %force is first variable
velocity = sol(:,:,2);   %velocity is first variable

figure(1);
hh=surf(x,t,pressure);
title(['pressure(x,t) for rho=' num2str(rho) ' kappa=' num2str(kappa) ' vp=' num2str(acousvp) ' Z0=' num2str(acousZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 1. Water acoustic transmission line simulation. ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Change the simulation parameters to the proper parameters for the')
disp('    water example in the lecture notes with rho=1000 kg/m^3 and ')
disp('    kappa=4.6 x 10^(-10) Pa^(-1).')
disp('    (Fig. 1 simulation should change accordingly.)')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. After making the changes in Q1,')
disp('    observe Fig. 1, it shows round-trip travel of a pulse after')
disp('    reflecting from the far end of the transmission line.')
disp('    Select the cursor tool (data tips) on the plot, ')
disp('    what is the initial X-position at time=Y=0 and voltage=Z=1? ')
disp('    What is the  X-position at far end of the line (largest X-value)? ')
disp('    This X at the far end should equal the linelength. ')
disp('    What is the round-trip time Y at position=X=0 and voltage Z=2 ?')
disp('    (the Z-values yoou observe should be within a few percent )')
disp(' ')
disp('Answer: initial X= ??? m    far-end X= ??? m  round-trip time= ??? s ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q1,')
disp('    Use the data in Q2 to measure the group velocity vg, ')
disp('    noting that the pulse did not start at X=0, so distance traveled is')
disp('    less than 2 time the linelength.')
disp(' ')
disp('Answer: vg= m/s')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1,')
disp('    Compare your Q3 measured group velocity vg, ')
disp('    to the theoretical group velocity computed from rho and kappa.')
disp(' ')
disp('Answer: theoretical vg= ??? m/s   same as measured vg (+/-5%): yes/no')
disp(' ')
disp(' ')



disp('PART 2: AIR Acoustic transmission line waves ==================')

disp(' ')

disp(' Air acoustic transmission line simulation parameters  ================')

%Acoustic line parameters
global rho kappa acoustermination acousvp acousZ0 acouslinelength 
rho=16  % volume mass density kg/m^3
kappa=1  % compressibility Pa^(-1)
acousvp=1/sqrt(rho*kappa) %phase velocity
acousZ0=sqrt(rho/kappa)   %impedance Z0
acouslinelength=acousvp/2  % length of line for simulation
% select one of these if you want short circuit or open circuit boundaries
%acoustermination='short' % short volts=0 or force=0
acoustermination='open' % open current=0 or velocity=0


% --------------------------------------------------------------------------
%see https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html
m = 0 % parameter for diff eqn solver, see https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html
disp(' solution points in space')
x=0;x = 0:(acouslinelength/200):acouslinelength;
disp(' solution points in time')
t=0;t=0:0.02:1.2
disp('now, solve differential equation')
sol = pdepe(m,@ACOUSwaveEQpde,@ACOUSwaveEQic,@ACOUSwaveEQbc,x,t);
pressure = sol(:,:,1);   %force is first variable
velocity = sol(:,:,2);   %velocity is first variable

figure(2);
hh=surf(x,t,pressure);
title(['pressure(x,t) for rho=' num2str(rho) ' kappa=' num2str(kappa) ' vp=' num2str(acousvp) ' Z0=' num2str(acousZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 2. Air acoustic transmission line simulation. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. Change the simulation parameters to the proper parameters for the')
disp('    air example in the lecture notes with rho=1.2 kg/m^3 and ')
disp('    kappa=7.14 x 10^(-6) Pa^(-1).')
disp('    (Fig. 2 simulation should change accordingly.)')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. After making the changes in Q5,')
disp('    observe Fig. 2, it shows round-trip travel of a pulse after')
disp('    reflecting from the far end of the transmission line.')
disp('    Select the cursor tool (data tips) on the plot, ')
disp('    what is the initial X-position at time=Y=0 and voltage=Z=1? ')
disp('    What is the  X-position at far end of the line (largest X-value)? ')
disp('    This X at the far end should equal the linelength. ')
disp('    What is the round-trip time Y at position=X=0 and voltage Z=2 ?')
disp('    (the Z-values yoou observe should be within a few percent )')
disp(' ')
disp('Answer: initial X= ??? m    far-end X= ??? m  round-trip time= ??? s ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. After making the changes in Q5,')
disp('    Use the data in Q6 to measure the group velocity vg, ')
disp('    noting that the pulse did not start at X=0, so distance traveled is')
disp('    less than 2 time the linelength.')
disp(' ')
disp('Answer: vg= ??? m/s')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. After making the changes in Q5,')
disp('    Compare your Q7 measured group velocity vg, ')
disp('    to the theoretical group velocity computed from rho and kappa.')
disp(' ')
disp('Answer: theoretical vg= ??? m/s   same as measured vg (+/-5%): yes/no')
disp(' ')
disp(' ')






disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. For acoustic pulse experiment in the movie pvcPipe35inchPulseVideo.mp4 ')
disp('    measure the round-trip time in the pipe (ONE cycle of observed audio ), ')
disp('    and determine the group velocity vg of the pulse in m/s.  ')
disp('   Hint: first convert inches to meters, round-trip is 2 times pipe length.')
disp(' ')
disp('Answer: Period of one cycle= ??? ms         vg= ??? m/s  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. Edit the paramaters rho and kappa and acouslinelength below here, ')
disp('    so that Fig. 3 below is a simulation of the pvc pipe wave,. ')
disp('    using the air parameters you used in Q5 and length of the pvc pipe. ')
disp(' ')
disp(' ')


%Acoustic line parameters
global rho kappa acoustermination acousvp acousZ0 acouslinelength 
rho=16  % volume mass density kg/m^3
kappa=1  % compressibility Pa^(-1)
acousvp=1/sqrt(rho*kappa) %phase velocity
acousZ0=sqrt(rho/kappa)   %impedance Z0
acouslinelength=1/4  % length of line for simulation
% select one of these if you want short circuit or open circuit boundaries
%acoustermination='short' % short volts=0 or force=0
acoustermination='open' % open current=0 or velocity=0



% --------------------------------------------------------------------------
%see https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html
m = 0 % parameter for diff eqn solver, see https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html
disp(' solution points in space')
x=0;x = 0:(acouslinelength/200):acouslinelength;
disp(' solution points in time')
t=0;t=0:0.1/1000:10/1000
disp('now, solve differential equation')
sol = pdepe(m,@ACOUSwaveEQpde,@ACOUSwaveEQic,@ACOUSwaveEQbc,x,t);
pressure = sol(:,:,1);   %force is first variable
velocity = sol(:,:,2);   %velocity is first variable

figure(3);
hh=surf(x,t,pressure);
title(['pressure(x,t) for rho=' num2str(rho) ' kappa=' num2str(kappa) ' vp=' num2str(acousvp) ' Z0=' num2str(acousZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 3. PVC pipe with air acoustic transmission line simulation. ')
disp(' ')
disp(' ')

 






  
 % THESE FUNCTIONS DEFINE THE DIFFERENTIAL EQUATIONS ETC
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % this is the pair of telegrapher's equations
%see https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html
%

%

function [c,f,s] = LCwaveEQpde(x,t,u,DuDx)
global LR CR termination vp Z0 linelength
c = [-CR ; -LR];
f = [u(2); u(1)] ;
% y = u(1) - u(2);
% F = exp(5.73*y)-exp(-11.47*y);
s = [0; 0];
end

% 
% v(x,0)=0, i(x,0)=0 initial conditions
function u0 = LCwaveEQic(x)
global LR CR termination vp Z0 linelength
% u0 = [(heaviside(x-0.4*linelength)-heaviside(x-0.6*linelength))*sin(pi*(x-0.4)/(0.2*linelength)); ...
%         (heaviside(x-0.4*linelength)-heaviside(x-0.6*linelength))*sin(pi*(x-0.4)/(0.2*linelength))/Z0];
% u0 = [(heaviside(x-0.4*linelength)-heaviside(x-0.6*linelength))*sin(pi*(x-0.4)/(0.2*linelength)); ...
%         (heaviside(x-0.4*linelength)-heaviside(x-0.6*linelength))*sin(pi*(x-0.4)/(0.2*linelength))/Z0];
u0 = [(heaviside(x)-heaviside(x-0.2*linelength)).*sin(pi*(x)/(0.2*linelength)); ...
        (heaviside(x)-heaviside(x-0.2*linelength)).*sin(pi*(x)/(0.2*linelength))/Z0];
end

% --------------------------------------------------------------------------
%Boundary Conditions
function [pl,ql,pr,qr] = LCwaveEQbc(xl,ul,xr,ur,t)
global LR CR termination vp Z0 linelength pulsetau
if strcmp(termination,'short')
  pl=[ul(1);  0  ]; ql = [0;  1 ]; pr = [ur(1); 0 ]; qr = [0; 1 ];
end
if strcmp(termination,'open')
  pl=[0;  ul(2)  ]; ql = [1;  0 ]; pr = [0; ur(2) ]; qr = [1; 0 ];
end
end
 
  
 % THESE FUNCTIONS DEFINE THE DIFFERENTIAL EQUATIONS ETC
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % this is the pair of telegrapher's equations
%see https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html
%
function [c,f,s] = ACOUSwaveEQpde(x,t,u,DuDx)
global rho kappa acoustermination acousvp acousZ0 acouslinelength 
c = [-kappa ; -rho];
f = [u(2); u(1)] ;
% y = u(1) - u(2);
% F = exp(5.73*y)-exp(-11.47*y);
s = [0; 0];
end

% 
% v(x,0)=0, i(x,0)=0 initial conditions
function u0 = ACOUSwaveEQic(x)
global rho kappa acoustermination acousvp acousZ0 acouslinelength 
% u0 = [(heaviside(x-0.4*linelength)-heaviside(x-0.6*linelength))*sin(pi*(x-0.4)/(0.2*linelength)); ...
%         (heaviside(x-0.4*linelength)-heaviside(x-0.6*linelength))*sin(pi*(x-0.4)/(0.2*linelength))/Z0];
% u0 = [(heaviside(x-0.4*linelength)-heaviside(x-0.6*linelength))*sin(pi*(x-0.4)/(0.2*linelength)); ...
%         (heaviside(x-0.4*linelength)-heaviside(x-0.6*linelength))*sin(pi*(x-0.4)/(0.2*linelength))/Z0];
u0 = [(heaviside(x)-heaviside(x-0.2*acouslinelength)).*sin(pi*(x)/(0.2*acouslinelength)); ...
        (heaviside(x)-heaviside(x-0.2*acouslinelength)).*sin(pi*(x)/(0.2*acouslinelength))/acousZ0];
end

% --------------------------------------------------------------------------
%Boundary Conditions
function [pl,ql,pr,qr] = ACOUSwaveEQbc(xl,ul,xr,ur,t)
global rho kappa acoustermination acousvp acousZ0 acouslinelength  
if strcmp(acoustermination,'short')
  pl=[ul(1);  0  ]; ql = [0;  1 ]; pr = [ur(1); 0 ]; qr = [0; 1 ];
end
if strcmp(acoustermination,'open')
  pl=[0;  ul(2)  ]; ql = [1;  0 ]; pr = [0; ur(2) ]; qr = [1; 0 ];
end
end