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


disp('PART 1: Parallel mechanical oscillator impedance ==================')
disp('Use matlab "diff(z,t)" to do symbolic derivative d/dt, and use  ')
disp('    "diff(z,t,t)" to do second derivative when defining a symbolic')
disp('    second-order differential equation.  ')
disp('   Then, use matlab dsolve to solve, given initial conditions.')
disp(' ')


clear z(t) z0 m c n r Fin(t) s Zs z0 dzdt0 Vs Fins Zmech sFins sVs  solFins solZmech
syms  z(t) z0 m c n r Fin(t) s Zs z0 dzdt0 Vs Fins Zmech sFins sVs solFins solZmech
z0=0   % initial condition for position z
dzdt0=0 % initial condition for velocity v=dzdt
dampedDiffEq = m*diff(z,t,t)+ r*diff(z,t)+ z/n  == Fin(t) %damped differential equation
dampedZs1 = laplace(dampedDiffEq,t,s);            % Laplace to get Z(s)
dampedZs2 =subs(dampedZs1 ,laplace(z(t), t, s),Zs); % Define L{z(t)}=Zs
dampedZs3 =subs(dampedZs2 ,laplace(Fin(t), t, s),Fins); % Define L{Fin(t)}=Fins
dampedZs4 =subs(dampedZs3 ,z(0),0);                   % Define z(0)=0
dampedZs5 =subs(dampedZs4 ,subs(diff(z(t), t), t, 0),0);  % Define dzdt(0)=0
dampedZs6 =subs(dampedZs5 ,Zs,Vs/s)  ;                % replace Zs by Vs/s where L{v(t)}=Vs
eqs=[ dampedZs6 , Zmech ==  Fins/Vs] ;
vars=[Fins Zmech];
[ solFins solZmech] = solve(eqs,vars,'Real',true);
disp('======= solution for mechanical impedance: ')
Zmech=simplify(solZmech)  % mechanical impedance
disp('======= solution for mechanical admittance: ')
Ymech=1/Zmech


disp(' simulation parameters  ================================')
g=9.8                 %gravitational constant
m=1                  %mass
n=1                 %spring complance
k=1/n                 %spring stiffness
r=1                %mechanical resistance of damper
Q=sqrt(m/n)/r         %Q
omega0=sqrt(1/(m*n))  %natural frequency
freqHz=omega0/2/pi    %natural freq. in Hz
tau=1/(freqHz)        %period
z0=-m*g*n-1           %initial position
dzdt0=0               %initial velocity

disp('convert Zmech(s) into a matlab system transfer function  ===============')
%convert Z(s) into a matlab system transfer function
clear s; 
TFstr = char(subs(Zmech));
s = tf('s');
eval(['ZmechSys = ',TFstr]);
disp('======= transfer function form for mechanical impedance: ')
ZmechSys=minreal(ZmechSys)
disp('======= transfer function form for mechanical admittance: ')
YmechSys=minreal(1/ZmechSys)
disp('======= poles and zeroes mechanical impedance: ')
[Zmechpoles,Zmechzeroes] = pzmap(ZmechSys)


hh=figure(1);
fmax=1000
fmin=0.001
bplt=bodeplot(ZmechSys,'r')
setoptions(bplt,'PhaseVisible','on')
setoptions(bplt,'FreqUnits','rad/s','Grid','on','YLim',{[-20 60],[-100 100]},'XLim',[fmin  fmax]);
title({'Spring/Mass/Damper  Zmech(s)=F(s)/v(s) (force/velocity)',[' m=' num2str(m) ',  n=' num2str(n) ',  r=' num2str(r) ',  Q=' num2str(Q)  '  f0 Hz=' num2str(freqHz) ]})
set(findall(hh,'Type','line'),'LineWidth',3);
figax=findall(hh,'Type','Axes');
set(figax(1),'LineWidth',3);set(figax(2),'LineWidth',3);set(figax(3),'LineWidth',3);
set(findall(hh,'Type','text'),'FontSize',16);
set(findall(hh,'Type','text'),'FontWeight','bold');
disp('         Fig. 1.  Mechanical impedance Zmech(s). ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Change the simulation parameters to match the damped example')
disp('    in the lecture notes with m=2, n=1/3,  and r=1.22.')
disp('    (Fig. 1 simulation should change accordingly.)')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. After making the changes in Q1,')
disp('    Select the cursor tool on the plot, ')
disp('    what is the magnitude in dB and phase in degrees of Zmech(s) ')
disp('    at a frequency of 10 rad/s?')
disp('    (if you cannot measure at that frequency, choose the closest )')
disp('     measurable frequency within 10% of this frequency) ')
disp(' ')
disp('Answer: ??? dB  ??? degrees   at frequency=??? rad/s')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q1,')
disp('    Observe the transfer function ZmechSys above which equals Zmech(s), ')
disp('    and use this to COMPUTE  Zmech(s) at the same frequency in Q2. ')
disp('    What is the magnitude and phase in degrees of this computed Zmech(s)?')
disp('    HINT: consider matlab functions freqresp, abs, and angle')
disp(' ')
disp('Answer: |Zmech(s)|=???  phase= ??? degrees')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1,')
disp('    What is 10log10 the SQUARE of your answer for |Zmech(s)| in Q3?')
disp('    And does the answer approx. equal the answer in Q2')
disp('    Hint: the inductor impedance should be purely imaginary.')
disp(' ')
disp('Answer: 10log10(|Zmech(s)|)=??? dB     aprrox. same as Q2: yes/no')
disp(' ')
disp(' ')

 
disp('PART 2: RLC analog of parallel mechanical oscillator impedance ==================')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q1,')
disp('    according to lecture notes, the parallel mechanical oscillator')
disp('    is equivalent to a series RLC circuit.  From lecture notes,')
disp('    what are the values of R, L, and C for the example above?')
disp('    (GIVE THE NUMERICAL VALUES BELOW, NOT expressed as variables.)')
disp(' ')
disp('Answer: R=??? ohm,    L=??? H,    C=??? F')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. Using the values in Q5,')
disp('    what is the impedance of the inductor L at the frequency used in Q2,')
disp('    and does this approoximately equal your answer in Q3?')
disp(' ')
disp('Answer: impedance=???,    approx. equal: yes/no')
disp(' ')
disp(' ')

 
disp('PART 3: Acoustic analog of parallel mechanical oscillator impedance ==================')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. Given the mass density of air is rho=1.22 kg per cubic meter,')
disp('    and given compressibility of air is kappa= 0.00001 =10^-5 Pa^{-1},')
disp('    and given an air-filled pipe of 0.02 m radius,')
disp('    and for a 0.01 m length of air in this pipe')
disp('    according to lecture notes, ')
disp('    compute the mechanical equivalent mass m and equivalent spring compliance n.')
disp('    NOTE: be careful that n=1/k in comparison to A/kappa/l in lecture notes.')
disp('    (GIVE THE NUMERICAL VALUES BELOW)')
disp(' ')
disp('Answer: m= ??? kg,    n=???  N/(m/s^2)')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. Using the values in Q7, and acoustic analogies in lecture notes,')
disp('    what are the equivalent values of L and C for Q7?')
disp('    (GIVE THE NUMERICAL VALUES BELOW, NOT expressed as variables.)')
disp(' ')
disp('Answer: L= ??? H,     C= ??? F')
disp(' ')
disp(' ')

 
disp('PART 4: Spring+mass experiment movie ===================================')
disp('Please download the experimental movie from the course website')
disp('     and use the observations from the experiment to determine')
disp('     the constants for the damped-mass-spring experiment,')
disp('     (You may assume no damping for the example.)')
disp(' ')
disp(' ')


disp(' experimental parameters  ================================')
g=9.8                 %gravitational constant
mexp=1            %experimental mass
kexp=1       %experimental spring stiffness
nexp=1/kexp                %experimental spring compliance
rexp=0                %mechanical resistance of damper
omega0=sqrt(1/(mexp*nexp))  %natural frequency
freqHz=omega0/2/pi    %natural freq. in Hz
tau=1/(freqHz)        %period



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. Observe the spring+mass movie slinkMovie1.mp4,')
disp('    and based on your observations in the movie experiment, ')
disp('    MAKE NECESSARY CHANGES ABOVE to:')
disp('    experimental mass = mexp and experimental compliance = nexp,')
disp('    and set experimental mechanical resistance rexp=0.')
disp('    HINT: use static equilibrium in beginning of movie to compute spring stiffness')
disp('    HINT: be careful to use force NOT MASS in computing spring stiffness k=1/n')
disp('    (GIVE THE NUMERICAL VALUES BELOW, NOT expressed as variables.)')
disp(' ')
disp('Answer: mexp= ??? kg,     nexp= ??? m/N')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. Observe the spring+mass movie,')
disp('    and determine the observed oscillation frequency in the movie, ')
disp('    and compare this to the natural frequency in Hz computed from')
disp('    the experimental parameters of Q9.')
disp('    HINT: MAKE SURE that both answers are in Hz.')
disp('    Hint: your 2 answers should be within 20%.')
disp(' ')
disp('Answer: observed freq.= ??? Hz,  computed freq.= ??? Hz')
disp(' ')
disp(' ')