add_geiger_ping Subroutine

    subroutine add_geiger_ping(tape, track, t1, t2, f, Amp)
        type(tape_recorder), intent(inout) :: tape
        integer, intent(in)  :: track
        real(wp), intent(in) :: t1, t2, f, Amp
        real(wp) :: b

        ! Pulsation (radians/second):
        real(wp) :: omega
        ! Time in seconds:
        real(wp) :: t
        integer  :: i

        omega = 2 * PI * f
        t = 0._wp
        do i = nint(t1*RATE), nint(t2*RATE)-1
            ! Bessel functions of the first kind: a short ping
            b = Amp * bessel_jn(1, omega*t) * bessel_jn(2, omega*t)
            tape%left( track, i) = tape%left( track, i) + b
            tape%right(track, i) = tape%right(track, i) + b

            t = t + dt
        end do
    end subroutine add_geiger_ping