% Author: Raymond G. Jacquot
% Modifications by: Eduardo M. O. Lopes
% (Excitações dos casos 5, 6 e 7 desabilitadas.)
% (Ainda em testes!)
% (20/07/22)

clc
clearvars % emol
close all % emol

tol=1e-8; % tolerance % emol

% Graphics position data
arrowx=[0 12 10 10 12]; arrowy=[0 0 .8  -.8 0];
basex=[0 0 -5 -2 -8 -2 10 20 38 65 0]; basey=[0 30 20 17 6 -4 -7 -4 -8 -.5 -.5];
massx=[32 52 52 32 32]; massy=[0 0 30 30 0];
sprx=[0 8 10 14 18 22 24 32]; spry=[22 22 24 20 24 20 22 22];
damp1x=[0 10]; damp1y=[8 8];
damp2x=[24 10 10 24]; damp2y=[10 10 6 6];
damp3x=[17 32]; damp3y=[8 8];
damp4x=[17 17]; damp4y=[10 6];

% Draw the system with labels
figure(1); clf;
axis off; hold on
axis([-10 70 -20 50]); hold on
patch(massx, massy,'r'); hold on
patch(basex,basey,'r'); hold on
plot(sprx,spry,'k','LineWidth',1.5); hold on
plot(damp1x,damp1y,'k','LineWidth',1.5); hold on
plot(damp2x,damp2y,'k','LineWidth',1.5); hold on
plot(damp3x,damp3y,'k','LineWidth',1.5); hold on
plot(damp4x,damp4y,'k','LineWidth',1.5); hold on
plot(arrowx+52*(ones(1,5)),arrowy+33*ones(1,5),'k','LineWidth',1.5); hold on
plot(1.2*arrowx+52*ones(1,5),1.2*arrowy+15*ones(1,5),'k','LineWidth',1.5); hold on;
% text(60,12,'p(t)');
% text(40,15,'m');text(17,5,'c');text(17,18,'k');text(55,30,'u(t)');
% text(5,35,'One Degree-of-Freedom Forced Vibration')
% text(5,-1.8,'Frictionless Plane')
% text(20,-15,'Press Enter to Continue');
text(60,12,'f(t)');
text(40,15,'m');text(17,5,'c');text(17,18,'k');text(55,30,'x(t)');
text(5,35,'Vibração forçada em sistema com 1 gdl') % emol
text(5,-1.8,'Plano sem atrito') % emol
text(12,-15,'>>> Usar a janela de comandos para entrada de parâmetros.'); % emol
% pause(3) % emol
% clc % emol
% disp('  This script considers the forced vibration of a sprung mass');
% disp('  with damping for a variety of forcing functions');
% disp(' ');
% disp('  The system parameters are:');
disp('  Esse programa considera a vibração de um sistema com 1 gdl'); % emol
disp('  para uma variedade de excitações.'); % emol
disp(' '); % emol
disp('  Os parâmetros do sistema são:'); % emol
k=[]; c=[]; m=[];
while isempty(m)
    % m=input('   Mass  m = ');
    m=input('   Massa  m (kg) = '); % emol
end
while isempty(c)
    % c=input('   Viscous Damping c = ');
    c=input('   Coeficiente de amortecimento viscoso c (kg/s) = '); % emol
end
while isempty(k)
    % k=input('   Stiffness k = ');
    k=input('   Constante de rigidez k (N/m) = '); % emol
end
disp(' '); % emol

% Compute natural frequency and damping ratio
wn=sqrt(k/m);
% disp(['    Undamped natural frequency, omegan = ' num2str(wn) ' r/s']);
disp(['Frequência natural = ' num2str(wn) ' rad/s']); % emol
Tn=2*pi/wn;
% disp(['    Undamped natural period Tn = ' num2str(Tn) ' s'])
% disp(['Período natural Tn = ' num2str(Tn) ' s']) % emol
z=c/(2*sqrt(k*m));
% disp(['Damping ratio, zeta = ' num2str(z)]);
disp(['Razão de amortecimento = ' num2str(z)]); % emol
if z < 1 % emol
    wD=wn*sqrt(1-z*z); % emol
    % disp(['Damped natural frequency, omegad = ' num2str(wD) ' r/s']);
    disp(['Frequência natural amortecida = ' num2str(wD) ' rad/s']); % emol
else % emol
    z=z*(1-tol); % emol
    wD=wn*sqrt(1-z*z); % emol
end % emol
disp(' '); % emol

dt=Tn/40;
t=0:dt:400*dt; % emol
Co=cos(wD*dt); Si=sin(wD*dt);
E=exp(-z*wn*dt); H=z/sqrt(1-z*z); J=sqrt(1-z*z);
PHI=zeros(2,2); % reserve space
GAMMA=zeros(2,2); % reserve space 
A=E*(Co+H*Si); B=E*Si/wD; Ap=-E*((wn/J)*Si); Bp=E*(Co-H*Si);
PHI=[A B;Ap Bp];% fill the transition matrix
C=(1/k)*(((2*z)/(wn*dt))+E*((((1-2*z*z)/(wD*dt))-H)*Si-(1+((2*z)/(wn*dt)))*Co));
D=(1/k)*(1-((2*z)/(wn*dt))+E*(((2*z*z-1)/(wD*dt))*Si+((2*z/(wn*dt)))*Co));
Cp=(1/k)*(-(1/dt)+E*(((wn/J)+(H/dt))*Si+(1/dt)*Co));
Dp=(1/(k*dt))*(1-E*(H*Si+Co));
GAMMA=[C D;Cp Dp];% fill the input distribution matrix

% Choose the excitation
% disp('  The forcing function takes one of the following forms')
% disp('        0- No Force, Initial displacement and velocity only')
% disp('        1- Sinusoidal Force starting at t = 0, zero initial conditions')
% disp('        2- Step Function at t = 0, zero initial conditions')
% disp('        3- Rectangular Pulse at t = 0, zero initial conditions')
% disp('        4- Half Sine Pulse at t = 0, zero initial conditions')
% disp('        5- Triangular Pulse Rising Linearly with a Sudden Drop at Td,')
% disp('                       zero initial conditions')
% disp('        6- Triangular Pulse with a Sudden Rise and a Linear Decay to Td,')
% disp('                       zero initial conditions')
% disp('        7- Symmetric Triangular Pulse, zero initial conditions')
disp('  A excitação pode assumir as seguintes formas:') % emol
disp('        0 - Apenas deslocamento e velocidade iniciais')
disp('        1 - Função harmônica com início em t = 0, condições iniciais nulas') % emol
disp('        2 - Função degrau em t = 0, condições iniciais nulas') % emol
disp('        3 - Pulso retangular em t = 0, condições iniciais nulas') % emol
disp('        4 - Pulso meio-seno em t = 0, condições iniciais nulas') % emol
% disp('        5 - Pulso triangular com crescimento linear e súbita queda em Td,') % emol
% disp('                condições iniciais nulas') % emol
% disp('        6 - Pulso triangular com crescimento repentino e decréscimo linear') % emol
% disp('                para Td, condições iniciais nulas') % emol
% disp('        7 - Pulso triangular simétrico, condições iniciais nulas') % emol
disp(' '); 

% Comandos para que apenas as opções de 0 a 7 possam ser digitadas.  % emol
% Caso contrário, a opção de excitação continuará a ser apresentada. % emol
% fopt=-1;
% while ((fopt~=0)&&(fopt~=1)&&(fopt~=2)&&(fopt~=3)&&(fopt~=4)&&(fopt~=5)&&(fopt~=6)&&(fopt~=7))||(isempty(fopt)) 
%     fopt=input('    Option Number = ');
% end

% Comando para entrada da opção de excitação. % emol
% Caso uma opção inválida seja digitada,      % emol
% o programa avisará e será encerrado.        % emol
fopt=input('    Opção de excitação = '); % emol
disp(' '); % emol

u0=[];v0=[];
switch fopt 
    case 0
       while isempty(u0)
            % u0=input('    Initial displacement u(0) = ');
            u0=input('    Deslocamento inicial x(0) (m) = '); % emol
       end
       while isempty(v0)
            % v0=input('    Initial velocity du(0)/dt = ');
            v0=input('    Velocidade inicial v(0) (m/s) = '); % emol
       end
       x=[u0; v0];U=[u0;v0];
    case 1
      A1=[];omega=[];
      while isempty(A1)
        % A1=input('    Amplitude of Sinusoidal Force = ');
        A1=input('    Amplitude da função harmônica (N) = '); % emol
      end
      while isempty(omega)||omega<=0
        % omega=input('    Radian Frequency of Sinusoidal Force = ');
        omega=input('    Frequência da função harmônica (rad/s) = '); % emol
      end
      f=A1*sin(omega*t);
      f1=[f(1,2:401) A1*sin(omega*(t(1,401)+dt))];   
    case 2 % emol
        A1=[];
        while isempty(A1)
            % A1=input('    Height of Step Function = ');
            A1=input('    Altura da função degrau (N) = ');
        end
        f=A1*ones(1,401);
        f1=f;
    case 3
        A1=[];Td=[];
        while isempty(A1)
            % A1=input('    Height of Pulse = ');
            A1=input('    Altura do pulso retangular (N) = ');
        end
        while isempty(Td)||Td<=0 % emol
            % Td=input('    Duration of Pulse in sec. = ');
            Td=input('    Duração do pulso retangular (s) = ');
        end
        nd=fix(Td/dt);
        f=[A1*ones(1,nd+1) zeros(1,401-(nd+1))];
        f1=f;
    case 4
        A1=[];Td=[];
        while isempty(A1)
            % A1=input('    Amplitude of Sine Pulse = ');
            A1=input('    Amplitude do pulso senoidal (N) = ');
        end
        while ((isempty(Td))||(Td<=0)) % emol
            % Td=input('    Duration of Pulse in sec. = ');
            Td=input('    Duração do pulso senoidal (s) (>= 0.02 s) = ');
        end
        nd=fix(Td/dt);
        t1=0:dt:nd*dt; % emol
        f=[A1*sin(pi*t1/Td) zeros(1,401-(nd+1))];
        f1=[f(1,2:(nd)) zeros(1,401-nd+1)];
%     case 5(cases 5 to 7 are disabled) % emol
%         A1=[];Td=[];
%         while isempty(A1)
%             A1=input('    Height of Triangular Pulse = ');
%             % A1=input('    Height of Triangular Pulse = ');
%         end
%         while ((isempty(Td))||(Td<=0))
%             Td=input('    Duration of Pulse in sec. = ');
%             % Td=input('    Duration of Pulse in sec. = ');
%         end
%         nd=fix(Td/dt);
%         t1=0:dt:nd*dt; % emol
%         f=[ A1*t1/Td zeros(1,401-(nd+1))];
%         f1=[f(1,2:nd+1) zeros(1,401-(nd))];
%     case 6 (cases 5 to 7 are disabled) % emol
%         A1=[];Td=[];
%         while isempty(A1)
%             A1=input('    Height of Triangular Pulse = ');
%             % A1=input('    Height of Triangular Pulse = ');
%         end
%         while ((isempty(Td))||(Td<=0))
%             Td=input('    Duration of Pulse in sec. = ');
%             % Td=input('    Duration of Pulse in sec. = ');
%         end
%         nd=fix(Td/dt);
%         t1=0:dt:nd*dt; % emol
%         f=[A1*(ones(1,nd+1)-t1/Td) zeros(1,401-(nd+1))];
%         f1=[f(1,2:nd+1) zeros(1,401-(nd))];
%     case 7 (cases 5 to 7 are disabled) % emol
%         A1=[];Td=[];
%         while isempty(A1)
%             A1=input('    Height of Triangular pulse = ');
%             % A1=input('    Height of Triangular pulse = ');
%         end
%         while ((isempty(Td))||(Td<=0))
%             Td=input('    Duration of Pulse in sec. = ');
%             % Td=input('    Duration of Pulse in sec. = ');
%         end
%         nd=fix(Td/dt);
%         t1=0:dt:nd*dt; % emol
%         if fix(nd/2)==nd/2 % if true, nd is even
%             n1=0:1:(nd/2); % emol
%             f3=A1*n1/((nd/2));
%             f=[f3 fliplr(f3(1,1:(nd/2))) zeros(1,401-(nd+1))];
%         else % nd is odd
%             n1=0:1:(nd-3)/2; % emol
%             f3=A1*n1*2/(nd-1);
%             f=[f3 A1 fliplr(f3) zeros(1,401-(nd))];
%         end
%         f1=[f(1,2:nd+1)  zeros(1,401-(nd))]; 
    otherwise % emol
        disp('Opção inválida. O programa será encerrado.') % emol
        disp(' ') % emol
        return % emol
end
disp(' '); % emol
if fopt~=0
    Fmax=max(max(f));
    Fmin=min(min(f));
    if Fmax==Fmin
         Fmin=0;
         F1=-.1;F2=1.1*Fmax;
    end
    if (Fmax==Fmin)&&(Fmax==0)
        F1=.1;F2=-.1;
        Frange=.2;
    else
        Frange=abs(Fmax-Fmin);
        F1=1.1*Fmax;
        F2=Fmax-1.1*Frange;
    end
end

if fopt==2||fopt==3||fopt==4||fopt==5||fopt==6||fopt==7 % emol
    figure(2); clf
    axis([-10*dt t(1,401) F2 F1]); hold on
    plot([dt*fliplr(-[1:1:10]) t],[zeros(1,10) f])
%     xlabel('Time t, sec.')
%     ylabel('Force p(t)')
%     text(t(1,200) ,F1-.2*Frange,'Press Enter to Continue')
    xlabel('Tempo (s)') % emol
    ylabel('Excitação (N)') % emol
    box on % emol
    grid on % emol
    pause(3) % emol
elseif fopt==1 % emol
    figure(2); clf % emol
    axis([0 t(1,401) Fmin Fmax]); hold on % emol
    plot(t,f) % emol
    xlabel('Tempo (s)') % emol
    ylabel('Excitação (N)') % emol 
    box on % emol
    grid on % emol
    pause(3) % emol
end
if fopt==0
    x=[u0;v0];U=[u0;v0];
else 
    x=[0;0];U=[0;0];
end

for k1=1:400
    if fopt~=0
        x=PHI*x+GAMMA*[f1(1,k1);f1(1,k1+1)];
    elseif fopt==0
        x=PHI*x;
    end
    U=[U x];
end
Umax=max(max(abs(U(1,:))));Umin=min(min(U(1,:)));
Udotmax=max(max(abs(U(2,:))));Udotmin=min(min(U(2,:)));

figure(3); clf % emol
axis([0 t(1,401) Umin Umax]); hold on
plot(t,U(1,:))
% xlabel('Time t, sec.')
% ylabel(' Displacement u(t)');
% text(t(1,200),Umax/2,'Press Enter to Continue');
xlabel('Tempo (s)') % emol
ylabel('Deslocamento (m)'); % emol
% text(t(1,200),Umax/2,'Pressione Enter para continuar.'); % emol
box on
grid on % emol
pause(3) % emol

figure(4); clf % emol
Udotmax=max(max(abs(U(2,:))));
axis([0 t(1,401) Udotmin Udotmax]); hold on
plot(t,U(2,:));
% xlabel('Time t, sec.')
% ylabel(' Velocity du(t)/dt');
% text(t(1,200),Umax/2,'Press Enter to Continue');
xlabel('Tempo (s)') % emol
ylabel('Velocidade (m/s)'); % emol
% text(t(1,200),Umax/2,'Pressione Enter para continuar.'); % emol
box on
grid on % emol
pause(3) % emol

figure(1); clf
axis off
axis([-10 70 -20 50]); hold on;
Mass=patch(massx, massy,'r'); hold on % emol
patch(basex,basey,'r'); hold on % emol
Spring=plot(sprx,spry,'k','LineWidth',1.5); hold on % emol
plot(damp1x,damp1y,'k','LineWidth',1.5); hold on % emol
plot(damp2x,damp2y,'k','LineWidth',1.5); hold on % emol
Damp3=plot(damp3x,damp3y,'k','LineWidth',1.5); hold on % emol
Damp4=plot(damp4x,damp4y,'k','LineWidth',1.5); hold on % emol

% text(5,-1.8,'Frictionless Plane'); hold on
% text(5,35,'One Degree-of-Freedom Free Vibration');
text(5,-1.8,'Plano sem atrito'); hold on % emol
text(5,35,'Vibração em sistema com 1 gdl'); hold on % emol
% txt1='Press Enter to Set Initial Displacement';
% txt1='Pressione Enter para ajustar o deslocamento inicial.'; % emol
txt1='Ajuste de deslocamento inicial.'; % emol
Txt1=text(10,-15,txt1); hold on
pause(3) % emol

% U1=U(1,:)*10/Umax;
U1=U(1,:)*6/Umax; % it should be 6 for motion within the cylinder % emol
figure(1); clf
axis off; hold on
axis([-10 70 -20 50]); hold on;
patch(basex,basey,'r'); hold on % emol
% patch(wallx,wally,'r','EraseMode','xor'); hold on
plot(damp1x,damp1y,'k','LineWidth',1.5); hold on % (rod of cylinder) % emol
plot(damp2x,damp2y,'k','LineWidth',1.5); hold on % (body of cylinder) % emol
Mass=patch(massx+U1(1,1)*ones(1,5), massy,'r'); hold on % emol
Damp3=plot(damp3x+U1(1,1)*ones(1,2),damp3y,'k','LineWidth',1.5); hold on % (rod of piston) % emol
Damp4=plot(damp4x+U1(1,1)*ones(1,2),damp4y,'k','LineWidth',1.5); hold on % (base of piston) %emol
Spring=plot(sprx+[0 0 .125*U1(1,1) .375*U1(1,1) .625*U1(1,1) .875*U1(1,1) U1(1,1) U1(1,1)],spry,'k','LineWidth',1.5); hold on % emol

% txt1=['Press Enter to Animate     umax =' num2str(Umax)];
% text(5,35,'One Degree-of-Freedom Free Vibration');
% txt1=['Pressione Enter para animar (umax = ' num2str(Umax) ').']; % emol
text(5,35,'Vibração em sistema com 1 gdl'); % emol
hold on
% text(5,-1.8,'Frictionless Plane');
text(5,-1.8,'Plano sem atrito'); % emol
hold on

% txt1=['Animação dentro de 3s (x_{máx} = ' num2str(Umax) 'm).']; % emol
% Txt1=text(10,-15,txt1);

txt1=['x_{máx} = ' num2str(Umax) 'm.']; % emol
Txt1=text(10,-15,txt1); % emol

% Txt2=text(10,-18,['t = ' num2str(t(1,1)) ' sec.']); hold on
Txt2=text(10,-18,['t = ' num2str(t(1,1)) ' s.']); hold on % emol
pause(3) % emol
for k=1:400
    set(Mass,'Xdata',massx+U1(1,k)*ones(1,5));
    set(Damp3,'Xdata',damp3x+U1(1,k)*ones(1,2));
    set(Damp4,'Xdata',damp4x+U1(1,k)*ones(1,2));
    set(Spring,'Xdata',sprx+[0 0 .125*U1(1,k) .375*U1(1,k) .625*U1(1,k) .875*U1(1,k) U1(1,k) U1(1,k)]);
    % set(Txt2,'String',['t = ' num2str(t(1,k)) ' sec.']);
    set(Txt2,'String',['t = ' num2str(t(1,k)) ' s.']); % emol
    pause(.04)
end