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


disp('PART 1: LC transmission line waves ==================')

disp(' ')

disp(' LC transmission line simulation parameters  ================')

%Electrical line parameters
global LR CR termination vp Z0 linelength pulsetau
LR=16  % indutance er unit length for right-handed line
CR=1  % capacitance er unit length for right-handed line
vp=1/sqrt(LR*CR)  %phase velocity
Z0=sqrt(LR/CR)   %impedance Z0
linelength=vp/2/10^6   % length of line for simulation
% select one of these if you want short circuit or open circuit boundaries
%termination='short' % short volts=0 or force=0
termination='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:(linelength/200):linelength;
disp(' solution points in time')
t=0:0.02/10^6:1.2/10^6
disp('now, solve differential equation')
sol = pdepe(m,@LCwaveEQpde,@LCwaveEQic,@LCwaveEQbc,x,t);
v = sol(:,:,1);   %voltage is first variable
i = sol(:,:,2);   %current is first variable

figure(1);
hh=surf(x,t,v);
title(['v(x,t) for LR=' num2str(LR) ' CR=' num2str(CR) ' vp=' num2str(vp) ' Z0=' num2str(Z0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 1. LC transmission line simulation. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Change the simulation parameters to match the RG-58 coaxial line')
disp('    example in the lecture notes with LR=259 nH/m, CR=94 pF/m.')
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= ??? ns ')
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 LR and CR.')
disp(' ')
disp('Answer: theoretical vg= ??? m/s      same as measured vg (+/-5%): yes/no')
disp(' ')
disp(' ')


 
disp('PART 2: mechanical transmission line waves ==================')
disp(' ')


disp(' mechanical transmission line simulation parameters  ================')


%Mechanical line parameters
global MR NR mechtermination mechvp mechZ0 mechlinelength 
MR=16  % mass per unit length for right-handed line
NR=1  % compliance per unit length for right-handed line
mechvp=1/sqrt(MR*NR)  %phase velocity
mechZ0=sqrt(MR/NR)   %impedance Z0
mechlinelength=mechvp/2  % length of line for simulation
% select one of these if you want short circuit or open circuit boundaries
%mechtermination='short' % short volts=0 or force=0
mechtermination='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:(mechlinelength/200):mechlinelength;
disp(' solution points in time')
t=0:0.02:1.2
disp('now, solve differential equation')
sol = pdepe(m,@MECHwaveEQpde,@MECHwaveEQic,@MECHwaveEQbc,x,t);
Force = sol(:,:,1);   %force is first variable
velocity = sol(:,:,2);   %velocity is first variable

figure(2);
hh=surf(x,t,Force);
title(['Force(x,t) for MR=' num2str(MR) ' NR=' num2str(NR) ' vp=' num2str(mechvp) ' Z0=' num2str(mechZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 2. Mechanical transmission line simulation. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. Change the simulation parameters to match the mass+spring line')
disp('    example in the lecture notes with MR=0.2 kg/m, and NR=5/3 .')
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 mechlinelength. ')
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.  Does this measured vg equal')
disp('     the theoretical vg? ')
disp(' ')
disp('Answer: vg= ??? m/s     equals theoretical (+/-5%): yes/no')
disp(' ')
disp(' ')


 
disp('PART 3: SLINKY transmission line waves ==================')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. Observe the experiment in the movie slinkWave1.mp4 ')
disp('    and USE THE SLOW MOTION VIDEO AT THE END of the movie ')
disp('    to measure the group velocity of the pulse vg.')
disp(' ')
disp('Answer: vg= ??? m/s     ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. Observe the experiment in the previous week movie slinkMovie1.mp4 ')
disp('    and determine the compliance n of the slinky. ')
disp(' ')
disp('Answer: compliance n= ??? m/N    ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. For the pulse-wave experiment in the movie slinkWave1.mp4 ')
disp('    with the slinky pulled to a length of 30 cm, ')
disp('    determine the mass-per-unit-length MR in kg/m  ')
disp('     and compliance=per-unit-length NR in N^(-1)')
disp('     and group velocity in m/s.')
disp('   Hint: first convert grams to kg, and cm to m.')
disp(' ')
disp('Answer: MR= ??? kg/m    NR= ??? N^(-1)    vg= ??? m/s  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q11. Edit the paramaters MR and NR and mechanicallength below here, ')
disp('    so that Fig. 3 below is a simulation of the slinky wave,. ')
disp('    using the measured parameters you derived in Q10. ')
disp(' ')
disp(' ')


disp(' mechanical transmission line simulation parameters  ================')


%Mechanical line parameters
global MR NR mechtermination mechvp mechZ0 mechlinelength 
MR=16  % mass per unit length for right-handed line
NR=1  % compliance per unit length for right-handed line
mechvp=1/sqrt(MR*NR)  %phase velocity
mechZ0=sqrt(MR/NR)   %impedance Z0
mechlinelength=0.8  % length of line for simulation
% select one of these if you want short circuit or open circuit boundaries
%mechtermination='short' % short volts=0 or force=0
mechtermination='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:(mechlinelength/200):mechlinelength;
disp(' solution points in time')
t=0:0.02:1
disp('now, solve differential equation')
sol = pdepe(m,@MECHwaveEQpde,@MECHwaveEQic,@MECHwaveEQbc,x,t);
Force = sol(:,:,1);   %force is first variable
velocity = sol(:,:,2);   %velocity is first variable

figure(3);
hh=surf(x,t,Force);
title(['Force(x,t) for MR=' num2str(MR) ' NR=' num2str(NR) ' vp=' num2str(mechvp) ' Z0=' num2str(mechZ0) ' ']);
xlabel('Distance x');
ylabel('Time t');
disp('         Fig. 3. SLINKY 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] = MECHwaveEQpde(x,t,u,DuDx)
global MR NR mechtermination mechvp mechZ0 mechlinelength 
c = [-NR ; -MR];
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 = MECHwaveEQic(x)
global MR NR mechtermination mechvp mechZ0 mechlinelength 
% 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*mechlinelength))*sin(pi*(x)/(0.2*mechlinelength)); ...
        (heaviside(x)-heaviside(x-0.2*mechlinelength))*sin(pi*(x)/(0.2*mechlinelength))/mechZ0];
end

% --------------------------------------------------------------------------
%Boundary Conditions
function [pl,ql,pr,qr] = MECHwaveEQbc(xl,ul,xr,ur,t)
global MR NR mechtermination mechvp mechZ0 mechlinelength  
if strcmp(mechtermination,'short')
  pl=[ul(1);  0  ]; ql = [0;  1 ]; pr = [ur(1); 0 ]; qr = [0; 1 ];
end
if strcmp(mechtermination,'open')
  pl=[0;  ul(2)  ]; ql = [1;  0 ]; pr = [0; ur(2) ]; qr = [1; 0 ];
end
end