clear all
close all
global L R CYC DELTA HALFDELTA ALL
I = 6.25; % interocualr distance cm
L = 1; R = 2; CYC = 3; DELTA = 4; HALFDELTA = 5; ALL = 6;

Hcdeg = -5;
Vcdeg = 2;
Tcdeg = 0;
HDdeg = 15;
VDdeg = 0;
TDdeg = 0;
eyes = GetEyesGivenCV(I,Hcdeg,Vcdeg,Tcdeg,HDdeg,VDdeg,TDdeg);

eyes(ALL).eyeballcolor = [0.95 0.95 1];
eyes(ALL).eyeballtransparency = 1;
% Colors in which to draw coordinate systems
eyes(ALL).coordcol = 'k' ;
eyes(ALL).nsurface = 50;

% Specify the coordinate system to be painted on the retina
coordsystem.name = 'azlongellat';
coordsystem.nlines = 13;
coordsystem.spacing = 15;
%%
coordsystem.linewidth = 0.5;
coordsystem.linewidththick = 2;
coordsystem.sca  = 0.995; % amount by which to exaggerate eyeball radius so coord lines are visible

figure('pos',[300   198   576   399])
hold on
% Generate a hemisphere to represent the retina
nn=eyes(ALL).nsurface;
[Xret Yret Zret] = sphere(nn);
Zret = Zret( 1:ceil((nn+1)/2),:);
Xret = Xret( 1:ceil((nn+1)/2),:);
Yret = Yret( 1:ceil((nn+1)/2),:);

% Define a surface

% Calculate exact disparity field
[ny,nx]= size(Xret);
for jx=1:nx
    for jy=1:ny
        % head-centered coords
        X = Xret(jy,jx);
        Y = Yret(jy,jx);
        Z = Zret(jy,jx);
        % Convert to coordinates on the cyclopean retina
        f = sqrt(X^2+Y^2+Z^2);
        kappa(jy,jx) = asin(Y/f);
        alpha(jy,jx) = asin(X/f/cos(kappa(jy,jx)));

    end
end
% define a surface:
S = 90 + 0.05 * sin(kappa) + 0.01 * cos(alpha);
U = -S.*tan(alpha);
W = -S.*tan(kappa)./cos(alpha);
Xsurf = U * eyes(CYC).m1(1) + W * eyes(CYC).m2(1) + S * eyes(CYC).m3(1);
Ysurf = U * eyes(CYC).m1(2) + W * eyes(CYC).m2(2) + S * eyes(CYC).m3(2);
Zsurf = U * eyes(CYC).m1(3) + W * eyes(CYC).m2(3) + S * eyes(CYC).m3(3);
[alphaDelapprox,kappaDelapprox] = ApproxDisparity(eyes,'azimuth-longitude/elevation-latitude','retina, S',alpha,kappa,S);

fvc = surf2patch(Xret,Yret,Zret,kappaDelapprox);
retina = patch(fvc,'edgecolor','none','facealpha',eyes(ALL).eyeballtransparency,'edgealpha',0);
%fvc = surf2patch(Xsurf,Ysurf,Zsurf,1);
%surface = patch(fvc,'facecolor',[1 1 1]*0.5,'edgecolor','none','facealpha',1,'edgealpha',0);
shading faceted; view(3)

% Show the lat & long

n = coordsystem.nlines;
gridlinepos = ([1:n]-(n+1)/2)*coordsystem.spacing;
gridline = [ min(gridlinepos) : range(gridlinepos)/100 : max(gridlinepos) ];
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
            if k==1
                % First draw lines of varying az, constant el:
                az = gridline;
                el = gridlinepos(j)*ones(1,length(gridline));
            else
                % Next draw lines of varying el, constant az:
                az = gridlinepos(j)*ones(1,length(gridline));
                el = gridline;
            end
            X = coordsystem.sca*cosd(el) .* sind(az);
            Y = coordsystem.sca*sind(el) ;
            Z = -coordsystem.sca*cosd(el) .* cosd(az);
            plot3(X,Y,Z,'col',eyes(ALL).coordcol,'linew',coordsystem.linewidth)
            if gridlinepos(j)==0
                plot3(X,Y,Z,'col',eyes(ALL).coordcol,'linew',coordsystem.linewidththick)
            end
        end
    end
end

camup([0 1 0])
camtarget([-0.0037    0.0503   -0.0199])
campos([-10 4 13])
camva(6)
set(gca,'Projection','persp')
axis equal off
set(gcf,'color','w')

set(retina,'BackFaceLighting','unlit','FaceLighting','phong','FaceColor','interp')

az = 0; 
el = 0;
azcol = 'k';
elcol = 'k';
fs = 20;
lw= 4;

% Draw arrows
% Differentiating wrt azimuth:
rlen = 20; % length of arrow to draw, in deg
azz = az+[0:0.01:1]*rlen;
ell = el*ones(size(azz));
r = 0.99;
X = r*cosd(ell) .* sind(azz);
Y = r*sind(ell) ;
Z = -r*cosd(ell) .* cosd(azz);
plot3(X,Y,Z,'col',azcol,'linew',lw)
x = 0.05; % length of arrow head
% Draw arrow head
plot3(X(end)-x*[1 0 1],Y(end)-x*[-1 0 1],Z(end)*[1 1 1],'col',azcol,'linew',lw)
% Label
text(X(end)-0.05,Y(end)+0.12,Z(end),'\partial/\partial\alpha','col',azcol,'fontw','bold','fonts',fs)
% Differentiating wrt elevation:
rlen = 20; % length of arrow to draw, in deg
ell = el+[0:0.01:1]*rlen;
azz = az*ones(size(ell));
X = r*cosd(ell) .* sind(azz);
Y = r*sind(ell) ;
Z = -r*cosd(ell) .* cosd(azz);
plot3(X,Y,Z,'col',elcol,'linew',lw)
% Draw arrow head
plot3(X(end)+x*[-1 0 1],Y(end)+x*[-1 0 -1],Z(end)*[1 1 1],'col',elcol,'linew',lw)
% Label
text(X(end)-0.05,Y(end)+0.1,Z(end)+0.15,'\partial/\partial\kappa','col',elcol','fontw','bold','fonts',fs)

% Plot gaze axis:
plot3([0 0 ],[0 0 ],[-1 16],'k','linew',lw)
zlim([-1 1])
plot3(0,0,0,'ko','markersize',12,'markerfacecol','k')


print -dtiff -r700 -opengl Fig_DifferentiatingAtFovea