      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
      character c*2
C n initial conditions, dimension 3, 4 magnitude, 2 particles. 
      parameter(n=4,ngnu=1000,iuo=6,iud0=10,iud1=11,iud2=12)
      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=.false.
      ltest=.true.
C
      write(iuo,'(/,"AngMon: Two particle system.")') 
      write(iuo,'(/,I3," initial conditions.")') n
      write(iuo,'(/,"  Angular momentum xl(), xl(4) magnitude.")')
      pi=acos(-one)
C Initial values from file:
      open(iud0,file="Kepler.txt",form='formatted',status='old')
      open(iud1,file="Orbits0.d",form='formatted',status='unknown')
      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,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
        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)
        write(iuo,'(4X,"r12: ",4F13.6)') r12
        write(iuo,'(4X,"v12: ",4F13.6)') v12
        write(iuo,'(4X,"vcm: ",4F13.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)
        write(iuo,'(4X,"xl,E:",5F13.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
        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
C Prepare plots of effective potentials:
        write(c,'(I2.2)') i
        open(iud2,file="Ueff"//c//".d",
     &       form='formatted',status='unknown')
          do ignu=1,ngnu
            rr=0.01d00*ignu
            Ueff=-xm(1)*xm(2)/rr+(xl(4)/rr)**2/(2*xmu)
            write(iud2,'(3F13.6)') rr,Ueff,E
          enddo
        close(iud2)
C Prepare orbit plots:
        alpha=xm(1)*xm(2)
        e=sqrt(one+2*E*xl(4)**2/alpha**2/xmu) ! Eccentricity.
        p=xl(4)**2/alpha/xmu             ! Half latus rectum.
        theta0=acos((p/r12(4)-one)/e)   ! Initial conditions.
        rdv=r12(1)*v12(1)+r12(2)*v12(2)+r12(3)*v12(3)
        if(rdv.lt.zero) theta0=-theta0
        x0=r12(4)*cos(theta0)
        y0=r12(4)*sin(theta0)
        thetamax=pi
        dtheta=2*thetamax/(2*ngnu)
        if(e.ge.one) then
          thetamax=acos(-one/e)
          dtheta=2*thetamax/(2*ngnu+1)
        endif
        open(iud2,file="Orbit"//c//".d",
     &       form='formatted',status='unknown')
          do ignu=-ngnu,ngnu
            theta=ignu*dtheta
            rr=p/(one+e*cos(theta))
            x=rr*cos(theta)
            y=rr*sin(theta)
            write(iud2,'(4F13.6)') x0,y0,x,y
          enddo
        close(iud2)
C
      enddo
      close(iud0)
      close(iud1)
      if(ltest) stop "AngMon: ltest."
      stop
      end
