boundary_edge Subroutine

subroutine boundary_edge(bplane, bedge, nb)
    !+ Returns 3 x `nb` array containing points along the BoundedPlane's boundary edge
    type(BoundedPlane), intent(in)             :: bplane
        !+ Bounded plane
    real(Float64), dimension(:,:), intent(out) :: bedge
        !+ Boundary edge points of bounded plane
    integer, intent(out)                       :: nb
        !+ Number of points in boundary edge

    integer :: i
    real(Float64) :: th, dth, x, y
    real(Float64), dimension(4) :: xx, yy

    select case (bplane%shape)
        case (1) !Rectangular boundary
            nb = 4
            if(nb.gt.size(bedge,2)) then
                if(inputs%verbose.ge.0) then
                    write(*,'("BOUNDARY_EDGE: Incompatible boundary edge array : ",i2," > ",i2)') nb, size(bedge,2)
                endif
                stop
            endif
            xx = [-bplane%hw,-bplane%hw,bplane%hw,bplane%hw]
            yy = [-bplane%hh,bplane%hh,bplane%hh,-bplane%hh]
            do i=1,nb
                bedge(:,i) = matmul(bplane%basis,[xx(i),yy(i),0.d0]) + bplane%origin
            enddo
        case (2)
            nb = 50
            if(nb.gt.size(bedge,2)) then
                if(inputs%verbose.ge.0) then
                    write(*,'("BOUNDARY_EDGE: Incompatible boundary edge array : ",i2," > ",i2)') nb, size(bedge,2)
                endif
                stop
            endif
            dth = 2*pi/nb
            do i=1,nb
                th = i*dth
                x = bplane%hw*cos(th)
                y = bplane%hh*sin(th)
                bedge(:,i) = matmul(bplane%basis,[x,y,0.d0]) + bplane%origin
            enddo
        case default
            if(inputs%verbose.ge.0) then
                write(*,'("BOUNDARY_EDGE: Unknown boundary shape: ",i2)') bplane%shape
            endif
            stop
    end select

end subroutine boundary_edge