program klatno
	real L
	call initial(theta1,v1,theta2,v2,w2,t,dt,dt2,L)
	do 10 i = 1,1000
	call output(theta1,v1,theta2,v2,t,x1,y1,x2,y2)
	call linear(theta1,v1,w2,t,dt,dt2,x1,y1,L)
	call nonlinear(theta2,v2,w2,t,dt,dt2,x2,y2,L)
	t = t + dt
10	continue
      stop
	end
*
	subroutine initial(theta1,v1,theta2,v2,w2,t,dt,dt2,L)
	real L
	t = 0
	theta1 = 0.25
	theta2 = theta1
*     pocetna brzina (m/s)
      v1 = 0.
	v2 = v1
*	 g/L  
      w2 = 9.
	L = 10./w2
      dt = 0.01
	dt2 = 0.5 * dt
	write (8,15)
15    format(3x,'vreme',10x,'theta1',10x,'brzina1',10x,'theta2',
     +   10x,'brzina2')
	return
	end
*
      subroutine linear(theta1,v1,w2,t,dt,dt2,x1,y1,L)
	real L
	a = - w2 * theta1
	v1mid = v1 + a * dt2
	theta1mid = theta1 + v1 * dt2
	amid = - w2 * theta1mid
	v1 = v1 + amid * dt
	theta1 = theta1 + v1mid * dt
	x1 = L * sin(theta1)
	y1 = L * (1 - cos(theta1))
	return
	end 
*	
      subroutine nonlinear(theta2,v2,w2,t,dt,dt2,x2,y2,L)
	real L
	a = - w2 * sin(theta2)
	v2mid = v2 + a * dt2
	theta2mid = theta2 + v2 * dt2
	amid = - w2 * sin(theta2mid)
	v2 = v2 + amid * dt
	theta2 = theta2 + v2mid * dt
	x2 = L * sin(theta2)
      y2 = L * (1 - cos(theta2))
	return
	end 
	subroutine output(theta1,v1,theta2,v2,t,x1,y1,x2,y2)
	write (8,*) t,theta1,v1,theta2,v2
	write (9,*) t,x1,y1,x2,y2
      return
	end