scatter Subroutine

        subroutine scatter(this, hgg, g2, dects)
            !! Scattering routine. Implments both isotropic and henyey-greenstein scattering
            !! taken from [mcxyz](https://omlc.org/software/mc/mcxyz/index.html)

            use constants, only : PI, TWOPI, wp
            use random,    only : ran2
            use detectors, only : dect_array

            class(photon), intent(inout) :: this
            !> g factor
            real(kind=wp), intent(in)    :: hgg
            !> g factor squared
            real(kind=wp), intent(in)    :: g2
            !> array of detectors. Only used if biased scattering is enabled.
            type(dect_array), optional, intent(in) :: dects(:)

            real(kind=wp) :: temp, uxx, uyy, uzz, a, p

            a = 0.9_wp
            p = 0.0_wp

            if(hgg == 0.0_wp)then
                !isotropic scattering
                this%cost = 2._wp * ran2() - 1._wp
            else
                !henyey-greenstein scattering
                if(ran2() < p .and. present(dects))then
                    !bias scattering
                    temp = ran2()*((1._wp / (1._wp - a)) - (1._wp / sqrt(a**2 + 1._wp))) + (1._wp/sqrt(a**2 + 1._wp))
                    temp = temp**(-2._wp)
                    this%cost = (1._wp/(2._wp*a)) * (a**2 +1._wp - temp)
                    this%nxp = dects(1)%p%pos%x - this%pos%x
                    this%nyp = dects(1)%p%pos%y - this%pos%y
                    this%nzp = dects(1)%p%pos%z - this%pos%z
                else
                    !unbiased
                    temp = (1.0_wp - g2) / (1.0_wp - hgg + 2._wp*hgg*ran2())
                    this%cost = (1.0_wp + g2 - temp**2) / (2._wp*hgg)
                end if
            end if

            this%sint = sqrt(1._wp - this%cost**2)

            this%phi = TWOPI * ran2()
            this%cosp = cos(this%phi)
            if(this%phi < PI)then
                this%sinp = sqrt(1._wp - this%cosp**2)
            else
                this%sinp = -sqrt(1._wp - this%cosp**2)
            end if

            if(1._wp - abs(this%nzp) <= 1e-12_wp)then ! near perpindicular
                uxx = this%sint * this%cosp
                uyy = this%sint * this%sinp
                uzz = sign(this%cost, this%nzp)
            else
                temp = sqrt(1._wp - this%nzp**2)
                uxx = this%sint * (this%nxp * this%nzp * this%cosp - this%nyp * this%sinp) / temp + this%nxp * this%cost
                uyy = this%sint * (this%nyp * this%nzp * this%cosp + this%nxp * this%sinp) / temp + this%nyp * this%cost
                uzz = -1._wp*this%sint * this%cosp * temp + this%nzp * this%cost
            end if

            this%nxp = uxx
            this%nyp = uyy
            this%nzp = uzz

        end subroutine scatter