!*********************************************************************
! Filename:      hw7.f
! Author:        Eddie Baron <baron@ou.edu>
! Created at:    Tue Sep  1 13:51:33 2009
! Modified at:   Wed Sep  2 13:27:46 2009
! Modified by:   Eddie Baron <baron@ou.edu>
! Version:       
! Description:   Solution to HW 7 problem particle oscillating through
!                the Earth
!*********************************************************************
      program hw7
      implicit none
      integer, parameter :: dp=kind(0.d0)
      real(kind=dp),parameter :: me=5.9742d27,re=6378.1d5,G=6.674d-8
      real (kind=dp) :: rho,pi,r,v,a,t
      real (kind=dp) :: dt,period
      real (kind=dp) :: rp1,rm1,am1,vm1
      real (kind=dp) :: pe,ke,te,omega
      integer :: i
!      
      pi = acos(-1.d0)
      rho = me/(4.d0*pi/3.d0*re**3)
      r = re
      v = 0.d0
      a = -4.d0/3.d0*g*pi*rho*r
      omega = sqrt(4.d0/3.d0*g*pi*rho)
      period = 2.d0*pi/omega
      dt = 1.d-4*period
      t = 0.d0
      pe = 0.5d0*(omega*r)**2
      ke = 0.5d0*v**2
      te = pe+ke
!--
!-- first step Leap frog
!--
      rp1 = r + v*dt + 0.5d0*a*dt**2
      rm1 = r
      am1 = a
      vm1 = v
      r = rp1
      a = -4.d0/3.d0*g*pi*rho*r
      v = 0.5d0*(a+am1)*dt
      t = t + dt
      pe = 0.5d0*(omega*r)**2
      ke = 0.5d0*v**2
!--
!-- ok now we can do Verlet
!--
      do i=1,10000
      rp1 = 2d0*r - rm1 + a*dt**2
      v = (rp1 - rm1)/(2.d0*dt)
      rm1 = r
      am1 = a
      r = rp1
      a = -4.d0/3.d0*g*pi*rho*r
      t = t + dt
      pe = 0.5d0*(omega*r)**2
      ke = 0.5d0*v**2
      write(*,'(i0,7es14.6)')i,t,r,v,a,((pe+ke)-te)/te,2.d0*pe/
     &   (omega*re*cos(omega*t))**2,2.d0*ke/(omega*re*sin(omega*t))**2
      enddo
      end program hw7