      program AngMom ! Bernd Berg, Oct 02 2011.
C Angular Momentum in CM from for given initial conditions.
C This program will be part of the complete program Kepler.f.
C LL: Landau-Lifshitz Mechanics.
      implicit real*8 (a-h,o-z)
      logical ltest
      character c*2
C n initial conditions, dimension 3, 4 magnitude, 2 particles. 
      parameter(n=4,ngnu=1000,iuo=6,iud0=10,iud1=11)
      parameter(zero=1.d00,one=1.d00)
      dimension r(4,2),v(4,2),xm(2), r12(4),v12(4),rcm(4),vcm(4),xl(4)
      ltest=.true.
      ltest=.false.
      if(ltest) stop "AngMon: ltest."
C
      write(iuo,'(/,"AngMon:",/,/,
     &  "  Angular momentum of 2-particles in their rest frame.")')
      write(iuo,'(/,I3," initial conditions.")') n
      write(iuo,'(/,"  Angular momentum xl(), xl(4) magnitude.")')
      write(iuo,'("  XMM total mass (xmu in Landau-Lifshitz).")')
      write(iuo,'("  Xmu reduced mass (m in Landau-Lifshitz).")')
      pi=acos(-one)
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,'(4X,"Mass, initial coordinates and velocities:")')
        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,2F13.6)') XMM,xmu
        do id=1,3
          r12(id)=r(id,1)-r(id,2)
          v12(id)=v(id,1)-v(id,2)
          vcm(id)=(xm(1)*v(id,1)+xm(2)*v(id,2))/XMM
          rcm(id)=(xm(1)*r(id,1)+xm(2)*r(id,2))/XMM ! at time zero.
        enddo
        r12(4)=sqrt(r12(1)**2+r12(2)**2+r12(3)**2)
        v12(4)=sqrt(v12(1)**2+v12(2)**2+v12(3)**2)
        vcm(4)=sqrt(vcm(1)**2+vcm(2)**2+vcm(3)**2)
        rcm(4)=sqrt(rcm(1)**2+rcm(2)**2+rcm(3)**2)
        write(iuo,'(4X,"r12: ",4F13.6)') r12
        write(iuo,'(4X,"v12: ",4F13.6)') v12
        write(iuo,'(4X,"vcm: ",4F13.6)') vcm
        write(iuo,'(4X,"rcm: ",4F13.6,3X,"at time zero.")') rcm
        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)
        write(iuo,'(4X,"xl,E:",5F13.6)') xl,E
        write(iuo,'(4X,"Second Calculation (control check):")')
C Note: The transformation of the coordinates does not effect
C the angular momentums, because of the parallel axis theorem
C and the fact that \vec{P} is zero in any rest frame. The 
C transformation of the velocities is essential.
        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
        xl(1)=xm(1)*(r(2,1)*v(3,1)-r(3,1)*v(2,1))
        r(4,1)=sqrt(r(1,1)**2+r(2,1)**2+r(3,1)**2)
        r(4,2)=sqrt(r(1,2)**2+r(2,2)**2+r(3,2)**2)
        v(4,1)=sqrt(v(1,1)**2+v(2,1)**2+v(3,1)**2)
        v(4,2)=sqrt(v(1,2)**2+v(2,2)**2+v(3,2)**2)
        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))
        E=xm(1)*v(4,1)**2/2+xm(2)*v(4,2)**2/2-xm(1)*xm(2)/r12(4)
        write(iuo,'(4X,"xl,E:",5F13.6)') xl,E
      enddo
      close(iud0)
      stop
      end
