#Exemplo - Exercício 3.3

clear

F0 = input("Digite o valor da amplitude da força (F0) em N: ");
w = input("Digite o valor da frequência angular da força harmônica (rad/s): ");
k = input("Digite o valor da rigidez da mola k em N/m: ");
m = input("Digite o valor da massa m em kg: ");
x0 = input("Digite o valor do deslocamento inicial x0 em m: ");
v0 = input("Digite o valor da velocidade inicial v0 em m/s: ");
wn = sqrt(k/m);
A = sqrt(x0^2+(v0/wn)^2);
if x0 == 0
  phi = pi;
else
  phi = atan(v0/(x0*wn));
end

delta_st = F0/k;
r = w/wn;

for i = 1: 1001
t(i) = (2*(1/wn)) * (i-1)/100;
F(i) = cos(w*t(i)); 
if (r <=1)
  if (r == 1)
    x(i) = x0*cos(wn*t(i))+(v0/wn)*sin(wn*t(i))+((delta_st*wn*t(i))/2)*sin(wn*t(i));
  else
  x(i) = A*cos(wn*t(i)-phi)+(delta_st/(1-r^2))*cos(w*t(i));
  z(i) = (x0 - (F0/(k-m*w^2)))*cos(wn*t(i)) +(v0/wn)*sin(wn*t(i))+ (F0/(k-m*w^2))*cos( w*t(i));
  end
else
  x(i) = A*cos(wn*t(i)-phi)-(delta_st/(1-r^2))*cos(w*t(i));
  z(i) = (x0 - (F0/(k-m*w^2)))*cos(wn*t(i)) +(v0/wn)*sin(wn*t(i))+ (F0/(k-m*w^2))*cos( w*t(i));
end
end

plot(t,x)
grid on