function Fig_VDispDefinition
% Figure to illustrate projecting the point P onto the two retina:
% show two different choices of retinal coordinates

clear all
close all

fs=14; % fontsize to use for axis labels
fsa=18; % fontsize to use for panel labels
fst = 18; % fontsize to use for titles

% Head-centered coordinates of viewed object
objectP = [-6 7 10];

% Coordinates of real eye:
interoculardistance = 6;
eyes.radius = 1.75;
eyes.eyeballcolor = [0.75 0.75 0.5];
eyes.eyeballtransparency = 1;
% Specify eye positions
eyes.eyepos(1).nodal = [1; 0 ; 0]*interoculardistance/2 ;
eyes.eyepos(1).azimuth = -15;
eyes.eyepos(1).elevation = 0;
eyes.eyepos(1).torsion = 0;
% Colors in which to draw coordinate systems
eyes.azcoordcol = 'k' ;
eyes.elcoordcol = 'k' ;
eyes.eyepos(2) = eyes.eyepos(1);
eyes.eyepos(2).azimuth = -eyes.eyepos(1).azimuth;
eyes.eyepos(2).nodal = [-1; 0 ; 0]*interoculardistance/2;

% Both retinae superimposed
jointeyes = eyes;
jointeyes = rmfield(jointeyes,'eyepos');
jointeyes.eyepos.nodal = [ 0 ; 0 ; 0];
jointeyes.eyepos.azimuth = 0;
jointeyes.eyepos.elevation = 0;
jointeyes.eyepos.torsion = 0;

% Specify the coordinate system to be painted on the retina
coordsystem.nlines = 11;
coordsystem.spacing = 15;
coordsystem.sca  = 1.01; % amount by which to exaggerate eyeball radius so coord lines are visible
coordsystem.linewidth = 0.5;
coordsystem.linewidththick = 1.5; % used for plotting meridians

figure('pos',[ 479   161   686   529])
for jcoordsystem = 1:2
    if jcoordsystem == 1
        coordsystem.name = 'azlongellong';
    elseif jcoordsystem == 2
        coordsystem.name = 'azlongellat';
    end
    ax(jcoordsystem)=subplot(2,2,jcoordsystem);

    arrowheadlength = 0.5;
    % Draw interocular axis
    Xaxis = [ [-1 1]*6;[0 0];[0 0] ] + [1 1 ; 0 0 ; 0 0];
    plotvector(Xaxis,'k-','linew',2)
    hold on
    % Add on arrowhead
    plot3(Xaxis(1,2) + [-1 0 -1]*arrowheadlength , Xaxis(2,2) * [1 1 1], [-1 0 1]*arrowheadlength,'col','k','linew',2)
    text(Xaxis(1,2),1,-2,'X','fonts',fs)
    % Mark Y axis
    Yaxis = [0 0;-1 1;0 0]*4 + [0 0; 1 1 ; 0 0];
    plotvector(Yaxis,'k-','linew',2)
    % Add on arrowhead
    plot3( [-1 0 1]*arrowheadlength , Yaxis(2,2) + [-1 0 -1]*arrowheadlength, Yaxis(3,2) * [1 1 1],'col','k','linew',2)
    text(0,Yaxis(2,2)+1,-2,'Y','fonts',fs)
    % Mark Z axis
    Zaxis = [0 0;0 0;-1 20];
    plotvector(Zaxis,'k-','linew',2)
    % Add on arrowhead
    plot3( [-1 0 1]*arrowheadlength , [0 0 0], Zaxis(3,2) + [-1 0 -1]*arrowheadlength,'col','k','linew',2)
    text(0,1,Zaxis(3,2),'Z','fonts',fs)

    % Mark on horiozntal plane (XZ plane)
    patch([-1 1 1 -1 -1]*(eyes.radius*1.5+interoculardistance/2),[0 0 0 0 0],[-1 -1 0 0 -1]*eyes.radius*1.5 + [0 0 1 1 0]*15,[0.8 0.8 0.8],'edgecol','none');

    % Draw the eyes
    DrawEyes(eyes);
    % Mark on coordinate system:
    eyecoordhandles = DrawEyeCoords(eyes,coordsystem);

    % Plot object in space, P:
    plotvector(objectP,'marker','o','markerfacecol','k','col','k')
    % Calculate
    imagecoords = ProjectToRetina(eyes,objectP');
    plotvector(imagecoords.headcentered{1},'ro','markerfacecol','r')
    plotvector(imagecoords.headcentered{2},'bo','markerfacecol','b')
    plotvector([imagecoords.headcentered{1} ;objectP]','r-')
    plotvector([imagecoords.headcentered{2} ;objectP]','b-')
    %
    camup([0 1 0])
    camtarget([ 0.105037957554784   1.443297682849908   1.973899533443569])
    campos([-10     5   -15])
    camva(32)
    lighting phong
    light('pos',Zaxis(:,2)/2+Yaxis(:,2)*2)
    camlight
    %    set(gca,'Projection','persp') % NB apparently some Matlab bug --
    %    persepctive rendering gets the intersections wrong, such that the
    %    horizontal plane no longer intersects the horizontal meridians. Leaves
    %    as orthogonal.
    axis equal off

    % Now show both retinae superimposed and viewed from behind:
    ax(jcoordsystem+2)=subplot(2,2,jcoordsystem+2);
    DrawEyes(jointeyes);
    % Mark on coordinate system:
    eyecoordhandles = DrawEyeCoords(jointeyes,coordsystem);
    % For some reason the vertical meridian isn't showing up, so mark it
    % on by hand:
    plot3([0 0 ],[-1 1]*eyes.radius,-[1 1]*eyes.radius,'col',eyes.elcoordcol,'linew',coordsystem.linewidththick);
    % Mark on positions of images: left eye
    images = [sin(imagecoords.AzLat) ; sin(imagecoords.ElLat) ; -cos(imagecoords.AzLat).*cos(imagecoords.ElLat)]*eyes.radius;
    plotvector(images(:,1),'ro','markerfacecol','r')
    plotvector(images(:,2),'bo','markerfacecol','b')

    camup([0 1 0])
    camtarget([ 0 0 0])
    campos([0   0 -25])
    camva(7)
    lighting phong
    camlight headlight
    camlight right
    axis equal off

    set(gcf,'color','w')
end

x1 = 0.05; x2 = 0.53;
y2 = 0.5;
y1 = 0.05;
w = 0.35;
h = 0.35;
axes(ax(1))
text(8,5,0,'a','fonts',fsa)
text(8,9,0,'Elevation-longitude','fonts',fst)
set(ax(1),'pos',[x1 y2 w h])
axes(ax(2))
text(8,5,0,'b','fonts',fsa)
text(8,9,0,'Elevation-latitude','fonts',fst)
set(ax(2),'pos',[x2 y2 w h])
axes(ax(3))
text(2,2,0,'c','fonts',fsa)
set(ax(3),'pos',[x1 y1 w h])
axes(ax(4))
text(2,2,0,'d','fonts',fsa)
set(ax(4),'pos',[x2 y1 w h])


print -dtiff -r800 Fig_VDispDefinition

function imagecoords = ProjectToRetina(eyes,point)
% 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.


for thiseye=1:length(eyes.eyepos);
    n = eyes.eyepos(thiseye).nodal;
    Rot = RotationMatrix(eyes.eyepos(thiseye));
    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) = atan(xonf);
    imagecoords.AzLat(thiseye) = atan(xonf/sqrt(yonf^2+1));
    imagecoords.ElLong(thiseye) = atan(yonf);
    imagecoords.ElLat(thiseye) = atan(yonf/sqrt(xonf^2+1));
    % Position of the image on the rotated eye, in head-centered
    % coords:
    lam = eyes.radius/norm(point-n);
    imagecoords.headcentered{thiseye} = n'-lam*(point-n)';
end

