//projection of rays onto plane
Di=-1; //[L] intersection of image plane with Z axis
ptx=-1.5:1.5; //[L] object point x-coordinates before rotation
pty=ptx; //[L] object point y-coordinates before rotation
gridx=repmat(ptx,1,length(pty));
gridy=repmat(pty,length(ptx),1);
gridy=gridy(:)';
ptz=10:10:100; //[L] object point z-coordinates before rotation

ry=0/180*%pi; //[rad] rotation of point about Y axis
rz=0/180*%pi; //[rad] further rotation of point about Z axis

c=cos(ry);
s=sin(ry);
my=[c,0,s;0,1,0;-s,0,c]; //rotation matrix about Y axis
c=cos(rz);
s=sin(rz);
mz=[c,-s,0;s,c,0;0,0,1]; //rotation matrix about Z axis

scp=%f; //plot 3D scatter figure?

clf();

numz=length(ptz);
cmap=rainbowcolormap(numz);
gcf().color_map = cmap;
colorbar(ptz(1),ptz($));

for jz=numz:-1:1 //1:numz
    gridz=ones(1,length(gridy))*ptz(jz);
    pt=[gridx;gridy;gridz]; //[L] 3D points location
    
    pt=my*pt;
    pt=mz*pt;
    pt(4,:)=sqrt(sum(pt.^2,1)); //store distance information for later
    
    if(scp)
        scatter3(pt(1,:),pt(2,:),pt(3,:),'markerEdgeColor',cmap(jz,:)); //plot 3D position
    else
        q=Di./pt(3,:); //project onto image plane
        pt(1,:)=pt(1,:).*q;
        pt(2,:)=pt(2,:).*q;
        //pt(3,:)=pt(3,:).*q; //unnecessary as answer is always Di
        plot(pt(1,:),pt(2,:),'o','MarkerEdgeColor',cmap(jz,:)); //plot projected points
    end
end

isoview on
