function tmp
clear all
close all
fs = 16; % font size
zerocol = [1 1 0.99]; % virtually white, but for some reason if you put 'w' it prints out black with "print" command
linew = 2; % color for plotting loci of zero

figure('pos',[5   278   873   420])

% Draw the cyclopean retina from behind
% (X,Y,Z) is the position on the cyclopean retina in Cartesian coordinates, i.e.
% X^2+Y^2+Z^2 = radius^2 of spherical retina
X=[-1:0.005:1];
Y=[-1:0.005:1];
for jX=1:length(X)
    for jY =1:length(Y)
        if X(jX)^2+Y(jY)^2<1
            % First convert (X,Y,Z) into lat and long on retina
            kappa  = asin(Y(jY));
            alpha = asin(X(jX)/cos(kappa));
            % Now apply expressions for expected vertical disparity:
            meankappaDel(jY,jX) = sin(kappa).*sin(alpha);
            meanetaDel(jY,jX) = tan(kappa).*sin(alpha)./(cos(alpha)^2+tan(kappa)^2);
        else
            % Disparity is NaN if we are off the retina
            meankappaDel(jY,jX) = NaN;
            meanetaDel(jY,jX) = NaN;
        end
    end
end

% Set colormap same for all plots, based on limits of
% elevation-latitude disparity (kappaDel) for sphere at R=70.
mn=min(min([meankappaDel]))
mx=max(max([meankappaDel]))

for j=1:2 % el-long, then el-lat vdisp
    ax(j)=subplot(1,2,j);
    % imagesc, except replace NaNs with white:
    if j==2
        im=meankappaDel;
    else
        im=meanetaDel;
    end
    colim = GetColor(im,mn,mx);
    image(X,Y,colim)
    hold on
    contour(X,Y,im,[-1:0.1:1],'col',zerocol)
    axis equal tight off
    % Mark fovea
    plot(0,0,'ko','markerfacecol','k','markersize',8)

    if j==1
        title('mean elevation longitude < \eta_\Delta >','fonts',fs)
    elseif j==2
        title('mean elevation latitude < \kappa_\Delta >','fonts',fs)
    end

    % Mark on lines of lat and long:
    hold on
    alpha=[-90:15:90]*pi/180;
    kappa=[-90:90]*pi/180;
    for ja=1:length(alpha)
        XlineAzlong = cos(kappa).*sin(alpha(ja));
        YlineAzlong = sin(kappa);
        plot(XlineAzlong,YlineAzlong,'k')
    end
    if j==2
        alpha=[-90:90]*pi/180;
        kappa=[-90:15:90]*pi/180;
        for jk=1:length(kappa)
            XlineEllat = cos(kappa(jk)).*sin(alpha);
            YlineEllat = sin(kappa(jk)).*ones(size(alpha));
            plot(XlineEllat,YlineEllat,'k')
        end
    elseif j==1
        alpha=[-90:90]*pi/180;
        eta=[-90:15:90]*pi/180;
        for j=1:length(eta)
            XlineEllong = tan(alpha)./sqrt(1+tan(alpha).^2+tan(eta(j)).^2);
            YlineEllong = tan(eta(j))./sqrt(1+tan(alpha).^2+tan(eta(j)).^2);
            plot(XlineEllong,YlineEllong,'k')
        end
    end
end
 set(ax(1),'pos',[0.005 0.05 0.4 0.8])
 set(ax(2),'pos',[0.41 0.05 0.4 0.8])
 % Do colorbar
axes('pos',[0.94 0.11 0.05 0.7876])
tmpim = [0:63]'/63*(mx-mn)+mn;
colim = GetColor(tmpim,mn,mx);
image(0,tmpim,colim);
set(gca,'ydir','normal','xticklabel','','xtick',[],'fonts',fs,'ytick',[-1:0.2:1])
ylabel('vertical disparity / <H_\Delta>','fonts',fs)
print -dtiff -r700 -zbuffer Fig_ExpectedVDisp

function colim = GetColor(im,mn,mx);
% Get pseudocolor
% NB I altered default color scale in order to display range better
map = colormap(jet);
nc = 64;
% Here scale the image in order to concentrate on |values|<0.2:
scledim = cumulativegaussian(im/mx,0,0.2);
%scledim = 0.5 + sign(im).*(abs(im)/mx).^0.42/2;
% FInd which color to use
indx = round(scledim*(nc-1))+1;
% Saturate at top and bottom of colorscale:
indx(indx<1)=1;
indx(indx>64)=64;
[ny,nx]=size(im);
colim = zeros(ny,nx,3);
nancol =  get(gcf,'color');
for jy=1:ny
    for jx=1:nx
        if isnan(indx(jy,jx))
            colim(jy,jx,:) = [1 1 1];
        else
            colim(jy,jx,:) = map(indx(jy,jx),:);
        end
    end
end


function g = cumulativegaussian(x,mu,sig)
g = 0.5 * (1 + erf((x-mu)./sqrt(2)./sig));
