function hd = coneAnalysis
  
%   datachoice = 'visual field';  % for main figure
  datachoice = 'retina';  % for supplementary figure


  warning('off','MATLAB:interp1:UsePCHIP')
  % deactivating the following cbrewer warning:
  % Warning: INTERP1(...,'CUBIC') will change in a future release. Use INTERP1(...,'PCHIP') instead. 
  % > In interp1>sanitycheckmethod (line 265)
  %   In interp1>parseinputs (line 401)
  %   In interp1 (line 80)
  %   In interpolate_cbrewer (line 31)
  %   In cbrewer (line 101)
  %   In heatmap20180823 (line 21) 
  
%   global globalparam
  close all
  
  cd(['..\..\data\Fig6 cone data\'])
  
    
  % Define which character sequences identify the different types of data file.
  filestring.all   = 'allCones';
  filestring.blue  = 'Blue';
  filestring.green = 'Green';
  filestring.red   = 'Red';
  filestring.uv    = 'UV';
  
  
  % Manually set max. cone density per square mm, as in Zimmermann et al.
  maxdensity = 35000;
  
  % Manually select data files, and automatically check their names.
  switch datachoice
    case 'visual field'
      fileselection = {'VisualFieldProjection_Red_8dpf.txt', ...
                       'VisualFieldProjection_Green_8dpf.txt','VisualFieldProjection_Blue_8dpf.txt', ...
                       'VisualFieldProjection_UV_8dpf.txt','VisualFieldProjection_allCones_8dpf.txt'};
    case 'retina'
      fileselection = {'DensityinRetina_Red_8dpf.txt', ...
                       'DensityinRetina_Green_8dpf.txt','DensityinRetina_Blue_8dpf.txt', ...
                       'DensityinRetina_UV_8dpf.txt','DensityinRetina_allCones_8dpf.txt'};
  end
  occurrence = struct('all',[],'blue',[],'green',[],'red',[],'uv',[]);
  for file = 1:numel(fileselection)
    occurrence.all   = [occurrence.all,   ~isempty(strfind(fileselection{file},filestring.all))];
    occurrence.blue  = [occurrence.blue,  ~isempty(strfind(fileselection{file},filestring.blue))];
    occurrence.green = [occurrence.green, ~isempty(strfind(fileselection{file},filestring.green))];
    occurrence.red   = [occurrence.red,   ~isempty(strfind(fileselection{file},filestring.red))];
    occurrence.uv    = [occurrence.uv,    ~isempty(strfind(fileselection{file},filestring.uv))];
  end
  
  
  % Check whether each type of file was selected exactly once.
  if sum(occurrence.all) == 0
    disp(['None of the file names contain the characters ''',filestring.all,'''.'])
  elseif sum(occurrence.all) > 1
    error(['More than one file name contains the characters ''',filestring.all,'''!'])
  else
    datafile.all = fileselection{find(occurrence.all)}; %#ok<FNDSB>
  end
  if sum(occurrence.blue) == 0
    disp(['None of the file names contain the characters ''',filestring.blue,'''.'])
  elseif sum(occurrence.blue) > 1
    error(['More than one file name contains the characters ''',filestring.blue,'''!'])
  else
    datafile.blue = fileselection{find(occurrence.blue)}; %#ok<FNDSB>
  end
  if sum(occurrence.green) == 0
    disp(['None of the file names contain the characters ''',filestring.green,'''.'])
  elseif sum(occurrence.green) > 1
    error(['More than one file name contains the characters ''',filestring.green,'''!'])
  else
    datafile.green = fileselection{find(occurrence.green)}; %#ok<FNDSB>
  end
  if sum(occurrence.red) == 0
    disp(['None of the file names contain the characters ''',filestring.red,'''.'])
  elseif sum(occurrence.red) > 1
    error(['More than one file name contains the characters ''',filestring.red,'''!'])
  else
    datafile.red = fileselection{find(occurrence.red)}; %#ok<FNDSB>
  end
  if sum(occurrence.uv) == 0
    disp(['None of the file names contain the characters ''',filestring.uv,'''.'])
  elseif sum(occurrence.uv) > 1
    error(['More than one file name contains the characters ''',filestring.uv,'''!'])
  else
    datafile.uv = fileselection{find(occurrence.uv)}; %#ok<FNDSB>
  end




%   % Fit centres from FitAzimElev20181217.txt:
%   leftnormmaxazim = 74.1675;
%   leftnormmaxelev = 5.9175;
%   leftnormmeanazim = 74.9488;
%   leftnormmeanelev = 5.7009;
%   leftnormratioazim = 74.9490;
%   leftnormratioelev = 5.7007;
%   rightnormmaxazim = -72.0633;
%   rightnormmaxelev = 7.3858;
%   rightnormmeanazim = -71.1931;
%   rightnormmeanelev = 7.5766;
%   rightnormratioazim = -71.1933;
%   rightnormratioelev = 7.5767;
  leftazim_body = -80.2647;
  leftelev_body = 6.0557;
  rightazim_body = +77.0306;
  rightelev_body = -1.9776;

  % Median eye position, pooled once, from positionAnalysis20181217.m:
  leftmedianazim =  +6.2450;
  rightmedianazim = -9.8760;
  leftmeanazim =  +5.2455;
  rightmeanazim = -9.8840;
  leftstdazim = 6.1703;
  rightstdazim = 6.5227;
  leftsemazim = 0.2472;
  rightsemazim = 0.2613;

  % Median vertical eye position, from figureVertical20190225.m:
  leftmeanelev = 3.4937;  % meanmeanposition.UUH(:,:,1);
  rightmeanelev = 4.8735; % meanmeanposition.UUH(:,:,2);
  
  % Corrected to eye-centred visual field angles:
  leftazim_eye = abs(leftazim_body - leftmeanazim);
  leftelev_eye = leftelev_body - leftmeanelev;
  rightazim_eye = abs(rightazim_body - rightmeanazim);
  rightelev_eye = rightelev_body - rightmeanelev;




  % Create one figure panel per file, and all plots on it.
  datapresent = fieldnames(datafile);
  % Reorder field names to reflect desired arrangement of panels on figure.
  datapresent = {datapresent{[4 3 2 5 1]}};
  
  for typeindex = 1:numel(datapresent)
    
    conetype = datapresent{typeindex};
    textfile = datafile.(conetype);

    switch datachoice
      case 'visual field'
        conedata = importdata(textfile) * maxdensity;
      case 'retina'
        conedata = importdata(textfile);
    end
    switch datachoice
      case 'visual field'
        conedata = rot90(conedata,3); % so up is at top, anterior to the left
      case 'retina'
        conedata = rot90(conedata,1); % so up is at bottom, anterior to the right  
    end
    conedata(conedata==0) = NaN; % consider nonzero data only, i.e. from the true area of the retina

    
    switch datachoice
      case 'visual field'
        % Undo approximate correction by Takeshi, then apply a more precise one:
        offsetAzimTakeshi = 8;
        offsetElevTakeshi = 5;
        offsetAzimAccurate = abs(leftmedianazim);
        offsetElevAccurate = leftmeanelev;
        resolutionAzim = size(conedata,2)/360;
        resolutionElev = size(conedata,1)/180;
        discreteShiftAzim = round((offsetAzimTakeshi-offsetAzimAccurate)/resolutionAzim);
        discreteShiftElev = round((offsetElevTakeshi-offsetElevAccurate)/resolutionElev);
        nandata = NaN(size(conedata));
        % Execute shift (min/max sanitise indices in case shift is negative):
        nandata(max(1,1-discreteShiftElev):min(end,end-discreteShiftElev), ...
                max(1,1-discreteShiftAzim):min(end,end-discreteShiftAzim)) = ...
          conedata(max(1,1+discreteShiftElev):min(end,end+discreteShiftElev), ...
                   max(1,1+discreteShiftAzim):min(end,end+discreteShiftAzim));
        % Overwrite original data:
        conedata = nandata;
      case 'retina'
        % Cannot apply correction, as VF filter used by Takeshi is unknown.
    end


  %   cd(['C:\Users\dehmelt\Dropbox\Documents\MATLAB\', ...
  %       '201808 SphereAnalysis'])
  % 
  %   numelevation = size(conedata,1);
  %   numazimuth = size(conedata,2);
  % 
  %   % The following four values represent the coordinates of the four corners
  %   % of the 2D data sets provided by Baden lab (and are currently unknown).
  %   elevation.max = 90;
  %   elevation.min = -90;
  %   azimuth.max = 180;
  %   azimuth.min = -180;

  %   % Build an (NX*NY x 2) array, where NX is the number of azimuth values,
  %   % NY is the number of elevation values, NX*NY is the total number of data
  %   % points, and for each data point (each row), the azimuth (column 1) and
  %   % elevation (column 2) are listed.
  %   elevation.step = (elevation.max-elevation.min)/(numelevation-1);
  %   azimuth.step = (azimuth.max-azimuth.min)/(numazimuth-1);
  %   elevationlist = elevation.max:-elevation.step:elevation.min;
  %   azimuthlist = azimuth.max:-azimuth.step:azimuth.min;
  %   predictor = [];
  %   for k = 1:numazimuth
  %     predictor = [predictor; ...
  %                  repmat(azimuthlist(k),[numelevation 1]), elevationlist'];
  %   end

  %   % One raw data value per data point, might have to be reformatted to match
  %   % the order in which the list of data azimuths/elevations in "predictor" is
  %   % being created. Predictor goes top-to-bottom, left-to-right.
  %   response = conedata(:);
  % 
  %   offset = 0;
  %   amplitude = 1;
  %   meanazimuth = 90;
  %   meanelevation = 0;
  %   concentration = 1;
  %   
  %   globalparam.offset = offset;
  %   globalparam.amplitude = amplitude;
  %   globalparam.concentration = concentration;
  % 
  % 
  %   figure(1)
  % %   imagesc(conedata)
  %   surf(conedata)
  %   shading flat
  %   axis square
  %   colormap(cbrewer('div','PRGn',100))
  %   hold on
  %   plot([1,numazimuth],[numelevation/2 numelevation/2],'--w','LineWidth',2)
  %   plot([numazimuth/2, numazimuth/2],[1,numelevation],'--w','LineWidth',2)
  %   hold off
  %   set(gca,'XTick',[1,numazimuth*[1/4,1/2,3/4,1]],'XTickLabel',{180 90 0 -90 -180})
  %   set(gca,'YTick',[1,numelevation*[1/4,1/2,3/4,1]],'YTickLabel',{90 45 0 -45 -90})
  %   xlabel('azimuth'),ylabel('elevation'),set(gcf,'Color',[1 1 1])
  %   dataclim = get(gca,'CLim');


    hd.cone.fg = figure(1);
    hd.cone.fg.Position = [400 160 1000 700];
%     hd.cone.fg.Position = [500 50 540 380];
    hd.cone.fg.Color = [1 1 1];

    hd.cone.ax(typeindex,1) = axes;
    if typeindex < 4
      hd.cone.ax(typeindex,1).Position = [.11+(typeindex-1)*.29 .61 hd.cone.fg.Position([4 3])/3500];
    else
      hd.cone.ax(typeindex,1).Position = [.19+(typeindex-4)*.42 .11 hd.cone.fg.Position([4 3])/3500];
    end
    
    az = 0:180/(size(conedata,1)-1):180;
    el = -90:180/(size(conedata,2)-1):90;
  %   conedatanonnan = fliplr(flipud(conedata));
    conedatanonnan = conedata;
    conedatanonnan(isnan(conedatanonnan)) = 0;
    density = conv2(conedatanonnan,fspecial('gaussian',[100 100],3),'same');
  %   surf(az,el,conv2(conedatanonnan,fspecial('gaussian',[100 100],10),'same'))
  %   shading flat
    hold on
    if strcmp(datachoice,'retina') && typeindex == 5 % FD20190829 added as simple correction
      density = fliplr(density);
    end
    imagesc(az,el,density)
%     if strcmp(conetype,'all')
%       contour(az,el,density,'LevelStep',1e-1,'Color',.5*[1 1 1])
%     else
%       contour(az,el,density,'LevelStep',3e-2,'Color',.5*[1 1 1])
%     end
    contourlevels = (.1:.1:1) * max(max(density));  % increments = 10%-of-max.
    contour(az,el,density,contourlevels,'Color',.5*[1 1 1]); 
    hold off
%     axis square
%     uistack(hd.cone.ax(typeindex,2),'bottom')
%     hd.cone.ax(2,2) = copyobj(hd.cone.ax(2,1),hd.cone.fg)

    hd.cone.ax(typeindex,1).XLim = [0,180];
    hd.cone.ax(typeindex,1).YLim = [-90,90];
    oldclim = hd.cone.ax(typeindex,1).CLim;
%     climcutoff = .8;
    climcutoff = .60;
%     climcutoff = 0;
    hd.cone.ax(typeindex,1).CLim = [min(oldclim)+climcutoff*(max(oldclim)-min(oldclim)), ...
                            max(oldclim)];
    switch datachoice
      case 'visual field'
        hd.cone.ax(typeindex,1).XLabel.String = 'azimuth';
        hd.cone.ax(typeindex,1).YLabel.String = 'elevation';
      case 'retina'
        hd.cone.ax(typeindex,1).XLabel.String = 'x  (a.u.)';
        hd.cone.ax(typeindex,1).YLabel.String = 'y  (a.u.)';
    end
    textfilesplit = strsplit(textfile,'_');
    switch textfilesplit{2}
      case 'allCones'
        densitylabel = 'all cones combined';
      case 'Green'
        densitylabel = 'green cones (mm^{-2})';
      case 'Red'
        densitylabel = 'red cones (mm^{-2})';
      case 'Blue'
        densitylabel = 'blue cones (mm^{-2})';
      case 'UV'
        densitylabel = 'ultraviolet cones';
    end
%     switch textfilesplit{2}
%       case 'allCones'
%         densitylabel = {'retinal density','(any cones)'};
%       case 'Green'
%         densitylabel = {'retinal density / mm^{-2}','(green cones)'};
%       case 'Red'
%         densitylabel = {'retinal density','(red cones)'};
%       case 'Blue'
%         densitylabel = {'retinal density','(blue cones)'};
%       case 'UV'
%         densitylabel = {'retinal density','(ultraviolet cones)'};
%     end
%     switch textfilesplit{2}
%       case 'allCones'
%         densitylabel = 'density (any cones)';
%       case 'Green'
%         densitylabel = 'density (green cones)';
%       case 'Red'
%         densitylabel = 'density (red cones)';
%       case 'Blue'
%         densitylabel = 'density (blue cones)';
%       case 'UV'
%         densitylabel = 'density (ultraviolet cones)';
%     end

    hd.cone.ax(typeindex,1).Title.String = densitylabel;
    hd.cone.ax(typeindex,1).Box = 'on';
%     hd.cone.ax(typeindex,1).XTick = 0:15:180;
%     hd.cone.ax(typeindex,1).YTick = -90:15:90;
%     hd.cone.ax(typeindex,1).XTickLabel = {0,[],[],45,[],[],90,[],[],135,[],[],180};
%     hd.cone.ax(typeindex,1).YTickLabel = {-90,[],[],-45,[],[],0,[],[],45,[],[],90};
%     hd.cone.ax(typeindex,1).XTick = 0:45:180;
%     hd.cone.ax(typeindex,1).YTick = -90:45:90;
%     hd.cone.ax(typeindex,1).XTickLabel = {0,45,90,135,180};
%     hd.cone.ax(typeindex,1).YTickLabel = {-90,-45,0,45,90};
    hd.cone.ax(typeindex,1).XTick = 0:30:180;
    hd.cone.ax(typeindex,1).YTick = -90:30:90;
    switch datachoice
      case 'visual field'
        hd.cone.ax(typeindex,1).XTickLabel = {0,[],[],90,[],[],180};
        hd.cone.ax(typeindex,1).YTickLabel = {-90,[],[],0,[],[],90};
      case 'retina'
        hd.cone.ax(typeindex,1).XTickLabel = {0,[],[],[],[],[],1};
        hd.cone.ax(typeindex,1).YTickLabel = {0,[],[],[],[],[],1};
    end
    hd.cone.ax(typeindex,1).FontSize = 14;
    hd.cone.ax(typeindex,1).TitleFontSizeMultiplier = 1;
    hd.cone.ax(typeindex,1).TitleFontWeight = 'bold';
%     hd.cone.ax(typeindex,1).Title.Interpreter = 'none';
    hd.cone.ax(typeindex,1).Title.Interpreter = 'tex';
  %   [hd.cone.ax(typeindex,1).XLabel.FontSize, ...
  %    hd.cone.ax(typeindex,1).YLabel.FontSize, ...
  %    hd.cone.ax(typeindex,1).Title.FontSize, ...
  %    ha.cone.ax(typeindex,1).FontSize,
    hd.cone.ax(typeindex,1).XGrid = 'on';
    hd.cone.ax(typeindex,1).YGrid = 'on';
    hd.cone.ax(typeindex,1).Layer = 'top';
    
    
    hold on
    plot([0,180],[0,0],'--k')
    plot([90,90],[-90,90],'--k')
    hold off

   
    
% %     leftnormmaxazimVF = abs(leftnormmaxazim - leftmedianazim);
% %     leftnormmeanazimVF = abs(leftnormmeanazim - leftmedianazim);
% %     leftnormratioazimVF = abs(leftnormratioazim - leftmedianazim);
% %     rightnormmaxazimVF = abs(rightnormmaxazim - rightmedianazim);
% %     rightnormmeanazimVF = abs(rightnormmeanazim - rightmedianazim);
% %     rightnormratioazimVF = abs(rightnormratioazim - rightmedianazim);



    hold on
    hd.cone.ob(1,1) = plot([153,165],[-65,-65],'k');
    hd.cone.ob(1,2) = plot([159,159],[-59,-71],'k');
    switch datachoice
      case 'visual field'
        hd.cone.ob(1,3) = text(153,-64.5,'A','HorizontalAlignment','right','VerticalAlignment','middle');
        hd.cone.ob(1,4) = text(166,-64.5,'P','HorizontalAlignment','left','VerticalAlignment','middle');
        hd.cone.ob(1,5) = text(159.2,-59.5,'U','HorizontalAlignment','center','VerticalAlignment','bottom');
        hd.cone.ob(1,6) = text(159.2,-69.5,'D','HorizontalAlignment','center','VerticalAlignment','top');
      case 'retina'
        hd.cone.ob(1,3) = text(153,-64.5,'P','HorizontalAlignment','right','VerticalAlignment','middle');
        hd.cone.ob(1,4) = text(166,-64.5,'A','HorizontalAlignment','left','VerticalAlignment','middle');
        hd.cone.ob(1,5) = text(159.2,-59.5,'D','HorizontalAlignment','center','VerticalAlignment','bottom');
        hd.cone.ob(1,6) = text(159.2,-69.5,'U','HorizontalAlignment','center','VerticalAlignment','top');
    end
    hold off
    
%     circle.radius = 40/2;
%     circle.xcentre = leftnormmaxazimVF; % pick preferred measure here
%     circle.ycentre = leftnormmaxelev; % pick preferred measure here
%     circle.x = -circle.radius:.1:circle.radius;
%     circle.y = sqrt(circle.x(1)^2-circle.x.^2);
%     hold on
%     plot([circle.xcentre + circle.x, fliplr(circle.xcentre + circle.x)], ...
%          [circle.ycentre + circle.y, circle.ycentre - circle.y],'-','LineWidth',3,'Color',0*[1 1 1])    
%     plot([circle.xcentre + circle.x, fliplr(circle.xcentre + circle.x)], ...
%          [circle.ycentre + circle.y, circle.ycentre - circle.y],'-','LineWidth',2,'Color',1*[1 1 1])    
%     hold off
    
%     boundary.angle = 64;
% %     centre.azimuth = leftnormmaxazimVF;
% %     centre.elevation = leftnormmaxelev;
%     centre.azimuth = rightazim;
%     centre.elevation = rightelev;
%     boundary.resolution = .01;
%     [boundary.azimuth,boundary.elevation] = ...
%       diskBoundary(centre.azimuth,centre.elevation,boundary.angle,boundary.resolution);
%     hold on
%     plot(boundary.azimuth,boundary.elevation,'-','LineWidth',3,'Color',0*[1 1 1])    
%     plot(boundary.azimuth,boundary.elevation,'-','LineWidth',2,'Color',1*[1 1 1])    
%     hold off

    boundary.angle = 40;
%     centre.azimuth = leftnormmaxazimVF;
%     centre.elevation = leftnormmaxelev;
%     centre.azimuth = leftazim_body;
%     centre.elevation = leftelev_body;
    centre.azimuth = abs(leftazim_body);
    centre.elevation = leftelev_body;
    switch datachoice
      case 'retina'
        centre.azimuth = 180-centre.azimuth;
        centre.elevation = -centre.elevation;
    end
    switch datachoice
      case 'visual field'
        boundary.resolution = .01;
        [boundary.azimuth,boundary.elevation] = ...
          diskBoundary(centre.azimuth,centre.elevation,boundary.angle,boundary.resolution);
        hold on
        plot(boundary.azimuth,boundary.elevation,'-','LineWidth',3,'Color',0*[1 1 1])    
        plot(boundary.azimuth,boundary.elevation,'-','LineWidth',2,'Color',1*[1 1 1])    
        hold off
    end


    mtp = 'o';
    msz = 100;
    mec = [0 0 0];
%     mfc = [.2 .8 .2];
%     mfa = .2;
    mfc = 1*[1 1 1];
    mfa = .5;
    hold on
%     hd.cone.ob(1,7) = scatter(rightazim,rightelev,msz,mtp, ...
%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
    hd.cone.ob(1,8) = scatter(centre.azimuth,centre.elevation,msz,mtp, ...
                               'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',.8);
%     hd.cone.ob(1,9) = scatter(leftnormratioazimVF,leftnormratioelev,msz,mtp, ...
%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
%     hd.cone.ob(1,10) = scatter(rightnormmaxazimVF,rightnormmaxelev,msz,mtp, ...
%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
%     hd.cone.ob(1,11) = scatter(rightnormmeanazimVF,rightnormmeanelev,msz,mtp, ...
%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
%     hd.cone.ob(1,12) = scatter(rightnormratioazimVF,rightnormratioelev,msz,mtp, ...
%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
    hold off
    
    
    
    if typeindex == 2 || typeindex == 3
      hd.cone.ax(typeindex,1).YLabel.String = [];
      hd.cone.ax(typeindex,1).YTickLabel = [];
    end
    switch conetype
      case 'blue'
        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Blues',100))
      case 'green'
        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Greens',100))
      case 'red'
        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Reds',100))
      case 'uv'
        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Purples',100))
      case 'all'
        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Greys',100))
    end
    
    hd.cone.ax(typeindex,2) = axes;
    hd.cone.ax(typeindex,2).Position = hd.cone.ax(typeindex,1).Position + ...
                                       hd.cone.ax(typeindex,1).Position(3) * .28*[1 0 0 0];
    hd.cone.cbar(typeindex) = colorbar;

    hd.cone.ax(typeindex,2).Colormap = hd.cone.ax(typeindex,1).Colormap;
    hd.cone.ax(typeindex,2).CLim     = hd.cone.ax(typeindex,1).CLim;
    hd.cone.ax(typeindex,2).FontSize = hd.cone.ax(typeindex,1).FontSize;
    hd.cone.ax(typeindex,2).Color = 'none';
    [hd.cone.ax(typeindex,2).XAxis.Visible,hd.cone.ax(typeindex,2).YAxis.Visible] = deal('off');
        
        
%     colormap(cbrewer('div','PRGn',100))



  end
  
  hd.cone.ax(1,3) = axes;
  hd.cone.ax(1,3).Position = [8 0 1 1];
  hd.cone.ax(1,3).Color = 'none';
  [hd.cone.ax(1,3).XAxis.Visible,hd.cone.ax(1,3).YAxis.Visible] = deal('off');
  for typeindex = 1:numel(datapresent)
    parentpos = hd.cone.ax(typeindex,1).Position;
    xpos = parentpos(1)+parentpos(3)/2;
    ypos = parentpos(2)+parentpos(4)*1.128;
    hd.cone.ob(3,typeindex) = text(xpos, ypos, 'retinal density / mm^{-2}');
%     hd.cone.ob(3,typeindex) = text(xpos, ypos, 'retinal density (µm^{-2})');
    hd.cone.ob(3,typeindex).FontWeight = 'normal';
    hd.cone.ob(3,typeindex).Interpreter = 'tex';
    hd.cone.ob(3,typeindex).FontSize = hd.cone.ax(typeindex,1).FontSize;
    hd.cone.ob(3,typeindex).HorizontalAlignment = 'center';
    hd.cone.ob(3,typeindex).VerticalAlignment = 'bottom';
  end
  uistack(hd.cone.ax(1,3),'bottom')

  hd.cone.ax(1,4) = axes;
  hd.cone.ax(1,4).Position = [0 0 1 1];
  hd.cone.ax(1,4).Color = 'none';
  [hd.cone.ax(1,4).XAxis.Visible,hd.cone.ax(1,4).YAxis.Visible] = deal('off');
  for typeindex = 1:numel(datapresent)
    parentpos = hd.cone.ax(typeindex,1).Position;
%     xpos = parentpos(1)+parentpos(3)*-.18;
    xpos = parentpos(1)+parentpos(3)*-.23;
    ypos = parentpos(2)+parentpos(4)*1.30;
%     hd.cone.ob(3,typeindex) = text(xpos, ypos, 'retinal density / mm^{-2}');
    switch typeindex
      case 1
        letterstring = 'a';
      case 4
        letterstring = 'b';
      case 5
        letterstring = 'c';
      otherwise
        letterstring = '';  
    end
    hd.cone.ob(4,typeindex) = text(xpos, ypos, letterstring);
    hd.cone.ob(4,typeindex).FontWeight = 'bold';
    hd.cone.ob(4,typeindex).Interpreter = 'tex';
    hd.cone.ob(4,typeindex).FontSize = hd.cone.ax(typeindex,1).FontSize * 18/14;
    hd.cone.ob(4,typeindex).HorizontalAlignment = 'right';
    hd.cone.ob(4,typeindex).VerticalAlignment = 'top';
  end
  uistack(hd.cone.ax(1,3),'bottom')

  [hd.cone.ax(:,1).TickLength] = deal([0 0]);

  
%   for k = 1:5, hd.cone.cbar(k,1).Position(1) = hd.cone.ax(k,1).Position(1) + .83*hd.cone.ax(k,1).Position(3); end
  
    runtime = clock;
    timestring = ['_', ...
                num2str(runtime(1),'%.4u'), ...
                num2str(runtime(2),'%.2u'), ...
                num2str(runtime(3),'%.2u'), ...
                '_', ...
                num2str(runtime(4),'%.2u'), ...
                num2str(runtime(5),'%.2u'), ...
                num2str(runtime(6),'%02.0f')];
                            
%   savefig(2,['conefigure/',textfilesplit{2},timestring])
%   figure(2),print(['conefigure/',textfilesplit{2},timestring],'-depsc')
% %   figure(2),print(['conefigure/',textfilesplit{2},timestring],'-dtiff')
%   figure(2),print(['conefigure/',textfilesplit{2},timestring],'-dtiffn')
%   figure(2),print(['conefigure/',textfilesplit{2},timestring],'-dpng')
  
end



function [azimuthOut, elevationOut] = diskBoundary(azimuthIn, elevationIn, angleIn, resolutionIn)

  centre.azimuth = azimuthIn;
  centre.elevation = elevationIn;
  border.angle = angleIn / 2; % convert solid angle ("diameter") to "radius"
  border.resolution = resolutionIn;

  % The following equation describes orthodromic (or "great circle") distance as an angle:
  % arccos(sin(centre.elevation)*sin(border.elevation) + ...
  %        cos(centre.elevation)*cos(border.elevation)*cos(border.azimuth-centre.azimuth)) = angle;

%   border.elevation = centre.elevation-border.angle/2 : .1 : centre.elevation+border.angle/2;
  intended  = -90 : border.resolution : 90;
%   intended  = (centre.elevation-border.angle) : border.resolution : (centre.elevation+border.angle);
  % "auxiliary" elevations are explicitly added to guarantee coverage near top/bottom of circle:
  auxiliary = centre.elevation + border.angle*[-1 1];
%   border.elevation = unique([intended,auxiliary]);
  
  limit = cosd([centre.elevation-border.angle centre.elevation+border.angle]);
%   intended = acosd(min(limit) : border.resolution : max(limit));

  %   % XXX Not sure why the following line leads to half of the curve segments not being computed...
%   intended = [-acosd(0:border.resolution:1),acosd(1-border.resolution:-border.resolution:0)];
%   auxiliary = centre.elevation + border.angle*[-1 1];

  border.elevation = unique([intended,auxiliary]);
%   figure(6)
%   plot(intended)

  border.azimuth = centre.azimuth + ...
                   acosd((cosd(border.angle) - sind(centre.elevation)*sind(border.elevation)) ./ ...
                         (cosd(centre.elevation)*cosd(border.elevation)));


  % Eliminate complex elements...
  border.azimuth(imag(border.azimuth)~=0) = NaN; % This statement sometimes fails horribly.
%   border.azimuth(imag(border.azimuth)>1e-5) = NaN; % This statement is silly, but more stable.

  % ...and compute corresponding values for the second half of the boundary.
  border.azimuth = [border.azimuth, 2*centre.azimuth-border.azimuth];
  border.elevation = [border.elevation, border.elevation];
  
  border.azimuth = mod(border.azimuth+180,360)-180;
%   jump = diff(border.azimuth) > 10*border.resolution;
%   border.azimuth(jump)   = NaN;
%   border.elevation(jump) = NaN;
  
  azimuthOut   = border.azimuth;
  elevationOut = border.elevation;

end



% function combinedvalue = vonmisesfisher2param_centre(param,predictor)
% 
%   global globalparam
%   
%   convention = 'CCW';
%   dimension = 3;
%   radius = 1;
%   if size(predictor,2) == 2
%     predictor(:,3) = radius*ones(size(predictor(:,1)));
%   end
%   [x,y,z] = geo2car(predictor(:,1),predictor(:,2),predictor(:,3),convention);
%   direction = [x,y,z];
%   
%   numberofdistributions = 1; % <---
%   value = NaN(size(predictor,1),numberofdistributions);
%   
%   for k = 1:numberofdistributions
%     
%     meanazimuth   = param(2*k-1);
%     meanelevation = param(2*k);
%     
%     offset        = globalparam.offset(k);
%     amplitude     = globalparam.amplitude(k);
%     concentration = globalparam.concentration(k);
%     
%     [meanx,meany,meanz] = geo2car(meanazimuth,meanelevation,radius,convention);
%     meandirection = [meanx,meany,meanz];
% 
%     bessel = concentration / ...
%              (2*pi * (exp(concentration) - exp(-concentration)));
% 
%     normalisation = concentration^(dimension/2 - 1) / ...
%                     ((2*pi)^(dimension/2) * bessel);
%                 
%     value(:,k) = normalisation * exp(concentration*meandirection*direction')';
%     value(:,k) = amplitude * value(:,k) + offset;
%    
%     badness = fitpenalty(meanazimuth,-180,180) + ...
%               fitpenalty(meanelevation,-90,90) + ...
%               fitpenalty(amplitude,0,Inf)      + ...
%               fitpenalty(concentration,0,Inf);
%             
%     value(:,k) = value(:,k) + 1e3*badness;
%     
%   end
%   
%   combinedvalue = sum(value,2);
%   
% end



% function penalty = fitpenalty(value,lowerlimit,upperlimit)
%   
%   factor = 1;
%   
%   if value < lowerlimit
%     
%     penalty = factor*(lowerlimit-value).^1;
%     
%   elseif value > upperlimit
%     
%     penalty = factor*(value-upperlimit).^1;
%     
%   else
%     
%     penalty = 0;
%     
%   end
%   
% end



% function combinedvalue = oneVonMisesFisher(params,predictor)
% %TWOVONMISESFISHER Evaluate two overlapping von Mises-Fisher distributions.
% %
% %   VALUE = VONMISESFISHER computes the value of only one(!) 
% %   spherical (p=3) von Mises-Fisher distributions at the chosen location.
% %
% %   It is formatted to match the input/output structure required by the
% %   built-in non-linear regression function nlinfit.
% 
%   convention = 'CW';
%   dimension = 3;
%   radius = 1;
%   if size(predictor,2) == 2
%     predictor(:,3) = radius*ones(size(predictor(:,1)));
%   end
%   [x,y,z] = geo2car(predictor(:,1),predictor(:,2),predictor(:,3),convention);
%   direction = [x,y,z];
%   
%   numberofdistributions = 1; % <-- This is the only difference to "twoVonMisesFisher.m".
%   value = NaN(size(predictor,1),numberofdistributions);
%   
%   for k = 1:numberofdistributions
%     
%     offset        = params(5*k-4);
%     amplitude     = params(5*k-3);
%     meanazimuth   = params(5*k-2);
%     meanelevation = params(5*k-1);
%     
%         meanelevation = 25; % HACK!
% 
%     
%     [meanx,meany,meanz] = geo2car(meanazimuth,meanelevation,radius,convention);
%     meandirection = [meanx,meany,meanz];
% %     meandirection = geo2car([params((-2:-1)+4*k),radius]);
%     concentration = params(5*k);
%     
%     bessel = concentration / ...
%              (2*pi * (exp(concentration) - exp(-concentration)));
% 
%     normalisation = concentration^(dimension/2 - 1) / ...
%                     ((2*pi)^(dimension/2) * bessel);
%                 
%     value(:,k) = normalisation * exp(concentration*meandirection*direction')';
%     value(:,k) = amplitude * value(:,k) + offset;
%     
%     % HACK to constrain values
%     if abs(meanazimuth) > 180  || ...
%        abs(meanelevation) > 90 || ...
%        amplitude < 0           || ...
%        concentration < 0
%      
%       value(:,k) = 1e10*ones(size(value(:,k)));
%       
%     end
% 
%   end
%   
%   combinedvalue = sum(value,2);
%   
% end