dpi1_0.for

!******************************************************************************
!
!                                 DPI Library
!                       Dynamical Plug-In for Mercury 6
!                                Release 1.0
!
!                    Copyright (C) 2002-2015  Diego Turrini
!
!     This program is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
! 
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
! 
!     You should have received a copy of the GNU General Public License
!     along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!------------------------------------------------------------------------------
!
!     DPI Library
!     Dynamical Plug-In for Mercury 6
!
!     Release: 1.0
!     Author: Diego Turrini
!     E-mail: diego.turrini_at_iaps.inaf.it
!     Last modified: January 2010
!     Modified by: Diego Turrini
!
!     Disclaimer: this library is supplied as a plug-in for the Mercury software
!                 developed by John E. Chamber: specifically, it is designed to 
!                 work with version 6 of Mercury. To facilitate its use in 
!                 conjunction with Mercury, the design of the subroutines was
!                 kept as similar as possible to that of analogous subroutines
!                 in Mercury. Mercury is developed by John E. Chambers and can
!                 be downloaded at the following URL:
!                 http://www.arm.ac.uk/~jec/home.html (see also the 
!                 Astrophysics Source Code Library, ascl id. 1201.008) 
!                 The DPI library contains some subroutines derived and adapted 
!                 from Mercury 6 and originally developed by John E. Chambers. 
!                 Information on the author and, where appropriate, the name of
!                 the original subroutines are reported in the headers.  
!
!     Acknowledgements: the author wishes to thank Patricia Eleanor Verrier
!                       for her assistance in debugging and restructuring the
!                       algorithm and the code through comparison with her
!                       MOIRAI software.
!
!     Bibliographic references:
!
!     Mercury 6 software:
!
!     Chambers J. E., 1999, Monthly Notices of the Royal Astronomical Society, 
!     vol. 304, pp. 793-799 
!
!     Symplectic mapping for S-type binary star systems:
!
!     Chambers J. E., Quintana E. V., Duncan M. J., Lissauer J. J., 2002, The
!     Astrophysical Journal, vol. 123, pp. 2884-2894
!
!     MOIRAI software:
! 
!     Verrier P. E., Evans N. W., 2007, Monthly Notices of the Royal 
!     Astronomical Society, vol. 382, pp. 1432-1446
!
!     DPI library:
!
!     Turrini D., Barbieri M., Marzari F., Thebault P., Tricarico P., 2005, 
!     Memorie della Societa' Astronomica Italiana Supplementi, vol. 6,
!     pp. 172-177
!
!     Turrini D., Barbieri M., Marzari F., Tricarico P., 2004, Memorie della
!     Societa' Astronomica Italiana Supplementi, vol. 5, pp. 127-130
!
!     Thebault P., Marzari F., Scholl H., Turrini D., Barbieri M., 2004, 
!     Astronomy & Astrophysics, vol. 427, pp. 1097-1104
!
!******************************************************************************
!
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to convert from the initial, inertial coordinates to
!     S-type binary coordinates as in Chambers et al. (2002)
!     N.B.: the subroutine computes positions and pseudovelocities, not the
!           state vectors
!     N.B.: the initial coordinates are centered on the initial position of the
!           central star yet are not heliocentric
!     Author: Diego Turrini
!     Last modified: May 2009
!******************************************************************************
      subroutine dpi_h2wb(nbig,nbod,m,xh,vh,xb,pvb)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 m(nbod)
      real*8 xh(3,nbod),vh(3,nbod)
      real*8 xb(3,nbod),pvb(3,nbod)
!     Local variables
      integer i,j
      real*8 mtot,mden,mvh(3),mxh(3)
!     Variables initialization
      do i=1,3
        mxh(i)=0.d0
        mvh(i)=0.d0
      end do
      do j=1,nbod
        do i=1,3
          xb(i,j)=0.d0
          pvb(i,j)=0.d0
        end do
      end do
!     Computing support variables
      mtot=m(1)
      if (nbig.gt.2) then
        do j=2,nbig-1
          mtot=mtot+m(j)
          do i=1,3
            mvh(i)=mvh(i)+m(j)*vh(i,j)
            mxh(i)=mxh(i)+m(j)*xh(i,j)
          end do
        end do
      end if
      mden=mtot
      mtot=mtot+m(nbig)
!     Computing transformed positions for S type binary systems
!     Computing position of central star
      do i=1,3
        xb(i,1)=(m(1)*xh(i,1)+m(nbig)*xh(i,nbig)+mxh(i))/mtot
      end do
!     If any, computing positions of massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            xb(i,j)=xh(i,j)-xh(i,1)
          end do
        end do
      end if
!     Computing position of binary star
      do i=1,3
        xb(i,nbig)=xh(i,nbig)-(m(1)*xh(i,1)+mxh(i))/mden
      end do
!     If any, computing positions of massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            xb(i,j)=xh(i,j)-xh(i,1)
          end do
        end do
      end if
!     Computing transformed pseudo-velocities for S type binary systems
!     Computing pseudo-velocity of central star
      do i=1,3
        pvb(i,1)=(m(1)*vh(i,1)+m(nbig)*vh(i,nbig)+mvh(i))/mtot
      end do
!     If any, computing pseudo-velocities of massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            pvb(i,j)=vh(i,j)-(m(1)*vh(i,1)+mvh(i))/mden
          end do
        end do
      end if
!     Computing pseudo-velocity of binary star
      do i=1,3
        pvb(i,nbig)=vh(i,nbig)-pvb(i,1)
      end do
!     If any, computing pseudo-velocities of massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            pvb(i,j)=vh(i,j)-(m(1)*vh(i,1)+mvh(i))/mden
          end do
        end do
      end if
      end subroutine dpi_h2wb
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to convert from the S-type binary coordinates to the initial,
!     inertial coordinates as in Chambers et al. (2002)
!     N.B.: the subroutine needs positions and pseudovelocities, not the
!           state vectors
!     N.B.: the initial coordinates are centered on the initial position of the
!           central star yet are not heliocentric
!     Author: Diego Turrini
!     Last modified: January 2010
!******************************************************************************
      subroutine dpi_wb2h(nbig,nbod,m,xh,vh,xb,pvb)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 m(nbod)
      real*8 xh(3,nbod),vh(3,nbod)
      real*8 xb(3,nbod),pvb(3,nbod)
!     Local variables
      integer i,j
      real*8 mtot,mden,mvb(3),mxb(3)
!     Variables initialization