      program phy4241orbits ! Bernd Berg, Feb 17 2009.
C Exact solution of the 3D gravitational 2-body problem: Orbits.
C Dimensional Units: G=1 for gravitational constant.
C Parameters in kepler.par: Fortran Units; nn=2 particles; 
C     nd dimension; G Newton's constant; Time Step dt.
      include 'implicit.sta'
      include 'constants.par'
      include 'kepler.par'
      parameter(norbits=4)
      include 'common_kepler.f'
      dimension r(nd,nn),v(nd,nn),xm(nn)
      ltest=.true.
      ltest=.false.
      write(iuo,'(/,"phy4241orbits: norbits =",I10,/)') norbits
C Set initial values:
      open(iud0,file="Korbits.dat",form='formatted',status='old')
      open(iud1,file="Kpar.d",form='formatted',status='unknown')
      write(iud1,'(/,8X,"L",12X,"Ecm",10X,"p",12X,"e",/)')
      do irpt=1,norbits
        read(iud0,*) i_in,xm(1),(r(id,1),id=1,3),(v(id,1),id=1,3)
        read(iud0,*) i_in,xm(2),(r(id,2),id=1,3),(v(id,2),id=1,3)
        call kepler(xm(1),xm(2),r,v)
        print*,"alpha =",alpha
        write(iuo,'(2X,"irpt, Ecm, xm:",I5,3G15.6)') irpt,Ecm,xm
        write(iuo,'(2X,"    e,   p: ",8X,2G15.6)') e,p
        U=-xm(1)*xm(2)/r12(4)
        write(iuo,'(2X,"  xmu,   U: ",8X,2G15.6)') xmu,U
        xkin=half*xmu*v12(4)**2
        write(iuo,'(2X,"  xkin, |xl|:",7X,2G15.6)') xkin,xl(4)
        write(iuo,'(4X,"r12:",4G15.6)') r12
        write(iuo,'(4X,"v12:",4G15.6)') v12
        write(iuo,'(4X,"r12dv12 =",G25.6)') r12dv12
        write(iuo,'(4X,"xl: ",4G15.6,/)') xl
        write(iud1,'(4F13.6)') xl(4),Ecm,p,e ! Write Kpar.d.
        call korbit(irpt,iud2) ! Creates data on k*.d files for gnuplot.
      enddo
      close(iud0)
      close(iud1)
C
      stop
      end

      include 'vec3d_cross.f'
      include 'vec_dot.f'
      include 'vec_norm.f'
      include 'fe_kin.f'
      include 'kepler.f'
      include 'korbit.f'
