      program SphTurn ! BB Oct 31 2011.
C Turning points of the spherical pendulum. Problem 22 PHY4936, 2011.
      implicit real*8 (a-h,o-z)
      logical ltest
      parameter(omega1=3.0d00,omega2=1.0d00,n=500)
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)
      ltest=.true.
c
      write(6,'(/,"Turning points of spherical pendulum:",/)')
      write(6,'("Special case of Problem 22 PHY4936, 2011.:",/)')
      pi=acos(-1.0d00)
      E=1.0d00
      write(6,'(" E =",1F6.1)') E
      thmax=pi
      open(8,file="SphTurn.d",form='formatted',status='unknown')
      do i=1,n
        theta=(i*thmax)/n
        y=E/sin(theta)**2+E*cos(theta)/2-E
        write(8,'(2G16.7,I9)') theta,y,i
      enddo
      close(8)
      if(ltest) stop "SphTurn: ltest."
c
      stop
      end
