program centralforce
C
C   CWJ SDSU April 2005
C
C   computes motion of a particle in a central force
C   integrates differential eqns using 4th order 
C   Runge-Kutta for 2 independent variable
C  (found in rungekuttalib.f )

      implicit none

      real r, theta

      real u,p  !  convenient variables; u = 1/r
                    !     p = du/dtheta

      real L      ! angular momentum (explicitly conserved)
      common/centralparams/L      

      real xn
      common/forceparam/xn

      real dtheta		! steps in theta	
      integer Nsteps	! # of steps taken

      real x,y		! convert r, theta at end

      integer n
      character*15 filename
      integer ilast

      print*,' Enter index n of force -1/r**n '
      print*,' (May be real number) '
      read*,xn

      print*,' Enter ang momentum L '
      read*,L

      print*,' Enter initial r '
      read*,r
      u = 1./r
      p = 0

      print*,' Enter dtheta '
      read*,dtheta

      print*,' Enter # of steps '
      read*,nsteps


      print*,' Enter name of .dat file for output '
      read*,filename
      ilast = index(filename,' ')-1

      open(unit=17,name=filename(1:ilast)//'.dat',status='unknown')
  
      do n = 1,nsteps
         call rungekutta2(dtheta,theta,u,p)
         r = 1./u
         print*,r,theta
         x = r*cos(theta)
         y = r*sin(theta)
         write(17,*)x,y
      enddo

      

      end

CCCCCCCCCCCCCCCCCCCCCCCCCC

      real function f(u,p,theta)
C
C  a function for the coupled first-order ODEs
C
C  u'(theta) = f(u,p,theta)
C  p'(theta) = g(u,p,theta)
C
C  to be fed into Runge-Kutta
C
C  for the central-force problem,
C  u = 1/r, and p = u'
C
C
      implicit none

      real u
      real p
      real theta

      f = p
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCC

      real function g(u,p,theta)
C
C  a function for the coupled first-order ODEs
C
C  u'(theta) = f(u,p,theta)
C  p'(theta) = g(u,p,theta)
C
C  to be fed into Runge-Kutta
C
C  for the central-force problem,
C  u = 1/r, and p = u'
C
C  Here I have an additional problem: I need to 
C  pass parameters MASS and L to this 
C  subroutine. Rather than passing through
C  subroutine and function calls, I will use 
C a "common block" which is Fortran 77's version
C of global variables
C

      implicit none

      real u
      real p
      real theta

C..............."global" variables.............

      real L
      common/centralparams/L

C............the force itself

      real force

      g = - u - 1./L**2/u**2*Force(1./u)

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCC

      real function force(r)
      implicit none

      real r
      real xn
      common/forceparam/xn

      Force = - 1./r**xn

      return
      end