dslit Subroutine

        subroutine dslit(this, spectrum, dict, seqs)
        !!sample from double slit to produce diff pattern
        ! todo add in user defined slit widths
        ! add correct normalisation
            use random,        only : ranu, ran2, randint, seq
            use sim_state_mod, only : state
            use tomlf,         only : toml_table, get_value
            use constants,     only : TWOPI
            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)

            integer       :: cell(3)
            real(kind=wp) :: x1, y1, z1, x2, y2, z2, a, b, tmp

            call spectrum%p%sample(this%wavelength, tmp)
            this%energy = 1._wp 
            this%fact = TWOPI/(this%wavelength)

            a = 60._wp*this%wavelength !distance between slits
            b = 20._wp*this%wavelength !2 slit width

            if(ran2() > 0.5_wp)then ! pick slit and sample x, y position
                x1 = ranu(a/2._wp,  a/2._wp + b)
                y1 = ranu(-b*0.5_wp, b*0.5_wp)
            else
                x1 = ranu(-a/2._wp,  -a/2._wp - b)
                y1 = ranu(-b*0.5_wp, b*0.5_wp)
            end if

            z2 = 5.0_wp - (1.e-5_wp*(2._wp*(5.0_wp/400._wp)))

            x2 = ranu(-5.0_wp, 5.0_wp)
            y2 = ranu(-5.0_wp, 5.0_wp)
            z1 = (10000._wp * this%wavelength) - 5.0_wp !screen location

            this%pos%x = x2
            this%pos%y = y2
            this%pos%z = z2

            this%phase = sqrt((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)

            this%nxp = (x2 - x1) / this%phase
            this%nyp = (y2 - y1) / this%phase
            this%nzp = -abs((z2 - z1)) / this%phase

            this%tflag = .false.
            this%cnts = 0
            this%bounces = 0
            this%weight = 1.0_wp

            !Set direction cosine/sine
            this%cost = this%nzp
            this%sint = sqrt(1._wp - this%cost**2)

            this%phi = atan2(this%nyp, this%nxp)
            this%cosp = cos(this%phi)
            this%sinp = sin(this%phi)

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

        end subroutine dslit