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

disp('PART 1: Basic Useful Functions ================================')
%create symbolic function
syms omega t x(t) y(t) g(t) h(t)
f=1000
omega= f*2*pi
g(t)=pulse(400*t)
x(t)=cos(omega*t)*exp(-500*t)*u(t)
%x(t)=pulse(f*(t-0.005))
y(t)=triangle(200*t)
h(t)=500*t*pulse(500*t)
hh=figure(1);
subplot(2,2,1);fplot(x,[-0.01 0.01],'LineWidth',3);grid on;
title('x(t)');xlabel('Time (s)'); ylabel('x(t)');
subplot(2,2,2);fplot(y,[-0.01 0.01],'LineWidth',3);grid on;
title('y(t)');xlabel('Time (s)'); ylabel('y(t)');
set(findall(hh,'Type','text'),'FontSize',16);
subplot(2,2,3);fplot(g,[-0.01 0.01],'LineWidth',3);grid on;
title('g(t)');xlabel('Time (s)'); ylabel('g(t)');
subplot(2,2,4);fplot(h,[-0.01 0.01],'LineWidth',3);grid on;
title('h(t)');xlabel('Time (s)'); ylabel('h(t)');
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 1.  Time functions. ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Change the function y(t) so the width of the triangle is 0.01 s')
disp('    in Fig. 1, and the height remains 1.')
disp('    Your turned-in plots should change accordingly.')
disp(' ')
disp(' ')



 



disp('PART 2: Undamped mechanical oscillator ============================')
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 F(t) g s Zs aZs dzdt(t) dzdt0
syms  z(t) z0 m c n r F(t) g s Zs aZs dzdt(t) dzdt0 
undampedDiffEq = m*diff(z,t,t)   + z/n  == -m*g %undamped differential equation
dzdt = diff(z,t);                               %dz/dt
cond = [z(0) == z0, dzdt(0) == dzdt0];          %initial conditions for z(0) and dz/dt(0)
undampedzt = dsolve(undampedDiffEq, cond)       %solve for z(t)
undampedZs = laplace(undampedzt,t,s)            % Laplace to get Z(s)
undampedvt = diff(undampedzt,t)                 %velocity v(t)=dz/dt
 
disp(' simulation parameters  ================================')
g=9.8                 %gravitational constant
m=1                   %mass
n=1                 %spring complance
k=1/n                 %spring stiffness
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
%video settings
maxtime=round(5*tau);
if tau<0.001 maxtime=round(5*10^6*tau)/10^6; end
if tau<0.1 maxtime=round(50*100*tau)/100; end  
maxtime   %display maxtime for plots in the printout
framerate=30;

% print out the results
defdigs=digits;    %save initial precision of printout
digits(4);         %lower precision
disp('==================================================')
disp('Solution for position z(t): ')
undampedzt=simplify(vpa(subs(undampedzt)),'Criterion','preferReal','Steps',100)
disp('Solution for velocity v(t): ')
undampedvt=simplify(vpa(subs(undampedvt)),'Criterion','preferReal','Steps',100)
disp('Solution for Z(s): ')
undampedZs=simplify(vpa(subs(undampedZs)),'Criterion','preferReal','Steps',100)
[numundampedZs,denundampedZs] = numden(vpa(subs(expand(undampedZs))));
digits(defdigs);   %restore initial precision of printout

for nt=1:200
    tt(nt)=(nt-1)*maxtime/200;
    t=tt(nt);
    zz(nt)=real(double(vpa(subs(undampedzt))));
    if nt>2 
        vv(nt)=(zz(nt)-zz(nt-1))/(tt(nt)-tt(nt-1));
    end
end
 vv(1)=vv(2);

hh=figure(2);
plot(tt,zz,'r',tt,vv,'b','LineWidth', 2);  grid('on');
title({'Undamped Spring+Mass z(t) in red and dz/dt in blue', [' m=' num2str(m)  ',  n=' num2str(n) ',  f0 Hz=' num2str(freqHz) ]})
xlabel('Time (s)'); ylabel('z (m)and dz/dt (m/s)')
gca.XGrid='on';gca.YGrid='on';gca.LineWidth=3;gca.FontWeight='bold';
set(findall(hh,'Type','line'),'LineWidth',3);
set(findall(hh,'Type','text'),'FontSize',14);
set(findall(hh,'Type','text'),'FontWeight','bold');
xlim([0, maxtime]);%ylim([-8, 2]);
disp('         Fig. 2.  Position z(t) and velocity v(t)=dz/dt. ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. Change the simulation parameters to match the undamped example')
disp('    in the lecture notes with m=2, k=3, z(0)=-7.53, and v(0)=0.')
disp('    (Fig. 2 simulation should change accordingly.)')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q2,')
disp('    does the time domain solution for position z(t) saved as undampedzt')
disp('    equal the time-domain soluion in the lecture notes example?')
disp('    (ignore any small errors in coefficients less than 5%)')
disp(' ')
disp('Answer:  yes/no ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q2,')
disp('    does the time domain solution for velocity v(t) saved as undampedvt')
disp('    equal the time-domain soluion in the lecture notes example?')
disp('    (ignore any small errors in coefficients less than 5%)')
disp(' ')
disp('Answer:  yes/no ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q2,')
disp('    does simulation in Fig. 2 show the correnct initial conditions')
disp('    at t=0 for the initial position z(0) and initial velocity v(0)?')
disp(' ')
disp('Answer:  yes/no ')
disp(' ')
disp(' ')



disp('PART 3: Damped mechanical oscillator ============================')
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 F(t) g s Zs aZs dzdt(t) dzdt0
syms  z(t) z0 m c n r F(t) g s Zs aZs dzdt(t) dzdt0 
dampedDiffEq = m*diff(z,t,t)+ r*diff(z,t)+ z/n  == -m*g %damped differential equation
dzdt = diff(z,t);                           %dz/dt
cond = [z(0) == z0, dzdt(0) == dzdt0];       %initial conditions for z(0) and dz/dt(0)
dampedzt = dsolve(dampedDiffEq, cond)       %solve for z(t)
dampedZs = laplace(dampedzt,t,s)            % Laplace to get Z(s)
dampedvt = diff(dampedzt,t)                 %velocity v(t)=dz/dt
 
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

%video settings
maxtime=round(5*tau);
if tau<0.001 maxtime=round(5*10^6*tau)/10^6; end
if tau<0.1 maxtime=round(50*100*tau)/100; end  
maxtime   %display maxtime for plots in the printout
framerate=30;

% print out the results
defdigs=digits;    %save initial precision of printout
digits(4);         %lower precision
disp('==================================================')
disp('Solution for position z(t): ')
dampedzt=simplify(vpa(subs(dampedzt)),'Criterion','preferReal','Steps',100)
disp('Solution for velocity v(t): ')
dampedvt=simplify(vpa(subs(dampedvt)),'Criterion','preferReal','Steps',100)
disp('Solution for Z(s): ')
dampedZs=simplify(vpa(subs(dampedZs)),'Criterion','preferReal','Steps',100)
[numdampedZs,dendampedZs] = numden(vpa(subs(expand(dampedZs))));
digits(defdigs);   %restore initial precision of printout

for nt=1:200
    tt(nt)=(nt-1)*maxtime/200;
    t=tt(nt);
    zz(nt)=real(double(vpa(subs(dampedzt))));
    if nt>2 
        vv(nt)=(zz(nt)-zz(nt-1))/(tt(nt)-tt(nt-1));
    end
end
 vv(1)=vv(2);

hh=figure(3);
plot(tt,zz,'r',tt,vv,'b','LineWidth', 2);  grid('on');
title({'Damped Spring+Mass z(t) in red and dz/dt in blue', [' m=' num2str(m)  ', n=' num2str(n)  ', r=' num2str(r) ', f0 Hz=' num2str(freqHz) ]})
xlabel('Time (s)'); ylabel('z (m)and dz/dt (m/s)')
gca.XGrid='on';gca.YGrid='on';gca.LineWidth=3;gca.FontWeight='bold';
set(findall(hh,'Type','line'),'LineWidth',3);
set(findall(hh,'Type','text'),'FontSize',14);
set(findall(hh,'Type','text'),'FontWeight','bold');
xlim([0, maxtime]);%ylim([-8, 2]);
disp('         Fig. 3.  Position z(t) and velocity v(t)=dz/dt. ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. Change the simulation parameters to match the damped example')
disp('    in the lecture notes with m=2, k=3, r=1.22, z(0)=-7.53, and v(0)=0.')
disp('    (Fig. 3 simulation should change accordingly.)')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. After making the changes in Q6,')
disp('    does the time domain solution for position z(t) ')
disp('    return to the static equilibrium position at t=20 s?')
disp(' ')
disp('Answer:  yes/no ')
disp(' ')
disp(' ')


disp('PART 4: Laplace analysis of damped mechanical oscillator ================')
disp('Use for the damped system with sinusoudal excitation, now consider  ')
disp('    Zm(s)/Fin(s) from sinusoudal Laplace analysis from lecture notes,')
disp(' ')


clear s  ZdampedZmOVERFins
syms  s  ZdampedZmOVERFins
dampedZmOVERFins = 1/(m*s^2+ r*s+ 1/n)     %Zm(s)/Fin(s) sinusoudal Laplace analysis from lecture notes

%convert Z(s) into a matlab system transfer function
clear s; 
TFstr = char(subs(dampedZmOVERFins));
s = tf('s');
eval(['dampedsys = ',TFstr]);
hh=figure(4);
pzmap(dampedsys);  % plot poles and szeroes
xlim([-2, 2]);ylim([-2, 2]);
title({'s-plane pole-zero plot of damped Spring+Mass ', [' m=' num2str(m)  ', n=' num2str(n)  ', r=' num2str(r) ', f0 Hz=' num2str(freqHz) ]})
disp('         Fig. 4. Poles of damped oscillator. ')
disp(' ')
disp(' ')
[dampedsyspoles,dampedsyszeroes] = pzmap(dampedsys)

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. After making the changes in Q6,')
disp('    looking at the pole-zero plot, is the system critically damped, ')
disp('    underdamped, overdamped, or undamped?')
disp(' ')
disp('Answer:  ???? ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. After making the changes in Q6,')
disp('    looking at the pole-zero plot, what are the coordinates of the  ')
disp('    upper pole in the plot?  And is this a stable system')
disp(' ')
disp('Answer: coords.= ??? +/- j ????     Stable: yes/no')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. After making the changes in Q6,')
disp('    what are the values of the damping factor (delta)  ')
disp('    and natural frequency (omega_0)?')
disp('    (Compare this to your answer in Q9)')
disp(' ')
disp('Answer:    delta = ???     omega_0= ???? ')
disp(' ')
disp(' ')


%  below, make a real-time video of the damped response
pause(2)
fig = uifigure('Visible','on');
selection = uiconfirm(fig,'Click OK to create a real-time movie','Confirm Close','CloseFcn',@(h,e) close(fig));
if selection(1)=='C'   return; end

initImg = imread('mechOscInit','png');
imsiz=size(initImg);
imsizvert=imsiz(1);imsizhor=imsiz(2);
top=45;lspr=203;lmass=106;ldash=lspr;bot=imsizvert-top-lspr-lmass;
imTop=imcrop(initImg,[1, 1, imsizhor, top]);
imSpring=imcrop(initImg,[1, top+1, imsizhor, lspr]);
imMass=imcrop(initImg,[1, top+lspr+1, imsizhor, lmass]);
imDash=imcrop(initImg,[1, top+lspr+lmass+1, imsizhor, ldash]);
imBot=imcrop(initImg,[1, top+lspr+lmass+ldash+1, imsizhor, bot]);

maxz=max(zz)
minz=min(zz)
plotzmiddle=-m*g*n
h100=figure(100);
dt=1/framerate;
vid = VideoWriter('mechOscMovie.mp4','MPEG-4');
vid.Quality = 95;vid.FrameRate=framerate ;
open(vid);
for nt=1:framerate*maxtime;
    tt(nt)=(nt-1)*dt;
    t=tt(nt);
    zz(nt)=double(vpa(subs(dampedzt)));
    z=zz(nt);
    sprnumrows=round((1-(z-plotzmiddle)/abs(maxz-minz))*lspr);
    if sprnumrows<1; sprnumrows=1; end
    if sprnumrows>(ldash+lspr-1); sprnumrows=(ldash+lspr-1); end
    imSpring2=imresize(imSpring,[ sprnumrows imsizhor]);
    imDash2=imresize(imDash,[ ldash+lspr-sprnumrows imsizhor]);
    imOut = cat(1,imTop,imSpring2,imMass,imDash2,imBot);
    h100=image(imOut);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
close(vid);
win=implay('mechOscMovie.mp4')

 

 disp(' you must play the movie to see the real-time behavior, the fames in Fig. 100 were not real time')