c   Simple pendulum october 18, 2006
      data g , Length  /9.80,2./
      x(t)=thetai*cos(omega*t)
      pi=2.*asin(1.)
      omega=sqrt(g/length)
      tscale=1./omega
      omega0=0.
      theta0=1.2*(pi/4.)
      thetai=theta0
      dt=tscale/100.
      tfinal=1.3*(2.*pi)*tscale
      nstep=int(tfinal/dt)
      kp=int(float(nstep)/70.)
      kount=kp
      theta1=theta0+omega0*dt
      print 100 ,0.,theta0, x(0.)
      do 10 i=2,nstep
      t=dt*float(i)
      theta2=2.*theta1-theta0-(dt**2)*omega**2*sin(theta1)
      if (i.eq.kount ) then
      print 100 ,t ,theta2, x(t)
      kount=kount+kp
      endif
      theta0=theta1
      theta1=theta2
10    continue
100   format(1x,'t,theta,x=',3(4x,e11.4))
      stop
      end