c      FOUCALT pendulum
      real lat
      data T ,lat /8.64E4, 45./
      ax(y0,y1,x)=2.*omegae*sin(lat)*(y1-y0)/dt - wp**2*x
      ay(x0,x1,y)=-2.*omegae*sin(lat)*(x1-x0)/dt - wp**2*y
      pi=2.*asin(1.)
      omegae=2.*pi/T
      wp=10.*omegae
      lat=lat*pi/180.
      x0=1.
      x1=x0
      y0=0.
      y1=y0
      tscale=amin1(T,1./wp)
      dt =tscale/200
      tf=20.*tscale
      nstep=int(tf/dt)
      kp=int(float(nstep)/100.)
      kount=kp
      print*,'      x                      y  '
      print*,'    '
      print 100 , x0,y0
      do 10 i=2,nstep
      t=dt*float(i)
      x2=2.*x1-x0 + (dt**2)*ax(y0,y1,x1)
      y2=2.*y1-y0 +(dt**2)*ay(x0,x1,y1)
      if(i.eq.kount)then
      print 100 , x2,y2
      kount=kount+kp
      endif
      x0=x1
      x1=x2
      y0=y1
      y1=y2
10    continue
100   format(1x,'x,y=',2(4x,e10.3))
      stop
      end