function Lint = intray(theta, r)
%Usage: Lint = intray(theta, r)
%Compute the line integral of lambda(x,y) over a line L(t) = p + t*np,
%where n = <cos(theta),sin(theta)>, np = <-sin(theta),cos(theta)>, p = r*n.

%Given theta and r set up normal vector n to line and point p on line
%nearest to the origin. Here plus or minus "a" is limits of integration.
n = [cos(theta), sin(theta)];
p = r*n;
a = sqrt(1-r^2);

%Set up a function "phi" defining attentuation along this line
phi = @(t)lambda(p(1)-t*n(2),p(2)+t*n(1));

%And integrate it---but due to nonsmooth phi, need to break up interval or
%quad skips right over compactly supported stuff. Here I'm using N=10
%intervals.
N = 10;
q = -a + 2*a*(0:N)/N;
Lint = 0.0;
for k = 1:N
    Lint = Lint + quad(phi,q(k),q(k+1),1.0e-4);
end
