get_vessels Function

    function get_vessels() result(array)
    !! setup blood vessel scene

        use opticalProperties
        use sdfs,             only : sdf, capsule, box
        use vector_class,     only : vector

        type(sdf), allocatable :: array(:)

        real(kind=wp), allocatable :: nodes(:, :), radii(:)
        integer, allocatable :: edges(:, :)
        integer :: io, edge_cnt, tmp1, tmp2, u, node_cnt, i
        real(kind=wp) :: x, y, z, radius, res, maxx, maxy, maxz
        real(kind=wp) :: musv, muav, gv, nv
        real(kind=wp) :: musd, muad, gd, nd
        type(vector) :: a, b

        type(opticalProp_t) :: opt(2)

        !MCmatlab: an open-source, user-friendly, MATLAB-integrated three-dimensional Monte Carlo light transport solver with heat diffusion and tissue damage
        muav = 231._wp
        musv = 94._wp
        gv = 0.9_wp
        nv = 1.37_wp

        muad = 0.458_wp
        musd = 357._wp
        gd = 0.9_wp
        nd = 1.37_wp

        opt(1) = mono(musv, muav, gv, nv)
        opt(2) = mono(musd, muad, gd, nd)

        !get number of edges
        open(newunit=u, file="res/edges.dat", iostat=io)
        edge_cnt = 0
        do
            read(u,*,iostat=io)tmp1, tmp2
            if(io /= 0)exit
            edge_cnt = edge_cnt + 1
        end do
        close(u)

        !get number of nodes and radii
        open(newunit=u, file="res/nodes.dat", iostat=io)
        node_cnt = 0
        do
            read(u,*,iostat=io)x, y, z
            if(io /= 0)exit
            node_cnt = node_cnt + 1
        end do
        allocate(edges(edge_cnt, 2), nodes(node_cnt, 3), radii(node_cnt))

        !read in edges
        open(newunit=u, file="res/edges.dat", iostat=io)
        do i = 1, edge_cnt
            read(u,*,iostat=io)edges(i, :)
            if(io /= 0)exit
        end do
        close(u)

        !read in nodes
        open(newunit=u, file="res/nodes.dat", iostat=io)
        do i = 1, edge_cnt
            read(u,*,iostat=io)nodes(i, :)
            if(io /= 0)exit
        end do
        close(u)

        !read in radii
        open(newunit=u, file="res/radii.dat", iostat=io)
        do i = 1, node_cnt
            read(u,*,iostat=io)radii(i)
            if(io /= 0)exit
        end do
        close(u)

        res = 0.001_wp!0.01mm
        maxx = maxval(abs(nodes(:, 1)))
        maxy = maxval(abs(nodes(:, 2)))
        maxz = maxval(abs(nodes(:, 3)))

        nodes(:, 1) = (nodes(:, 1) / maxx) - 0.5_wp
        nodes(:, 2) = (nodes(:, 2) / maxy) - 0.5_wp
        nodes(:, 3) = (nodes(:, 3) / maxz) - 0.5_wp
        nodes(:, 1) = nodes(:, 1) * maxx * res
        nodes(:, 2) = nodes(:, 2) * maxy * res
        nodes(:, 3) = nodes(:, 3) * maxz * res

        allocate(array(edge_cnt+1))

        do i = 1, edge_cnt
            a = vector(nodes(edges(i, 1), 1), nodes(edges(i, 1), 2), nodes(edges(i, 1), 3))
            b = vector(nodes(edges(i, 2), 1), nodes(edges(i, 2), 2), nodes(edges(i, 2), 3))
            radius = radii(edges(i, 1)) * res
            array(i) = capsule(a, b, radius, opt(1), 1)
        end do

        array(i) = box(vector(.32_wp, .18_wp, .26_wp), opt(2), 2)

    end function get_vessels