Jacobian Subroutine

  subroutine Jacobian ( n, vx, r, dr, x, xr, rx )
  !! Routine that sets up the physical coordinates as well as the Jacobian and
  !! inverse Jacobian to convert derivatives between reference element and
  !! physical element.
  !!
  !! This is almost exactly the same as the routine of the same name in
  !! [[submodule_DG_implementation.f90]] except it also returns xr. Should
  !! probably be combined. I don't think xr is used anywhere so maybe this is
  !! not needed.

    integer(ip), intent(in) :: n
    !! The order of the element.
    real(wp), dimension(2), intent(in) :: vx
    !! The physical coordinates of the left and right boundary of the element.
    real(wp), dimension(n+1), intent(in) :: r
    !! The node locations inside the reference element, \(r_i\).
    real(wp), dimension(n+1,n+1), intent(in) :: dr
    !! The derivative matrix for the reference element, \(D_{ij}\).
    real(wp), dimension(n+1), intent(out) :: x
    !! On output the physical coordinates of the nodes of the element, \(x_i\).
    real(wp), dimension(n+1), intent(out) :: xr
    !! On output contains \((xr)_i=\sum_{i,j}^{n+1} D_{ij} x_j\).
    real(wp), dimension(n+1), intent(out) :: rx
    !! On output contains \((rx)_i=\frac{1}{(xr)_i}\).
    integer(ip) :: i,j

    do i = 1, n+1
      x(i) = vx(1)+0.5_wp*(1.0_wp+r(i))*(vx(2)-vx(1))
    end do
    xr = 0.0_wp
    do i = 1, n+1
      do j = 1, n+1
        xr(i) = xr(i) + dr(i,j) * x(j)
      end do
    end do
    rx = 1.0_wp / xr
    
  end subroutine Jacobian