      program DoublePendulum ! BB Jan 20 2010.
C Purpose: Calculation of eigenvectors for given eigenvalues.
C o2      omega+^2 and omega-^2 eigenvalues with g/l=1.
C op,om   omega+ and omega- eigenvalues
C z=zero  consistency check for correct eigenvalues o2.
      implicit real*8 (a-h,o-z)
      parameter(n=500,tmax=50)
c     parameter(phi0=0.d00,phid0=+0.1d00,psi0=0.d00,psid0=-0.1d00)
      parameter(phi0=0.d00,phid0=+1.0d00,psi0=0.d00,psid0=-1.0d00)
c
      write(6,'(/,"Constants for calculation of eigenvectors:",/)')
      o2=2.d0+sqrt(2.d0)
      op=sqrt(o2)
      a=-2*o2+2
      b=-o2
      c=-o2+1
      z=o2**2-4*o2+2
      write(6,'("a+, b+, c+ =",4F14.6)') a,b,c,z
      psip=-a/b
      chkp=-b/c
c
      o2=2.d0-sqrt(2.d0)
      om=sqrt(o2)
      a=-2*o2+2
      b=-o2
      c=-o2+1
      z=o2**2-4*o2+2
      write(6,'("a-, b-, c- =",4F14.6)') a,b,c,z
      psim=-a/b
      chkm=-b/c
      write(6,'("psi+, psi- =",2F14.6)') psip,psim
      write(6,'("chk+, chk- =",2F14.6)') chkp,chkm
c
      write(6,'(/,"Eigenvalues op,om:",4X,2F9.5)') op,om
      write(6,'("Initial conditions:")')
      write(6,'("phi0,phid0,psi0,psid0:",4F9.5)') 
     &           phi0,phid0,psi0,psid0
      write(6,'(/,"Calculated time series on ts.d.",/)')
      Ap=phi0/2-psi0/2/sqrt(2.d00)
      Am=phi0/2+psi0/2/sqrt(2.d00)
      Bp=(phid0/2-psid0/2/sqrt(2.d00))/op
      Bm=(phid0/2+psid0/2/sqrt(2.d00))/om
      write(6,'("Ap,  Am,   Bp,  Bm:   ",4F9.5)') Ap,Am,Bp,Bm
      open(8,file="ts.d",form='formatted',status='unknown')
      do i=0,n
        t=(i*tmax)/n
        phi= +Ap*cos(op*t)+Bp*sin(op*t)+Am*cos(om*t)+Bm*sin(om*t)
        psi=(-Ap*cos(op*t)-Bp*sin(op*t)+Am*cos(om*t)+Bm*sin(om*t)
     &      )*sqrt(2.d00)
        write(8,'(I6,3F14.6)') i,t,phi,psi
      enddo
      close(8)
c
      stop
      end
