import courses.structuraldynamics.*

k=100.0; m=1.0; xsi=0.1
theSDOF = new sdof(k,m)
theSDOF.setCritDampRatio(xsi)

omegas=[]
om_0=theSDOF.getNaturalFrequency(); println ("om_0= "+om_0);omegas.add(om_0)
om_ext_1=1.5*om_0; a_1=10.0; println ("om_ext_1= "+om_ext_1);omegas.add(om_ext_1)
om_ext_2=8.0*om_0; a_2=100.0; println ("om_ext_2= "+om_ext_2);omegas.add(om_ext_2)
om_ext_3=12.1*om_0; a_3=15.0; println ("om_ext_1= "+om_ext_1);omegas.add(om_ext_3)
om_ext_4=17.0*om_0; a_4=22.0; println ("om_ext_2= "+om_ext_2);omegas.add(om_ext_4)

f={
	(a_1*sin(om_ext_1*it)+a_2*sin(om_ext_2*it)
	+a_3*sin(om_ext_3*it)+a_4*sin(om_ext_4*it))
}

Nt=2**14-1
dt=0.001d
Tot=Nt*dt
/////////////////////////////////////////////////////////////////////////
// Do the FFT of f, F(ω) and compute U(ω)=F(ω)H(ω)
////////////////////////////////////////////////////////////////////////
FFT_f_Plot= new PlotFrame()
FFT_h_Plot= new PlotFrame()
FFT_u_Plot= new PlotFrame()
fext=[]
(0..Nt).each{fext.add(f(it*dt))}

fft= new FastFourierTransformer(DftNormalization.STANDARD)

F_f = fft.transform(fext as double[], TransformType.FORWARD)
Lp=F_f.size()

freq=[]
freqplt=[]
F_fplt=[]

F_h=[]
F_hplt=[]
F_u=[]
F_uplt=[]

Fs=1.0/dt
(0..(Lp)/2).each{
	freqplt.add(2.0*PI*it*Fs/Lp)
	F_fplt.add(F_f[it].abs()/(Lp/2.0))

	freq.add(freqplt[it])
	F_h.add(theSDOF.TransferFunction(freqplt[it]))
	F_hplt.add(F_h[it].abs())

	F_u.add(F_h[it]*F_f[it])
	F_uplt.add(F_u[it].abs())
}

(2..(Lp)/2).each{
	freq.add(-freqplt[-it])
	F_h.add(theSDOF.TransferFunction(-freqplt[-it]))
	F_u.add(theSDOF.TransferFunction(-freqplt[-it])*F_f[-it])
}

omegas.each{
	FFT_f_Plot.vline(it,Color.green)
	FFT_h_Plot.vline(it,Color.green)
	FFT_u_Plot.vline(it,Color.green)
}

FFT_f_Plot.addFunction(new plotfunction(freqplt as double[], F_fplt))
FFT_f_Plot.setMarker(true)
FFT_f_Plot.Title("FFT, F[f(t)]=F(ω) of excitation function")
FFT_f_Plot.show()

FFT_h_Plot.addFunction(new plotfunction(freqplt as double[], F_hplt))
FFT_h_Plot.setMarker(true)
FFT_h_Plot.Title("FFT, F[h(t)]=H(ω) of transfer function")
FFT_h_Plot.show()

pf=new plotfunction(freqplt as double[], F_uplt)
pf.setName("F[u(t)]=F(ω)H(ω)")
FFT_u_Plot.addFunction(pf)
FFT_u_Plot.Title("FFT, F[u(t)]=F(ω)H(ω) of response")
FFT_u_Plot.show()
/////////////////////////////////////////////////////////////////////////
// Do the IFFT of U(ω) to get u(t)
////////////////////////////////////////////////////////////////////////
u = fft.transform(F_u as Complex[], TransformType.INVERSE)
ResponsePlot= new PlotFrame()
uplt=[]
(0..<u.length).each{
	uplt.add(2.0*u[it].re())
}
pf=new plotfunction(dt,uplt as double[])
pf.setName("IFFT[U(ω)]=u(t)")
ResponsePlot.addFunction(pf)

theSDOF.setRHS(f as DoubleFunction)
theSDOF.duhamel(Tot,dt)

pf=new plotfunction(dt,theSDOF.Disp())
pf.setName("Duhamel u(t)")
ResponsePlot.addFunction(pf)
ResponsePlot.setMarker(false)
ResponsePlot.Title("IFFT, F^(-1)[U(ω)]=u(t) of response")
ResponsePlot.setAutoColor(true)
ResponsePlot.makeLegend(true)
ResponsePlot.hline(0.0)
ResponsePlot.show()
/////////////////////////////////////////////////////////////////////////
// Do the FFT of u(t), U(ω) and compare with H(ω)F(ω)
////////////////////////////////////////////////////////////////////////
F_un = fft.transform(theSDOF.Disp(), TransformType.FORWARD)
Lp=F_un.size()

F_unplt=[]

Fs=1.0/dt
(0..(Lp)/2).each{
	F_unplt.add(F_un[it].abs())
}

pf=new plotfunction(freqplt as double[], F_unplt)
pf.setMarkerStyle(1)
pf.setMarkerFill(true)
pf.setName("F[u(t)]")
FFT_u_Plot.addFunction(pf)
FFT_u_Plot.setMarker(true)
FFT_u_Plot.setAutoColor(true)
FFT_u_Plot.makeLegend(true)
FFT_u_Plot.show()

