annulus Subroutine

        subroutine annulus(this, spectrum, dict, seqs)
            !! annular source
            use constants,     only : TWOPI 
            use utils,         only : deg2rad
            use tomlf,         only : toml_table, get_value
            use random,        only : ran2, rang, seq
            use sim_state_mod, only : state
            use piecewiseMod

            class(photon) :: this
            type(spectrum_t), intent(in) :: spectrum
            type(toml_table), optional, intent(inout) :: dict
            type(seq), optional, intent(inout) :: seqs(2)

            character(len=:), allocatable :: beam_type
            real(kind=wp) :: beta, rlo, rhi, radius, tmp, mid, angle, x, y, z, phi, sinp, cosp
            type(vector)  :: pos
            integer       :: cell(3)

            call get_value(dict, "beta", beta)
            call get_value(dict, "radius", rlo)
            call get_value(dict, "radius_hi", rhi)
            call get_value(dict, "annulus_type", beam_type)

            if(beam_type == "tophat")then
                radius = rlo + (rhi - rlo) * sqrt(ran2())
            elseif(beam_type == "gaussian")then
                mid = (rhi - rlo) / 2.
                call rang(radius, tmp, mid, 0.04_wp)
            else
                error stop "No such beam type!"
            end if

            phi = TWOPI * ran2()

            angle = deg2rad(beta)

            cosp = cos(phi)
            sinp = sin(phi)
            x = radius * cosp
            y = radius * sinp
            z = state%grid%zmax - 1e-8_wp! just inside surface of medium. TODO make this user configurable?
            pos = vector(x, y, z)
            this%pos = pos

            this%nxp = sin(angle) * cosp
            this%nyp = sin(angle) * sinp
            this%nzp = -cos(angle)

            this%phi = phi
            this%cosp = cosp
            this%sinp = sinp
            this%cost = this%nzp
            this%sint = sqrt(1._wp - this%cost**2)

            this%tflag = .false.
            this%bounces = 0
            this%cnts = 0
            this%weight = 1.0_wp
            call spectrum%p%sample(this%wavelength, tmp)

            ! Linear Grid 
            cell = state%grid%get_voxel(this%pos)
            this%xcell = cell(1)
            this%ycell = cell(2)
            this%zcell = cell(3)
        end subroutine annulus