program temp
c      Analysis of Kapitza pendulum. a*Cos(om*t) is the base excitation,
c      g is the gravitational body force, l is the length of the 
c      pendulum, th0, v0 are initial theta and $\dot{\theta}$.
       integer i,j,nsteps,niter,printstep
       real*8 a,om,g,l,th0,v0,thn,thnp1k,thmk,vn,vm,vnp1,thmdelta,
     & tn,tnp1,thnp1delta,tdelta,fdelta1,fdelta2,xi1,tprint,
     & pi,m11,m12,m21,m22,detm,tend,uf(300000),vf(300000),
     & time(300000)
       open(11,file='shape',status='unknown')
c      ************************************************
       pi=4.0d0*datan(1.0d0)
       xi1=dsqrt(0.6d0)
       g=9.81d0
       l=2.0d0
       a=2.0d0
       om=3.0d0
       th0=3.00d0
       v0=0.00d0
       tend=10.0d0
       niter=10
       tdelta=0.001d0
       tprint=0.1d0
c      prints after every tprint times
       thmdelta=0.0d0
       thnp1delta=0.0d0
       nsteps=nint(tend/tdelta)
       printstep=nint(tprint/tdelta)
       if(printstep.lt.1)printstep=1
       if(nint(1.0d0*nsteps/printstep).gt.300000)then
        write(11,*) 'Increase size of arrays uf, vf etc.'
        stop
       endif
       print*,'nsteps,niter ',nsteps,niter
       time(1)=0.0d0
       thn=th0
       vn=v0
       uf(1)=th0
       vf(1)=v0
       thmk=th0
       thnp1k=th0
       do 1 i=2,nint(1.0d0*nsteps/printstep)+1
        time(i)=printstep*(i-1)*tdelta
1      continue
       do 100 i=1,nsteps
       tn=(i-1)*tdelta
       tnp1=i*tdelta
       do 10 j=1,niter
      m11= (8*(-16*l - tdelta**2*Cos(thmk)*
     &  (3*g - 2*a*om**2*Cos((om*(tn + tnp1))/2.))))/9. + 
     & (5*((-1 - xi1)*(-4*l*(-4 + 12*xi1) + 
     &    tdelta**2*(1 - xi1**2)*
     & (3*g - 2*a*om**2*Cos((om*(tn + tnp1 + tn*xi1 - tnp1*xi1))/2.))*
     &     Cos(thmk - thmk*xi1**2 - 
     &       (xi1*(-thn + thnp1k - (thn + thnp1k)*xi1))/2.)) + 
     & (-1 + xi1)*(-4*l*(-4 - 12*xi1) + 
     &    tdelta**2*(1 - xi1**2)*
     & (3*g - 2*a*om**2*Cos((om*(tn + tnp1 - tn*xi1 + tnp1*xi1))/2.))*
     &     Cos(thmk - thmk*xi1**2 + 
     &       (xi1*(-thn + thnp1k + (thn + thnp1k)*xi1))/2.))))/9.

         m12=(64*l)/9. + (5*((-1 - xi1)*(-4*l*(2 - 3*xi1) - 
     &    (tdelta**2*(1 - xi1)*xi1*
     &       (3*g - 2*a*om**2*
     &          Cos((om*(tn + tnp1 + tn*xi1 - tnp1*xi1))/2.))*
     &       Cos(thmk - thmk*xi1**2 - 
     &         (xi1*(-thn + thnp1k - (thn + thnp1k)*xi1))/2.))/2.) + 
     & (-1 + xi1)*(-4*l*(2 + 3*xi1) + 
     &    (tdelta**2*xi1*(1 + xi1)*
     &       (3*g - 2*a*om**2*
     &          Cos((om*(tn + tnp1 - tn*xi1 + tnp1*xi1))/2.))*
     &       Cos(thmk - thmk*xi1**2 + 
     &         (xi1*(-thn + thnp1k + (thn + thnp1k)*xi1))/2.))/2.)))/9.

        m21= (8*(-16*l + tdelta**2*Cos(thmk)*
     &  (-3*g + 2*a*om**2*Cos((om*(tn + tnp1))/2.))))/9. + 
     & (5*((1 - xi1)*(-16*l + 48*l*xi1 + 
     &    (tdelta**2*(2 - 2*xi1**2)*
     &       (-3*g + 2*a*om**2*
     &          Cos((om*(tn + tnp1 + tn*xi1 - tnp1*xi1))/2.))*
     &       Cos((2*thmk - (-thn + thnp1k)*xi1 + 
     &           (-2*thmk + thn + thnp1k)*xi1**2)/2.))/2.) + 
     & (1 + xi1)*(-16*l - 48*l*xi1 + 
     &    (tdelta**2*(2 - 2*xi1**2)*
     &       (-3*g + 2*a*om**2*
     &      Cos((om*(tn + tnp1 - tn*xi1 + tnp1*xi1))/2.))*
     &       Cos((2*thmk + (-thn + thnp1k)*xi1 + 
     &     (-2*thmk + thn + thnp1k)*xi1**2)/2.))/2.)))/9.

       m22=(64*l)/9. + (5*((1 - xi1)*(8*l - 12*l*xi1 + 
     &    (tdelta**2*(-xi1 + xi1**2)*
     &       (-3*g + 2*a*om**2*
     &          Cos((om*(tn + tnp1 + tn*xi1 - tnp1*xi1))/2.))*
     &       Cos((2*thmk - (-thn + thnp1k)*xi1 + 
     &           (-2*thmk + thn + thnp1k)*xi1**2)/2.))/2.) + 
     & (1 + xi1)*(8*l + 12*l*xi1 + 
     &    (tdelta**2*(xi1 + xi1**2)*
     &       (-3*g + 2*a*om**2*
     &          Cos((om*(tn + tnp1 - tn*xi1 + tnp1*xi1))/2.))*
     &       Cos((2*thmk + (-thn + thnp1k)*xi1 + 
     &           (-2*thmk + thn + thnp1k)*xi1**2)/2.))/2.)))/9.

        detm=m11*m22-m12*m21

      fdelta1= (-8*(8*l*(-2*thmk + thn + thnp1k) - 
     & tdelta**2*(3*g - 2*a*om**2*Cos((om*(tn + tnp1))/2.))*
     & Sin(thmk)))/9. - 
     & (5*((-1 - xi1)*(-4*l*(2*(-2*thmk + thn + thnp1k) - 
     &       3*(-4*thmk + 3*thn + thnp1k + tdelta*vn)*xi1) + 
     &    tdelta**2*(3*g - 2*a*om**2*
     &        Cos((om*(tn + tnp1 + tn*xi1 - tnp1*xi1))/2.))*
     &     Sin(thmk - thmk*xi1**2 - 
     &       (xi1*(-thn + thnp1k - (thn + thnp1k)*xi1))/2.)) + 
     & (-1 + xi1)*(-4*l*(2*(-2*thmk + thn + thnp1k) + 
     &       3*(-4*thmk + 3*thn + thnp1k + tdelta*vn)*xi1) + 
     &    tdelta**2*(3*g - 2*a*om**2*
     &        Cos((om*(tn + tnp1 - tn*xi1 + tnp1*xi1))/2.))*
     &     Sin(thmk - thmk*xi1**2 + 
     &   (xi1*(-thn + thnp1k + (thn + thnp1k)*xi1))/2.))))/9.

        fdelta2=(-8*(8*l*(-2*thmk + thn + thnp1k) + 
     & tdelta**2*(-3*g + 2*a*om**2*
     & Cos((om*(tn + tnp1))/2.))*Sin(thmk)))/9.
     & - (5*((1 - xi1)*(8*l*(-2*thmk + thn + thnp1k) - 
     &    12*l*(-4*thmk + 3*thn + thnp1k + tdelta*vn)*xi1 + 
     &    tdelta**2*(-3*g + 2*a*om**2*
     &        Cos((om*(tn + tnp1 + tn*xi1 - tnp1*xi1))/2.))*
     &     Sin((2*thmk - (-thn + thnp1k)*xi1 + 
     &         (-2*thmk + thn + thnp1k)*xi1**2)/2.)) + 
     & (1 + xi1)*(8*l*(-2*thmk + thn + thnp1k) + 
     &    12*l*(-4*thmk + 3*thn + thnp1k + tdelta*vn)*xi1 + 
     &    tdelta**2*(-3*g + 2*a*om**2*
     &        Cos((om*(tn + tnp1 - tn*xi1 + tnp1*xi1))/2.))*
     &     Sin((2*thmk + (-thn + thnp1k)*xi1 + 
     &         (-2*thmk + thn + thnp1k)*xi1**2)/2.))))/9.

       thmdelta=(m22*fdelta1-m12*fdelta2)/detm
       thnp1delta=(-m21*fdelta1+m11*fdelta2)/detm
       thmk=thmk+thmdelta
       thnp1k=thnp1k+thnp1delta
       if(dsqrt(thmdelta**2+thnp1delta**2).lt.1.0d-10)go to 20
10     continue
       print*,'Failed to converge'
       stop
20     continue
       vm=(4*thmk + thnp1k -5*thn - tdelta*vn)/(2*tdelta)
       vnp1=4*(thnp1k+thn-2*thmk)/tdelta + vn
       vn=vnp1

       if(i.ge.printstep)then
       if(mod(i,printstep).eq.0)then
       uf(nint(1.0d0*i/printstep)+1)=thnp1k
       vf(nint(1.0d0*i/printstep)+1)=vnp1
       endif
       endif
       thn=thnp1k
100    continue

       do 200 i=1,nint(1.0d0*nsteps/printstep)+1
       write(11,500) time(i),uf(i),vf(i)
200    continue
500    format(1x,e14.7,1x,e14.7,1x,e14.7)
       end