function S = rtsMonteCarloSphere(sqrtN, startAngle, endAngle)

%allocate space for the samples
S = zeros(sqrtN^2, 2);

if sqrtN == 1
    S = [0 0];
    return;
end
cosStart = cos(startAngle);
cosEnd = cos(endAngle);

i=1; % array index
oneoverN = 1.0/sqrtN;
for a = 0:sqrtN-1
    for b = 0:sqrtN-1
        %generate unbiased distribution of spherical coords
        U1 = rand(1, 1);
        U2 = rand(1, 1);
        %x = 1.0 - (a + U1) * oneoverN * (1.0 - cosEnd)
        x = cosStart - (a + U1) * oneoverN * (1.0 - cosEnd);
        y = (b + U2) * oneoverN;
        theta = acos(x);
        phi = 2.0 * pi * y;
        S(i, 1) = theta;
        S(i, 2) = phi;
        i = i+1;
    end
end