! 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) :: x1,y1,z1,l1, 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-3 ! 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
!   push x's forward by dt/2 before start
	z = z +  dt/2*vz
	
!  integration in time
	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, start with rolling position back by dt/2 
		  x1 = x - 0.5d0*dt*vx
		  y1 = y - 0.5d0*dt*vy
		  z1 = z - 0.5d0*dt*vz
		  l1 = sqrt(x1*x1+y1*y1+z1*z1)
		  en = 0.5d0*( k_m*(l1-l0)**2 + vx*vx+vy*vy+vz*vz ) +g*z1 
		  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