! elastic pendulum with force f = F/m  = -k_m (l-l0)

	implicit none
	real(8) :: l, l0, k_m, f, g, t, dt, T_end, x,y,z, vx,vy,vz
	real(8) :: en, en_err, e0
	integer :: i,j
	T_end = 1000d0 ! s
	l0 = 0.1d0 	! m 
	k_m = 33d0 	! Hz**2
	g = 9.81d0	! m/s**2

   do j = 0,2 ! 3 variants will be computed
	dt = 1d-4 ! s
	x = -0.1d0 ; y = 0.2d0; z = 0.1d0 
	vx = 0d0;  vy = 0d0;  vz = -0.1d0
	if (j==2) y = y +0.01d0
	if (j==1) dt = dt*0.5d0
!   E0
	l = sqrt(x*x+y*y+z*z); print*,' init l=',l
	e0 = 0.5d0*( k_m*(l-l0)**2 + vx*vx+vy*vy+vz*vz ) +g*z
 
!  integration in time [0,T_end]
	do i = 1, T_end/dt
		t = i*dt ! end of the timestep time
		l = sqrt(x*x+y*y+z*z)
		f = k_m*(l0-l)/l
		vx = vx + dt*f*x
		vy = vy + dt*f*y
		vz = vz + dt*(f*z - g)
		x = x + dt*vx
		y = y + dt*vy
		z = z + dt*vz
		if (i==T_end/dt) then
!	 Energy computation
		  l  = sqrt(x*x+y*y+z*z)
		  en = 0.5d0*( k_m*(l-l0)**2 + vx*vx+vy*vy+vz*vz ) +g*z 
		  en_err = en/e0-1d0
		  print'(i8,a,4f9.4,e11.3)',i,' i,t,r,l',x,y,z,l,en_err
		endif
	end do
   end do
 end