Program ORBIT
c
c  This program is listed in Appendix G of both "An Introduction
c  to Modern Astrophysics," Bradley W. Carroll and Dale A. Ostlie,
c  Addison-Wesley Publishing Company, copyright 1996, and "An 
c  Introduction to Modern Stellar Astrophysics," Dale A. Ostlie
c  and Bradley W. Carroll, Addison-Wesley Publishing Company,
c  copyright 1996.
c
      real*8 a, aau, e, time, tyears, dt, LoM, period, Pyears
      real*8 Mstrsun, Mstar, Msun, theta, dtheta, r, x, y
      real*8 G, AU, spyr
      integer*4 n, k, kmax, i
c  define basic constants
      data G, Msun, AU, spyr
     1 / 6.67259d-08, 1.989d33, 1.4960d13, 3.155815d7 /
c  open output file for orbital parameters
      open(unit=10,file='orbit.dat',form='formatted',status='unknown')
c  enter physical parameters for the system
      write(*,*) ' Enter the mass of the parent star (in solar masses):'
      read(*,*) Mstrsun
      write(*,*) ' Enter the semimajor axis of the orbit (in AU):'
      read(*,*) aau
      write(*,*) ' Enter the eccentricity:'
      read(*,*) e
c  calculate orbital period using Kepler's third law (Eq. 2.35)
      Pyears = sqrt(aau**3/Mstrsun)
      write(*,*) ' '
      write(*,100) Pyears
c  enter the number of time steps and the time interval to be printed
      write(*,*) ' '
      write(*,*) ' '
      write(*,*) ' You may now enter the number of time steps to ',
     1 'be calculated and the'
      write(*,*) ' frequency with which you want time steps printed.'
      write(*,*) ' Note that taking too large a time step during the',
     1 ' calculation will'
      write(*,*) ' produce inaccurate results.'
      write(*,*) ' '
      write(*,*) ' Enter the number of time steps desired for the ',
     1 'calculation: '
      read(*,*) n
      write(*,*) ' '
      write(*,*) ' How often do you want time steps printed?'
      write(*,*) '     1 = every time step'
      write(*,*) '     2 = every second time step'
      write(*,*) '     3 = every third time step'
      write(*,*) '                etc.'
      read(*,*) kmax
c  convert to cgs units
      period = Pyears*spyr
      dt = period/float(n-1)
      a = aau*AU
      Mstar = Mstrsun*Msun
c  initialize print counter, angle, and elapsed time
      k = 0
      theta = 0.0d0
      time = 0.0d0
c  print output header
      write(10,200) Mstrsun, aau, Pyears, e
c  start main time step loop
      do 10 i = 1,n
c  increment print counter
          k = k+1
c  use Eq. (2.3) to find r
          r = a*(1.0d0-e**2)/(1.0d0+e*cos(theta))
c  convert to cartesian coordinates and print time and position
c  also print last point to close ellipse
          if (k.eq.1 .or. i.eq.n) then
              x = r*cos(theta)/AU
              y = r*sin(theta)/AU
              tyears = time/spyr
              write(10,300) tyears, x, y
          end if
c  calculate the angular momentum per unit mass L/m (Eq. 2.32)
          LoM = sqrt(G*Mstar*a*(1.0d0-e**2))
c  calculate the next value of theta by combining Eqs. (2.27)
c  and (2.29) (Kepler's second law)
          dtheta = LoM/r**2*dt
          theta = theta + dtheta
c  update the elapsed time
          time = time + dt
c  check print counter then return to the top of the loop
          if (k.eq.kmax) k = 0
   10 continue
      write(*,*) ' '
      write(*,*) ' The calculation is finished and listed in ',
     1  'orbit.dat'
      stop ' '
c  formats
  100 format(1x,' The period of this orbit is ',f10.3,' years')
  200 format(1x,25x,'Elliptical Orbit',//,
     1 26x,'Mstar = ',f10.3,' Mo',/,
     2 26x,'a     = ',f10.3,' AU',/,
     3 26x,'P     = ',f10.3,' yr',/,
     4 26x,'e     = ',f10.3,///,
     5 11x,'t (yr)',14x,'x (AU)',12x,'y (AU)')
  300 format(7x,f10.3,10x,f10.3,8x,f10.3)
      end