      program AngMom ! Bernd Berg, Feb 2 2011.
C Angular Momentum in CM from fro given initial conditions.
C LL: Landau-Lifshitz Mechanics.
      implicit real*8 (a-h,o-z)
      logical ltest
C n initial conditions, dimension 3, magnitude in 4, 2 particles.
      parameter(n=4,iuo=6,iud0=10)
      dimension r(3,2),v(3,2),xm(2)
      dimension r12(4),v12(4),rcm(4),vcm(4),xl(4)
c     data r/8*0.d00/
      ltest=.true.
      ltest=.false.
      if(ltest) stop "AngMon: ltest."
C
      write(iuo,'(/,"AngMon: Two particle system.")')
      write(iuo,'(/,I3," initial conditions.")') n
      write(iuo,'(/,"  Angular momentum xl(), xl(4) magnitude.")')
C Initial values from file:
      open(iud0,file="Kepler.txt",form='formatted',status='old')
      do i=1,n
        write(iuo,'(/,"Set",I2,":",/)') i
        read(iud0,*) i_in,xm(1),(r(id,1),id=1,3),(v(id,1),id=1,3)
        write(iuo,'(7F9.4)') 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)
        write(iuo,'(7F9.4)') xm(2),(r(id,2),id=1,3),(v(id,2),id=1,3)
        XMM=xm(1)+xm(2)     ! Total mass (mu in LL).
        xmu=xm(1)*xm(2)/XMM ! Reduced mass (m in LL).
        write(iuo,'(/,4X,"XMM,xmu:",10X,2F14.6)') XMM,xmu
        r12(4)=0
        v12(4)=0
        vcm(4)=0
        do id=1,3
          r12(id)=r(id,1)-r(id,2)
          r12(4)=r12(4)+r12(id)**2
          v12(id)=v(id,1)-v(id,2)
          v12(4)=v12(4)+v12(id)**2
          vcm(id)=(xm(1)*v(id,1)+xm(2)*v(id,2))/XMM
          vcm(4)=vcm(4)+vcm(id)**2
        enddo
        r12(4)=sqrt(r12(4))
        vcm(4)=sqrt(vcm(4))
        write(iuo,'(4X,"r12:",4F14.6)') r12
        write(iuo,'(4X,"v12:",4F14.6)') v12
        write(iuo,'(4X,"vcm:",4F14.6)') vcm
        xl(1)=xmu*(r12(2)*v12(3)-r12(3)*v12(2))
        xl(2)=xmu*(r12(3)*v12(1)-r12(1)*v12(3))
        xl(3)=xmu*(r12(1)*v12(2)-r12(2)*v12(1))
        xl(4)=sqrt(xl(1)**2+xl(2)**2+xl(3)**2)
        E=xmu*v12(4)**2/2+xm(1)*xm(2)/r12(4) ! Energy
        write(iuo,'(4X,"xl,E:",F13.6,4F14.6)') xl,E
        write(iuo,'(4X,"Second Calculation (contol check):")')
        do id=1,3 ! CM frame coordinates.
          r(id,1)=r(id,1)-rcm(id)
          r(id,2)=r(id,2)-rcm(id)
          v(id,1)=v(id,1)-vcm(id)
          v(id,2)=v(id,2)-vcm(id)
        enddo
        print*,v
        E=xmu*v12(4)**2/2+xm(1)*xm(2)/r12(4) ! Energy
        xl(1)=xm(1)*(r(2,1)*v(3,1)-r(3,1)*v(2,1))
     &       +xm(2)*(r(2,2)*v(3,2)-r(3,2)*v(2,2))
        xl(2)=xm(1)*(r(3,1)*v(1,1)-r(1,1)*v(3,1))
     &       +xm(2)*(r(3,2)*v(1,2)-r(1,2)*v(3,2))
        xl(3)=xm(1)*(r(1,1)*v(2,1)-r(2,1)*v(1,1))
     &       +xm(2)*(r(1,2)*v(2,2)-r(2,2)*v(1,2))
        print*,xl
        write(iuo,'(4X,"xl,E:",F13.6,4F14.6)') xl,E
      enddo
      close(iud0)
C
      stop
      end