function imagecoords = ProjectToRetina(eyes,points)
% Calculate retinal coordinates of a point in space.
% "eyes" is a structure giving the binocular eye position.
% "points" is a 3-row vector, [X;Y;Z]. Each column is assumed to specify a
% different point in space.
% NB this is exact; no approximations of any kind are used.
global L R CYC DELTA HALFDELTA 

[nrow,npts]=size(points);
if nrow~=3
    error('error in ProjectToRetina: points must have 3 rows')
end    

for thiseye=[L R]
    n = eyes(thiseye).nodal;
    Rot = eyes(thiseye).M;
    for jpt=1:npts
        % Go through each point and project it onto the retina
        point = points(:,jpt);
        nrm = dot(point-n , Rot*[0;0;1]);
        xonf = - dot(point-n , Rot*[1;0;0]) / nrm;
        yonf = - dot(point-n , Rot*[0;1;0]) / nrm;

        imagecoords.AzLong(thiseye,jpt) = atan(xonf);
        imagecoords.AzLat(thiseye,jpt) = atan(xonf/sqrt(yonf^2+1));
        imagecoords.ElLong(thiseye,jpt) = atan(yonf);
        imagecoords.ElLat(thiseye,jpt) = atan(yonf/sqrt(xonf^2+1));
    end
end

% Work out retinal disparity and cyclopean location of each point:
flds=fieldnames(imagecoords);
for j=1:length(flds)
    imagecoords.(flds{j})(DELTA,:) = imagecoords.(flds{j})(R,:) - imagecoords.(flds{j})(L,:);
    imagecoords.(flds{j})(CYC,:) = (imagecoords.(flds{j})(R,:) + imagecoords.(flds{j})(L,:))/2;
end
% NB cyclopean eccentricity is the mean eccentricity in left and right eyes
% Work out retinal eccentricity of each point:
imagecoords.Eccentricity = acos(cos(imagecoords.AzLat).*cos(imagecoords.ElLong));

