##leer librerias
library(optimbase)

T=1 #periodo de la fuerza
wt=2*pi/T
t=seq(0,5*T,T/100)

a0=0;

#construir Fuerza
Ft=a0/2*ones(length(t),1)
for (n in 1:100){
an=-8/(((2*n-1)*pi)^2)
Ft=Ft+an*cos((2*n-1)*wt*t)
}

#graficar Fuerza
plot(t,Ft,type="l",xlab="Tiempo (s)",ylab="Fuerza (N)")
grid(NULL,NULL)


#parametros del sistema
k=1e4
m=100
wn=sqrt(k/m)
chi=0.01

#construir respuesta
Xt=a0/(2*k)*ones(length(t),1)
for (n in 1:100){
an=-8/(((2*n-1)*pi)^2)
thetan=atan((2*chi*n*wt)/(wn^2-(n*wt)^2))
xcn=(an/m)/sqrt((wn^2-(n*wt)^2)^2 + (2*chi*wn*n*wt)^2)*cos(n*wt*t-thetan)
Xt=Xt+xcn
}

#graficar Respuesta
dev.new()
plot(t,Xt,type="l",xlab="Tiempo (s)",ylab="Respuesta (m)")
grid(NULL,NULL)