pascal_2D_quad Module Function

pure module function pascal_2D_quad(N, x, y) result(row)
!*
! Generates an array of points related to a quadrilateral using Pascal's
! triangle in 2D, where rows are 0-indexed
!
! Pascal's triangle in 2D looks like this, with points used in bi-quadratic quadrilateral in bold:
!   \[ [\mathbf{1}] \]
!   \[ [\mathbf{x},~ \mathbf{y}] \]
!   \[ [\mathbf{x^2},~ \mathbf{x y},~ \mathbf{y^2}] \]
!   \[ [x^3,~ \mathbf{x^2y},~ \mathbf{xy^2},~ y^3] \]
!   \[ [x^4,~ x^3y,~ \mathbf{x^2y^2},~ xy^3, y^4] \]
!   \[ \vdots \]
!   \[ [x^N,~ x^{N-1}y,~ \cdots ~,~ xy^{N-1},~ y^N] \]

integer,  intent(in)  :: N            !! Order of the qaudrilateral
real(wp), intent(in)  :: x            !! X-coordinate of node used in calculation
real(wp), intent(in)  :: y            !! Y-coordinate of node used in calculation
real(wp), dimension((N+1)**2) :: row  !! Output row

integer :: ii
real(wp), dimension(:), allocatable :: temp_pre, temp_post

row = 0.d0

! Collects the first N rows of a 2D pascal triangle as function of x and y
temp_pre = [( pascal_row_2D(ii, x, y), ii = 0, N )]
temp_post = [( pascal_2D_quad_post(N, ii, x, y), ii = N+1, 2*N )]

row = [temp_pre, temp_post]

return
  contains
      pure function pascal_2D_quad_post(N_, ii_, x_, y_) result(row_)
          integer,  intent(in)  :: N_, ii_
          real(wp), intent(in)  :: x_, y_
          real(wp), dimension(2*N_-ii_+1) :: row_

          integer :: start, finish
          real(wp), dimension(:), allocatable :: temp

          temp = pascal_row_2D(ii_, x_, y_)
          start = ii_-N_+1
          finish = N_+1

          row_ = temp( start:finish )

          return
      end function pascal_2D_quad_post
  end function pascal_2D_quad