%  
%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: Strings and Transverse Waves')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')


disp(' ')

disp(' Initialize guitar string simulation parameters  ================')

% string parameters
global rho tension stringtermination stringvp stringZ0 stringlength mode 
rho=0.02 % linear mass density of the string kg/m
tension=80  % compressibility Pa^(-1)
stringvp=1/sqrt(rho/tension) %phase velocity
stringZ0=sqrt(rho*tension)   %impedance Z0
stringlength=stringvp/2  % length of line for simulation
stringlength=0.66
stringFundamentalHz=stringvp/stringlength/2
% select one of these if you want short circuit or open circuit boundaries
%stringtermination='short' % short volts=0 or force=0
stringtermination='open' % open current=0 or velocity=0

 
disp('=========================================================== ')
disp(' ')
disp('PART 1: Pulse traveling wave on a string ==================')
disp(' ')
disp('=========================================================== ')


% --------------------------------------------------------------------------
%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:(stringlength/250):stringlength;
disp(' solution points in time')
t-0;t=0:4*stringlength/stringvp/50:4*stringlength/stringvp;
disp('now, solve differential equation')
%use pulse traveling wave initial conditions
sol = pdepe(m,@STRINGwaveEQpde,@STRINGpulsewaveEQic,@STRINGwaveEQbc,x,t);
%use standing wave initial conditions
%sol = pdepe(m,@STRINGwaveEQpde,@STRINGstandingwaveEQic,@STRINGwaveEQbc,x,t);
force = sol(:,:,1);   %force is first variable
velocity = sol(:,:,2);   %velocity is first variable

figure(1);
hh=surf(x,t,velocity);
title(['velocity(x,t) for rho=' num2str(rho) ' tension=' num2str(tension) ' vp=' num2str(stringvp) ' Z0=' num2str(stringZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 1. Velocity for string traveling wave pulse simulation. ')
disp(' ')
disp(' ')


figure(2);
hh=surf(x,t,force);
title(['force(x,t) for rho=' num2str(rho) ' tension=' num2str(tension) ' vp=' num2str(stringvp) ' Z0=' num2str(stringZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 2. Force for string traveling wave pulse simulation. ')
disp(' ')
disp(' ')


hh=figure(3);
sp=subplot(2,1,1);
plot(x,velocity(1,:),'b','LineWidth', 2);  grid('on');ylim([-2, 2]);xlim([0 max(x)]);
title({'String Pulse Traveling Wave', ['rho=' num2str(rho) ' tension=' num2str(tension) ' vp=' num2str(stringvp) ' Z0=' num2str(stringZ0)  ],'Velocity'});
xlabel('Distance x (m)'); ylabel('Velocity v_y (m/s)');
sp=subplot(2,1,2);plot(x,force(1,:),'r','LineWidth', 2);  grid('on');ylim([-2, 2]);xlim([0 max(x)]);
title('Force');xlabel('Distance x (m)'); ylabel('Force F_y (N)');
disp('         Fig. 3. Initial condition for string traveling wave pulse simulation. ')
disp(' ')
disp(' ')


%create a video of the string standing wave
vid = VideoWriter('stringPulseTravelingWaveMovie.mp4','MPEG-4');
vid.Quality = 95;vid.FrameRate=10 ;
open(vid);
hh=figure(100);
for nt=1:length(t)
    sp=subplot(2,1,1);
    plot(x,velocity(nt,:),'b','LineWidth', 2);  grid('on');ylim([-2, 2]);xlim([0 max(x)]);
    title({'String Pulse Traveling Wave', ['rho=' num2str(rho) ' tension=' num2str(tension) ' vp=' num2str(stringvp) ' Z0=' num2str(stringZ0)  ],'Velocity'});
    xlabel('Distance x (m)'); ylabel('Velocity v_y (m/s)');
    sp=subplot(2,1,2);plot(x,force(nt,:),'r','LineWidth', 2);  grid('on');ylim([-2, 2]);xlim([0 max(x)]);
    title('Force');xlabel('Distance x (m)'); ylabel('Force F_y (N)');
   frame = getframe(gcf);
    writeVideo(vid,frame);
end
close(vid);
% win=implay('stringPulseTravelingWaveMovie.mp4')
% disp('         Fig. IGNORE: This figure is only used to create videos. ')
% disp(' ')
% disp(' ')
close(figure(100))


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Change the simulation parameters to the proper parameters for the')
disp('    low-E guitar string example in the lecture notes with rho=0.007 kg/m ')
disp('    and tension Ts=84 N.')
disp('    (Figs. 1-3 simulation should change accordingly.)')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. After making the changes in Q1,')
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 at force Fy approx. 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 you observe should be within a few percent )')
disp('    Also record the initial force Fy at t=0 at the initial pulse location.')
disp(' ')
disp('Answer: initial X=??? m    far-end X=??? m  round-trip time=???? s ')
disp('        initial Fy=???')
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 times 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 tension.')
disp(' ')
disp('Answer: theoretical vg=??? m/s   same as measured vg (+/-5%): yes')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. 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, and meaasure the')
disp('    velocity vy at time=0 and at the same X-position as in Q2')
disp('     (velocity is z-value in Fig. 1 at measured at time=0) ')
disp('    What is the measured chanracteristic impedance Zo=Fy/vy using ')
disp('    the measured values at time=0 for Fy in Q2 and vy in Q5 ?')
disp('    What is the theorretical Zo?')
disp(' ')
disp('Answer: at t=0, vy=??? m/s   Zo=Fy/vy= ??? N/(m/s)  theory Zo=???? N/(m/s)')
disp(' ')
disp(' ')



disp('=========================================================== ')
disp(' ')
disp('PART 2: Standing wave on a string ==================')
disp(' ')
disp('=========================================================== ')

disp(' set desired mode of standing wave')
mode=3

% --------------------------------------------------------------------------
%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:(stringlength/250):stringlength;
disp(' solution points in time')
t-0;t=0:4*stringlength/stringvp/50:4*stringlength/stringvp;
disp('now, solve differential equation')
%use pulse traveling wave initial conditions
%sol = pdepe(m,@STRINGwaveEQpde,@STRINGpulsewaveEQic,@STRINGwaveEQbc,x,t);
%use standing wave initial conditions
sol = pdepe(m,@STRINGwaveEQpde,@STRINGstandingwaveEQic,@STRINGwaveEQbc,x,t);
force = sol(:,:,1);   %force is first variable
velocity = sol(:,:,2);   %velocity is first variable

figure(4);
hh=surf(x,t,velocity);
title(['velocity(x,t) for rho=' num2str(rho) ' tension=' num2str(tension) ' vp=' num2str(stringvp) ' Z0=' num2str(stringZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 4. Velocity for string standing wave simulation. ')
disp(' ')
disp(' ')


figure(5);
hh=surf(x,t,force);
title(['force(x,t) for rho=' num2str(rho) ' tension=' num2str(tension) ' vp=' num2str(stringvp) ' Z0=' num2str(stringZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 5. Force for string standing wave simulation. ')
disp(' ')
disp(' ')

hh=figure(6);
sp=subplot(2,1,1);
plot(x,velocity(1,:),'b','LineWidth', 2);  grid('on');ylim([-2, 2]);xlim([0 max(x)]);
title({'String Standing Wave', ['rho=' num2str(rho) ' tension=' num2str(tension) ' vp=' num2str(stringvp) ' Z0=' num2str(stringZ0)  ],'Velocity'});
xlabel('Distance x (m)'); ylabel('Velocity v_y (m/s)');
sp=subplot(2,1,2);plot(x,force(1,:),'r','LineWidth', 2);  grid('on');ylim([-2, 2]);xlim([0 max(x)]);
title('Force');xlabel('Distance x (m)'); ylabel('Force F_y (N)');
disp('         Fig. 6. Initial condition for string standing wave simulation. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. Change the standing wave simulation parameters mode n=1')
disp('    to observe the fundamental mode. ')
disp('    (Fig. 4 simulation should change accordingly.)')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. After making the changes in Q6,')
disp('    observe Fig. 4, it shows the velocity vy standing wave,')
disp('    with boundary conditions vy=0 at x=0 and x=0.66.')
disp('    Select the cursor tool (data tips) on the plot, ')
disp('    what is the initial velocity vy at time=Y=0 and at the midpoint X=0.33? ')
disp('    As the velocity at the midpoint X=0.33 decreases, then peaks again. ')
disp('    At what time does the second peak occur? ')
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 vy=??? m/s    second peak at t=???? s ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. After making the changes in Q6,')
disp('    Use the data in Q7 calculate the frequency of the standing wave. ')
disp('    Also, calculate the theoretical fundamental frequency of the fundamental')
disp('    mode of the standing wave.')
disp(' ')
disp('Answer: measured frequency=??? Hz     theory=??? Hz')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. After making the changes in Q6,')
disp('    In Fig. 4, measure the velocity vy at X=0 and for various times. ')
disp('    At x=0, is the velocity always zero or very close to zero?.')
disp(' ')
disp('Answer: velocity fr X=0 is nearly zero for all times: yes/no')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. After making the changes in Q6,')
disp('    In Fig. 5, measure the Forcf vy at the midpoint X=0.33 for various times. ')
disp('    At x=0/33, is the force always zero or very close to zero?.')
disp(' ')
disp('Answer: force at X=0.33 is nearly zero for all times: yes/no')
disp(' ')
disp(' ')



%create a video of the string standing wave
vid = VideoWriter('stringStandingWaveMovie.mp4','MPEG-4');
vid.Quality = 95;vid.FrameRate=10 ;
open(vid);
hh=figure(100);
for nt=1:length(t)
    sp=subplot(2,1,1);
    plot(x,velocity(nt,:),'b','LineWidth', 2);  grid('on');ylim([-2, 2]);xlim([0 max(x)]);
    title({'String Standing Wave', ['rho=' num2str(rho) ' tension=' num2str(tension) ' vp=' num2str(stringvp) ' Z0=' num2str(stringZ0)  ],'Velocity'});
    xlabel('Distance x (m)'); ylabel('Velocity v_y (m/s)');
    sp=subplot(2,1,2);plot(x,force(nt,:),'r','LineWidth', 2);  grid('on');ylim([-2, 2]);xlim([0 max(x)]);
    title('Force');xlabel('Distance x (m)'); ylabel('Force F_y (N)');
   frame = getframe(gcf);
    writeVideo(vid,frame);
end
close(vid);
% disp('         Fig.  IGNORE: This figure is only used to create videos. ')
% disp(' ')
% disp(' ')
close(figure(100))


win=implay('stringPulseTravelingWaveMovie.mp4')
win=implay('stringStandingWaveMovie.mp4')

 

  
 % 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

 % 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] = STRINGwaveEQpde(x,t,u,DuDx)
global rho tension stringtermination stringvp stringZ0 stringlength 
c = [-1/tension ; -rho];
f = [u(2); u(1)] ;
% y = u(1) - u(2);
% F = exp(5.73*y)-exp(-11.47*y);
s = [0; 0];
end

% 
% pulse traveling wave initial conditions
function u0 = STRINGpulsewaveEQic(x)
global rho tension stringtermination stringvp stringZ0 stringlength 
%traveling wave initial conditions
u0 = [(heaviside(x)-heaviside(x-0.2*stringlength)).*sin(pi*(x)/(0.2*stringlength)); ...
        (heaviside(x)-heaviside(x-0.2*stringlength)).*sin(pi*(x)/(0.2*stringlength))/stringZ0];
%standing wave initial conditions, 
%divide one by 1000 since when voltages reinforce, currents go to zero in a standing wave
% u0 = [(heaviside(x)-heaviside(x-0.99*stringlength)).*sin(pi*(x)/(0.99*stringlength))/1000; ...
%         (heaviside(x)-heaviside(x-0.99*stringlength)).*sin(pi*(x)/(0.99*stringlength))/stringZ0];
end


% 
% standing wave initial conditions
function u0 = STRINGstandingwaveEQic(x)
global rho tension stringtermination stringvp stringZ0 stringlength mode
%traveling wave initial conditions
% u0 = [(heaviside(x)-heaviside(x-0.2*stringlength)).*sin(pi*(x)/(0.2*stringlength)); ...
%         (heaviside(x)-heaviside(x-0.2*stringlength)).*sin(pi*(x)/(0.2*stringlength))/stringZ0];
%standing wave initial conditions, 
%divide one by 1000 since when voltages reinforce, currents go to zero in a standing wave
u0 = [(heaviside(x)-heaviside(x-0.99*stringlength)).*sin(mode*pi*(x)/(0.99*stringlength))/1000; ...
        (heaviside(x)-heaviside(x-0.99*stringlength)).*sin(mode*pi*(x)/(0.99*stringlength))/stringZ0];
end

% --------------------------------------------------------------------------
%Boundary Conditions
function [pl,ql,pr,qr] = STRINGwaveEQbc(xl,ul,xr,ur,t)
global rho tension stringtermination stringvp stringZ0 stringlength  
if strcmp(stringtermination,'short')
  pl=[ul(1);  0  ]; ql = [0;  1 ]; pr = [ur(1); 0 ]; qr = [0; 1 ];
end
if strcmp(stringtermination,'open')
  pl=[0;  ul(2)  ]; ql = [1;  0 ]; pr = [0; ur(2) ]; qr = [1; 0 ];
end
end