c
c Linear Systems Introduction
c	Mass-Spring Response to Unit Impulse
c
c	(Solving ODE via Fourier Transforms)
c
c Oct.-23, 2019 Andreas Muenchow
c
	real h,k,m,f,df,arg
	real xi,ff,hmag,hdir
	integer i,nf
	complex h2
	pi = atan2(1.,1.)*4.
c
999	write(6,*)'Enter Friction: 0-None, 1-Tiny, 2-Normal, 3-Large'
	read(5,*)ifriction
c
c Spring Constant (kg/s^2)
c
	k = 10.
c
c Mass m (kg)
c
	m = 2.
c
c Friction r, kg/s
c
	if (ifriction .eq. 0) then
	   xi = 0.
	else if (ifriction .eq. 1) then
	   xi = 0.01
	else if (ifriction .eq. 2) then
	   xi = 0.1
	else if (ifriction .eq. 3) then
	   xi = 1.
	else
	   goto 999
	end if
	r = xi*sqrt(k*m)
c
c Resonance frequency
c
	fres = sqrt(k/m)/(2.*pi)
	write(6,*)'Resonance frequency: ',fres,' s'
	write(6,*)'Friction  : ',r,' kg/s'
c
	df = 0.005
	nf = 100
c
	do 1 i=1,nf
	   f = float(i)*df
	   ff = f/fres
	   arg = 2.*pi*f
	   den = k-m*arg*arg
c
c H Transfer Function s^2/kg
c
	   h = 1./den
	   h2 = (1./k)/(1.-ff*ff+complex(0.,2.*xi*ff))
	   if (h .gt. 0) then
	      write(8,*)ff,h,real(h2),imag(h2)
	   end if
	   x = real(h2)
	   y = imag(h2)
	   hmag = sqrt(x*x+y*y)
	   hdir = atan2(y,x)*180./pi
	      write(9,*)ff,hmag,hdir
1	continue
	stop
	end