function eyecoordhandles = DrawEyeCoords(eyes,coordsystem)
% Draw an appropriate coordinate system on the spherical eyeballs.

LEFT=1;RIGHT=2;
h=[];
for thiseye=1:length(eyes.eyepos)
    Rot{thiseye}=RotationMatrix(eyes.eyepos(thiseye));
end
% Fields of coordsystem:
%    .nlines = number of gridlines
%    .spacing = spacing in between each line
n = coordsystem.nlines;
gridlinepos = ([1:n]-(n+1)/2)*coordsystem.spacing;
gridline = [ min(gridlinepos) : range(gridlinepos)/100 : max(gridlinepos) ];
switch coordsystem.name
    case 'Cartesian'
        for k=1:2 % First draw lines of varying x, constant y; Next draw lines of varying y, constant x;
            if length(gridline>0)
                for j=1:n
                    % First get positions on cyclopean eye centered at origin, in primary position:
                    nrm = sqrt(gridline.^2 + gridlinepos(j)^2 + eyes.radius^2 );
                    if k==1
                        % First draw lines of varying x, constant y:
                        X = eyes.radius * gridline./ nrm;
                        Y = eyes.radius * gridlinepos(j)./ nrm;
                        col = eyes.azcoordcol;
                    else
                        % Next draw lines of varying y, constant x:
                        Y = eyes.radius * gridline./ nrm;
                        X = eyes.radius * gridlinepos(j)./ nrm;
                        col = eyes.elcoordcol;
                    end
                    Z = -eyes.radius^2./ nrm;
                    if sind(gridlinepos(j))==0
                        lw = coordsystem.linewidththick;
                    else
                        lw = coordsystem.linewidth;
                    end
                    % Now convert to actual eye position:
                    for thiseye=1:length(eyes.eyepos)
                        tmp = coordsystem.sca * Rot{thiseye} * [X;Y;Z] + repmat(eyes.eyepos(thiseye).nodal,1,length(X));
                        hold on
                        h = [h plot3(tmp(1,:),tmp(2,:),tmp(3,:),'col',col,'linew',lw)];
                    end
                end
            end
        end

    case 'CartesianVirtualPlane'
        for k=1:2 % First draw lines of varying x, constant y; Next draw lines of varying y, constant x;
            if length(gridline>0)
                for j=1:n
                    % First get positions on cyclopean eye centered at origin, in primary position:
                    nrm = sqrt(gridline.^2 + gridlinepos(j)^2 + eyes.radius^2 );
                    if k==1
                        % First draw lines of varying x, constant y:
                        X = eyes.radius * gridline./ nrm;
                        Y = eyes.radius * gridlinepos(j)./ nrm;
                        col = eyes.azcoordcol;
                    else
                        % Next draw lines of varying y, constant x:
                        Y = eyes.radius * gridline./ nrm;
                        X = eyes.radius * gridlinepos(j)./ nrm;
                        col = eyes.elcoordcol;
                    end
                    if sind(gridlinepos(j))==0
                        lw = coordsystem.linewidththick;
                    else
                        lw = coordsystem.linewidth;
                    end
                    % Plot lines on virtual plane
                    hold on
                    Z = -eyes.radius^2./ nrm;
                    lam = -eyes.radius./Z;
                    plot3(lam.*X,lam.*Y,-eyes.radius*ones(size(gridline)),'col',col,'linew',lw)
                    % Now convert to actual eye position:
                    for thiseye=1:length(eyes.eyepos)
                        tmp = coordsystem.sca * Rot{thiseye} * [X;Y;Z] + repmat(eyes.eyepos(thiseye).nodal,1,length(X));
                        hold on
                        h = [h plot3(tmp(1,:),tmp(2,:),tmp(3,:),'col',coordsystem.retcol,'linew',lw)];
                    end
                end
            end
        end

    case 'azlongellong'
        for k=1:2 % First draw lines of varying x, constant y; Next draw lines of varying y, constant x;
            if length(gridline>0)
                for j=1:n
                    % First get positions on cyclopean eye centered at origin, in primary position:
                    nrm = sqrt(tand(gridline).^2 + tand(gridlinepos(j))^2 + 1 );
                    if k==1
                        % First draw lines of varying x, constant y:
                        X = eyes.radius * tand(gridline)./ nrm;
                        Y = eyes.radius * tand(gridlinepos(j))./ nrm;
                        col = eyes.azcoordcol ;
                    else
                        % Next draw lines of varying y, constant x:
                        Y = eyes.radius * tand(gridline)./ nrm;
                        X = eyes.radius * tand(gridlinepos(j))./ nrm;
                        col = eyes.elcoordcol ;
                    end
                    Z = -eyes.radius./ nrm * sign(cosd(gridlinepos(j)));
                    if sind(gridlinepos(j))==0
                        lw = coordsystem.linewidththick;
                    else
                        lw = coordsystem.linewidth;
                    end
                    % Now convert to actual eye position:
                    for thiseye=1:length(eyes.eyepos)
                        tmp = coordsystem.sca * Rot{thiseye} * [X;Y;Z] + repmat(eyes.eyepos(thiseye).nodal,1,length(X));
                        hold on
                        h = [h plot3(tmp(1,:),tmp(2,:),tmp(3,:),'col',col,'linew',lw)];
                    end
                end
            end
        end
    case 'azlongellat'
        for k=1:2 % First draw lines of varying x, constant y; Next draw lines of varying y, constant x;
            if length(gridline>0)
                for j=1:n
                    % First get positions on cyclopean eye centered at origin, in primary position:
                    if k==1
                        % First draw lines of varying az, constant el:
                        az = gridline;
                        el = gridlinepos(j)*ones(1,length(gridline));
                        col = eyes.azcoordcol;
                    else
                        % Next draw lines of varying el, constant az:
                        az = gridlinepos(j)*ones(1,length(gridline));
                        el = gridline;
                        col = eyes.elcoordcol;
                    end
                    X = eyes.radius * cosd(el) .* sind(az);
                    Y = eyes.radius * sind(el) ;
                    Z = -eyes.radius * cosd(el) .* cosd(az);
                    if abs(sind(gridlinepos(j)))<1e-10
                        lw = coordsystem.linewidththick;
                    else
                        lw = coordsystem.linewidth;
                    end
                    % Now convert to actual eye position:
                    for thiseye=1:length(eyes.eyepos)
                        tmp = coordsystem.sca * Rot{thiseye} * [X;Y;Z] + repmat(eyes.eyepos(thiseye).nodal,1,length(X));
                        hold on
                        h = [h plot3(tmp(1,:),tmp(2,:),tmp(3,:),'col',col,'linew',lw)];
                    end
                end
            end
        end
    case 'azlatellong'
        for k=1:2 % First draw lines of varying x, constant y; Next draw lines of varying y, constant x;
            if length(gridline>0)
                for j=1:n
                    % First get positions on cyclopean eye centered at origin, in primary position:
                    if k==1
                        % First draw lines of varying az, constant el:
                        az = gridline;
                        el = gridlinepos(j)*ones(1,length(gridline));
                        col = eyes.azcoordcol;
                    else
                        % Next draw lines of varying el, constant az:
                        az = gridlinepos(j)*ones(1,length(gridline));
                        el = gridline;
                        col = eyes.elcoordcol;
                    end
                    X = eyes.radius .* sind(az);
                    Y = eyes.radius * cosd(az).*sind(el) ;
                    Z = -eyes.radius * cosd(az) .* cosd(el);                    
                    if abs(sind(gridlinepos(j)))<1e-10
                        lw = coordsystem.linewidththick;
                    else
                        lw = coordsystem.linewidth;
                    end
                    % Now convert to actual eye position:
                    for thiseye=1:length(eyes.eyepos)
                        tmp = coordsystem.sca * Rot{thiseye} * [X;Y;Z] + repmat(eyes.eyepos(thiseye).nodal,1,length(X));
                        hold on
                        h = [h plot3(tmp(1,:),tmp(2,:),tmp(3,:),'col',col,'linew',lw)];
                    end
                end
            end
        end
    case 'azlatellat'
        for k=1:2 % First draw lines of varying x, constant y; Next draw lines of varying y, constant x;
            if length(gridline>0)
                for j=1:n
                    % First get positions on cyclopean eye centered at origin, in primary position:
                    if k==1
                        % First draw lines of varying az, constant el:
                        az = gridline;
                        el = gridlinepos(j)*ones(1,length(gridline));
                        col = eyes.azcoordcol;
                    else
                        % Next draw lines of varying el, constant az:
                        az = gridlinepos(j)*ones(1,length(gridline));
                        el = gridline;
                        col = eyes.elcoordcol;
                    end
                    X = eyes.radius .* sind(az);
                    Y = eyes.radius .* sind(el) ;
                    Z = -eyes.radius * sqrt(cosd(az).^2 - sind(el).^2);        
                    if abs(sind(gridlinepos(j)))<1e-10
                        lw = coordsystem.linewidththick;
                    else
                        lw = coordsystem.linewidth;
                    end
                    % Remove bits that go out of range
                    nrm = X.^2+Y.^2+Z.^2;
                    jjj = find(nrm<=eyes.radius^2 + 1e-6 & abs(X)<=eyes.radius & abs(Y)<=eyes.radius & abs(Z)<eyes.radius & imag(Z)==0);
                    X=X(jjj);
                    Y=Y(jjj);
                    Z=Z(jjj);
                    % Now convert to actual eye position:
                    for thiseye=1:length(eyes.eyepos)
                        tmp = coordsystem.sca * Rot{thiseye} * [X;Y;Z] + repmat(eyes.eyepos(thiseye).nodal,1,length(X));
                        hold on
                        h = [h plot3(tmp(1,:),tmp(2,:),tmp(3,:),'col',col,'linew',lw)];
                    end
                end
            end
        end


end

% Mark fovea
if isfield(coordsystem,'markfovea')
    for thiseye=1:length(eyes.eyepos)
        tmp = coordsystem.sca * Rot{thiseye} * [0;0;-eyes.radius] + eyes.eyepos(thiseye).nodal;
        h = [h plot3(tmp(1,:),tmp(2,:),tmp(3,:),'col',eyes.azcoordcol,'marker','o','markerfacecol',eyes.azcoordcol,'markersize',12)];
    end
end
if isfield(coordsystem,'markpupil')
    for thiseye=1:length(eyes.eyepos)
        tmp = coordsystem.sca * Rot{thiseye} * [0;0;eyes.radius] + eyes.eyepos(thiseye).nodal;
        h = [h plot3(tmp(1,:),tmp(2,:),tmp(3,:),'col',eyes.azcoordcol,'marker','o','markerfacecol',eyes.azcoordcol,'markersize',12)];
    end
end

eyecoordhandles = h;

