function [X, Y, Z] = quadraytrace(M, Xlim, Ylim, Zlim, pnum)
%QUADRAYTRACE generates three vectors X, Y, Z suitable for 3-D plotting
% using scatter3. It can plot ellipsoids, paraboloids and hyperboloids.
% M is the function representing the ellipsoids, paraboloids or
% hyperboloids and has to be invertible.
% Xlim, Ylim, Zlim are vectors with two element that specify the
% upper and lower limits of the corresponding axes and pnum is
% the number of points (size of vectors X, Y, Z) returned.
%   Minas Spetsakis, Nov. 2024

X = zeros(1, pnum);
Y = zeros(1, pnum);
Z = zeros(1, pnum);

i = 1;
while i<=pnum
    p = [Xlim(1) + rand*(Xlim(2)-Xlim(1))
        Ylim(1) + rand*(Ylim(2)-Ylim(1))
        Zlim(1) + rand*(Zlim(2)-Zlim(1))];
    d = randn(3,1);
    Mi = inv(M);
    dMd = d'*Mi*d;
    pMp = p'*Mi*p;
    pMd = p'*Mi*d;
    discrim = sqrt(pMd^2-dMd*(pMp-1));
    if ~isreal(discrim)
        continue;
    end
    lam1 = (-pMd+discrim)/dMd;
    lam2 = (-pMd-discrim)/dMd;
    p1 = p+lam1*d;
    if isinside(p1,Xlim, Ylim, Zlim)
        X(i) = p1(1);
        Y(i) = p1(2);
        Z(i) = p1(3);
        i = i+1;
        if i>pnum
            break;
        end
    end
    p2 = p+lam2*d;
    if isinside(p2,Xlim, Ylim, Zlim)
        X(i) = p2(1);
        Y(i) = p2(2);
        Z(i) = p2(3);
        i = i+1;
        if i>pnum
            break;
        end
    end
end

end

function [res] = isinside(p, Xlim, Ylim, Zlim)
%ISINSIDE returns true if p is within the X, Y, Z limits
% of Xlim, Ylim, Zlim. These are 2-element vectors.

res = 0;
if p(1)<Xlim(1)
    return
end
if p(1)>Xlim(2)
    return
end
if p(2)<Ylim(1)
    return
end
if p(2)>Ylim(2)
    return
end
if p(3)<Zlim(1)
    return
end
if p(3)>Zlim(2)
    return
end
res = 1;
return;
end
