coord_trans Subroutine

  subroutine coord_trans ( a, b, xp, xip, dxpdt, d2xpdt2, xi, &
                           x, dxdt, dxdxi, d2xdt2, d2xdxi2, d2xdtdxi )
  !! Routine that calculates the Tortoise coordinates \((t,r_*)\) from time
  !! dependent coordinates \((\lambda,\xi)\) where the particle is kept at a
  !! fixed coordinate location as well as the informtation needed to transform
  !! the wave equation to time dependent coordinates.

    implicit none

    real(wp), intent(in) :: a
    !! The lower boundary location of the time dependent coordinate region.
    real(wp), intent(in) :: b
    !! The upper boundary location of the time dependent coordinate region.
    real(wp), intent(in) :: xp
    !! The current particle location in Tortoise coordinates,
    !! \(r^{\mathrm{p}}_*\).
    real(wp), intent(in) :: xip
    !! The constant particle location in time dependent coordinates,
    !! \(\xi^{\mathrm{p}}\).
    real(wp), intent(in) :: dxpdt
    !! The current time derivative of the particle location in Tortoise
    !! coordinates, \(\frac{d r_*^{\mathrm{p}}}{dt}\).
    real(wp), intent(in) :: d2xpdt2
    !! The current second time derivative of the particle location in Tortoise
    !! coordinates, \(\frac{d^2 r_*^{\mathrm{p}}}{dt^2}\).
    real(wp), dimension(:), intent(in) :: xi
    !! A 1d array containing coordinate values \(\xi_i\).
    real(wp), dimension(:), intent(out) :: x
    !! A 1d array that on output contains the Tortoise coordinates \((r_*)_i\).
    real(wp), dimension(:), intent(out) :: dxdt
    !! A 1d array that on outout contains \(\frac{d r_*}{dt}\).
    real(wp), dimension(:), intent(out) :: dxdxi
    !! A 1d array that on outout contains \(\frac{d r_*}{d\xi}\).
    real(wp), dimension(:), intent(out) :: d2xdt2
    !! A 1d array that on outout contains \(\frac{d^2 r_*}{dt^2}\).
    real(wp), dimension(:), intent(out) :: d2xdxi2
    !! A 1d array that on outout contains \(\frac{d^2 r_*}{d\xi^2}\).
    real(wp), dimension(:), intent(out) :: d2xdtdxi
    !! A 1d array that on outout contains \(\frac{d^2 r_*}{dt d\xi}\).
    integer(ip) :: i
    real(wp) :: xpma, xipma, xima, bmxp, bmxip, bmxi, bma, ximxip, xipmxp
    real(wp) :: xipmainv, xipmamulbmxipinv, dtfac

    do i = 1, size(xi)
      xpma = xp - a
      xipma = xip - a
      xima = xi(i) - a
      bmxp = b - xp
      bmxi = b - xi(i)
      bmxip = b - xip
      bma = b - a
      ximxip = xi(i) - xip
      xipmainv = 1.0_wp/xipma
      xipmamulbmxipinv = xipmainv/bmxip
      xipmxp = xip - xp
      x(i) = a + xpma*xipmainv*xima &
           +( bmxp*xipma-xpma*bmxip )*xipmamulbmxipinv/bma*xima*ximxip
      dtfac = xima*bmxi*xipmamulbmxipinv
      dxdt(i) = dtfac*dxpdt
      dxdxi(i) = ( (2.0_wp*xi(i)-xip-a)*xipmxp + xpma*bmxip ) &
                 * xipmamulbmxipinv
      d2xdt2(i) = dtfac*d2xpdt2
      d2xdxi2(i) = 2.0_wp*xipmxp*xipmamulbmxipinv
      d2xdtdxi(i) = ( a + b - 2.0_wp * xi(i) ) * xipmamulbmxipinv * dxpdt
    end do
  end subroutine coord_trans