function test( N, ubound )

    rs = linspace( 0.0, ubound, N )
    dr = ubound / (N - 1)

    rfunc = rs .* exp( - rs )

    integral = - imag( fft( rfunc ) ) * dr

    file = open( "test.dat", "w" )

    for m = 1 : N

        if m == 1
            q = 0.0
            H = 0.0
        elseif m < div( N, 2 )
            q = (2 * pi / ubound) * (m - 1)
            H = integral[ m ]
        else
            q = - (2 * pi / ubound) * ( N - m + 1 )  # correct?
            H = integral[ m ]
        end

        H_exact = 2 * q / ( q^2 + 1 )^2

        if 0 <= q <= 200.0
            @printf( file, "%20.10e %20.10e %20.10en",
                     q, H, H_exact )
        end
    end

    close( file )
end

test( 10^6, 50.0 )