function [hd,area,freq,areastim,freqstim] = figureFreqArea() %#ok<FNDEF>
  
  close all
  
  [freq,~,freqstim] = gainAnalysisFreq;
  [area,~,areastim] = gainAnalysisArea;
  freqoutline = displayStimuli('freq');
  areaoutline = displayStimuli('area');
  freqpattern = displayPattern('freq',freqstim.uniquespatialfreq);
  areapattern = displayPattern('area',freqstim.uniquespatialfreq(4)*ones(1,7));
  
  hd.collage.fg = figure(7);
  hd.collage.fg.Name = 'Sphere manuscript (Fig 7); OKR gain dependence on stimulus size, frequency';
  hd.collage.fg.Color = [1 1 1];
  hd.collage.fg.Position = [0 40 1070 800];
  
  
%   hd.collage.ax(2,1) = axes;
  
  
  % Populate the combined figure from previously generated partial figures...
  hd.collage.ax(1,1) = copyobj(freqpattern.pattern.ax(1), hd.collage.fg);
  hd.collage.ax(1,2) = copyobj(freqoutline.default.ax(1), hd.collage.fg);
  hd.collage.ax(2,1) = copyobj(areapattern.pattern.ax(1), hd.collage.fg);
  hd.collage.ax(2,2) = copyobj(areaoutline.default.ax(1), hd.collage.fg);
  for subpanel = 1:3
    % elements of panel c (location of stimulus centres):
    hd.collage.ax(3,subpanel) = copyobj(area.gain.ax(1), hd.collage.fg);
    % elements of panel d (frequency dependence of gain):
    hd.collage.ax(4,subpanel) = copyobj(freq.gain.ax(1+subpanel), hd.collage.fg);
    % elements of panel e (area dependence of gain):
    hd.collage.ax(5,subpanel) = copyobj(area.gain.ax(1+subpanel), hd.collage.fg);
  end
%   close([1 2 3])
  % ...then adjust the size, position and other properties of those objects.
  
  % Shared properties of all axis objects anywhere:
  [hd.collage.ax(1:2,1:2).Units,hd.collage.ax(3:5,:).Units] = deal('pixels');
  [hd.collage.ax(1:2,1:2).FontSize,hd.collage.ax(3:5,:).FontSize] = deal(15);
  [hd.collage.ax(1:2,1:2).LabelFontSizeMultiplier, ...
   hd.collage.ax(1:2,1:2).TitleFontSizeMultiplier, ...
   hd.collage.ax(3:5,:).LabelFontSizeMultiplier, ...
   hd.collage.ax(3:5,:).TitleFontSizeMultiplier] = deal(1);

  % Positions and sizes within panel a (stimulus frequency):
  hd.collage.ax(1,1).Units = 'pixels';
  hd.collage.ax(1,1).Position = [ 20 600 260 170];
  hd.collage.ax(1,2).Units = 'pixels';
  hd.collage.ax(1,2).Position = [300 600 170 170];

  % Positions and sizes within panel b (stimulus sizes):
  hd.collage.ax(2,1).Units = 'pixels';
  hd.collage.ax(2,1).Position = [540 600 270 170];
  hd.collage.ax(2,2).Units = 'pixels';
  hd.collage.ax(2,2).Position = [820 600 170 170];

  % Positions and sizes within panel c (location of stimulus centres):
  for subpanel = 1:3
    hd.collage.ax(3,subpanel).Position(1) = 30;  % x offset
    hd.collage.ax(3,subpanel).Position(2) = 335 - (subpanel-1)*180;  % y offset
    hd.collage.ax(3,subpanel).Position(3) = 200;  % x width
    hd.collage.ax(3,subpanel).Position(4) = 200;  % y height
  end
  hd.collage.ax(3,2).Position(2) = hd.collage.ax(3,2).Position(2) - 8;  % shift single y offset
  
  % Positions and sizes within panel d (frequency dependence of gain):
  for subpanel = 1:3
    hd.collage.ax(4,subpanel).Position(1) = hd.collage.ax(3,subpanel).Position(1) + ...
                                            hd.collage.ax(3,subpanel).Position(3) + ...
                                            160 + (subpanel-1)*230;  % x offset
    hd.collage.ax(4,subpanel).Position(2) = 350;  % y offset
    hd.collage.ax(4,subpanel).Position(3) = 210;  % x width
    hd.collage.ax(4,subpanel).Position(4) = 142;  % y height
  end

  % Positions and sizes within panel e (area dependence of gain):
  for subpanel = 1:3
    hd.collage.ax(5,subpanel).Position(1) = hd.collage.ax(4,subpanel).Position(1);  % x offset
    hd.collage.ax(5,subpanel).Position(2) = 70;  % y offset
    hd.collage.ax(5,subpanel).Position(3) = hd.collage.ax(4,subpanel).Position(3);  % x width
    hd.collage.ax(5,subpanel).Position(4) = hd.collage.ax(4,subpanel).Position(4);  % y height
  end
  
  % Set camera angles for panels a and b:
  [hd.collage.ax(1:2,2).CameraPosition] = deal([-1 -1.5 0]);
  [hd.collage.ax(1:2,2).CameraViewAngle] = deal(60);
  
  % Set camera angles for panel c:
  hd.collage.ax(3,1).CameraPosition = [0.1 0 20];
  hd.collage.ax(3,2).CameraPosition = [20 0 0];
  hd.collage.ax(3,3).CameraPosition = [-16 -9 6];
  [hd.collage.ax(3,:).CameraViewAngle] = deal(9.5);
  
  % Add azimuth and elevation labels to the plots on panel c:
  axes(hd.collage.ax(3,1))
  labelx = [hd.collage.ax(3,1).Children(1:6).XData] * 1.15;
  labely = [hd.collage.ax(3,1).Children(1:6).YData] * 1.15;
  labelz = [hd.collage.ax(3,1).Children(1:6).ZData] * 0;
  stringlist = areastim.azimuth(7*(1:6)); % unexpected order here, fix by quick hack
  stringlist([1 3 4 6]) = stringlist([3 1 6 4]);
  for location = 1:6
    hold on
    hd.collage.tx(3,1,location) = text(1.1*labelx(location), labely(location), labelz(location), ...
                                       [num2str(stringlist(location),'%+.0f'),'°']);
    hold off
  end
  [hd.collage.tx(3,1,1:3).HorizontalAlignment] = deal('right');
  [hd.collage.tx(3,1,4:6).HorizontalAlignment] = deal('left');
  [hd.collage.tx(3,1,:).VerticalAlignment] = deal('middle');
  axes(hd.collage.ax(3,2))
  labelx = [hd.collage.ax(3,2).Children(1:6).XData] * 0;
  labely = [hd.collage.ax(3,2).Children(1:6).YData] * 1.22;
  labelz = [hd.collage.ax(3,2).Children(1:6).ZData] * 1.22;
  for location = [2 3 5 6]
    hold on
    hd.collage.tx(3,2,location) = text(labelx(location), labely(location), labelz(location), ...
                                       [num2str(areastim.elevation(7*location),'%+.0f'),'°']);
    hold off
  end
  [hd.collage.tx(3,2,2:3).HorizontalAlignment] = deal('right');
  [hd.collage.tx(3,2,5:6).HorizontalAlignment] = deal('left');
  [hd.collage.tx(3,2,[2 3 5 6]).VerticalAlignment] = deal('middle');
  [hd.collage.tx(3,1,:).FontSize, hd.collage.tx(3,2,[2 3 5 6]).FontSize] = ...
    deal(hd.collage.ax(3,1).FontSize);

  
  % Recolour the background spheres on panels a and b, and remove visible axes:
  greyscale = (.8:.01:.99)' * [1 1 1];
  [hd.collage.ax(1:2,2).Colormap] = deal(greyscale);
  for panel = [1 2]
    for subpanel = [1 2]
      [hd.collage.ax(panel,subpanel).XAxis.Visible, ...
       hd.collage.ax(panel,subpanel).YAxis.Visible, ...
       hd.collage.ax(panel,subpanel).ZAxis.Visible] = deal('off');
    end
  end
  
  % Add and format titles to panels a and b:
  for panel = 1:2
    hd.collage.ax(panel,1).Title.String = 'stimulus pattern';
    hd.collage.ax(panel,2).Title.String = 'stimulus crop';
    for subpanel = 1:2
      hd.collage.ax(panel,subpanel).Title.FontWeight = 'normal';
      hd.collage.ax(panel,subpanel).Title.Units = 'pixels';
    end
  hd.collage.ax(panel,1).Title.Position(2) = hd.collage.ax(panel,1).Position(4) - 14 + 14;
  hd.collage.ax(panel,2).Title.Position(2) = hd.collage.ax(panel,2).Position(4) + 2;
  end



  % Reformat the meridians, equators etc. on panel c:
  for subpanel = 1:3
    [hd.collage.ax(3,subpanel).Children(7:9).LineWidth] = deal(1);
  end
  
  % Add and format titles to panel c:
  hd.collage.ax(3,1).Title.String = 'centre azimuth';
  hd.collage.ax(3,2).Title.String = 'centre elevation';
  hd.collage.ax(3,3).Title.String = 'oblique view';
  for subpanel = 1:3
    hd.collage.ax(3,subpanel).Title.FontWeight = 'normal';
    hd.collage.ax(3,subpanel).Title.Units = 'pixels';
  end
  hd.collage.ax(3,1).Title.Position(2) = hd.collage.ax(3,1).Position(4) - 70;
  hd.collage.ax(3,2).Title.Position(2) = hd.collage.ax(3,2).Position(4) - 70;
  hd.collage.ax(3,3).Title.Position(2) = hd.collage.ax(3,3).Position(4) - 45;

  
%   % Remove redundant x axis labels on panels d and e:
%   for panel = [4 5]
%     for subpanel = [1 3]
%       hd.collage.ax(panel,subpanel).XLabel.String = [];
%     end
%   end
  
  % Remove redundant titles on panel e:
  for panel = 5
    for subpanel = [1 2 3]
      hd.collage.ax(panel,subpanel).Title.String = [];
    end
  end
  
  % Reposition X and Y axis labels on panels d and e:
  for panel = [4 5]
    for subpanel = 1:3
      hd.collage.ax(panel,subpanel).XLabel.Units = 'pixels';
    end
    hd.collage.ax(panel,1).YLabel.Units = 'pixels';
  end
  for subpanel = 1:3
    hd.collage.ax(4,subpanel).XLabel.Position(2) = -33;
    hd.collage.ax(5,subpanel).XLabel.Position(2) = -37;
  end
  hd.collage.ax(4,1).YLabel.Position(1) = -58;
  hd.collage.ax(5,1).YLabel.Position(1) = -46;
  
  % Reformat X axis on panel e:
  for subpanel = 1:3
    hd.collage.ax(5,subpanel).XLim(2) = 1.35;
    hd.collage.ax(5,subpanel).XTick = [.1 1];
  end

  
  % Reposition and reformat titles on panel d:
  for subpanel = 1:3
    hd.collage.ax(4,subpanel).Title.Units = 'pixels';
    hd.collage.ax(4,subpanel).Title.Position(2) = hd.collage.ax(4,subpanel).Position(4) + 25;
    hd.collage.ax(4,subpanel).Title.FontWeight = 'bold';
  end
  
  % Overwrite X axis ticks and tick labels on panel d (frequency dependence of gain):
  [hd.collage.ax(4,:).XTick] = deal([.02 .05 .1]);
  [hd.collage.ax(4,:).XTickLabel] = deal({'0.02','0.05','0.10'});

  % Overwrite Y axis ticks and tick labels on panel d (frequency dependence of gain):
  [hd.collage.ax(4,:).YTick] = deal([0 .05 .1]);
  [hd.collage.ax(4,:).YTickLabel] = deal({'0','0.05','0.10'});
  
  % Display panel lettering:
  hd.collage.ax(6,1) = axes;
  hd.collage.ax(6,1).Color = 'none';
  hd.collage.ax(6,1).Position = [0 0 1 1];
  hd.collage.ax(6,1).Units = 'pixels';
  hd.collage.ax(6,1).Visible = 'off';
  hold on
  hd.collage.lt(1) = text(hd.collage.ax(1,1).Position(1) - 0 , ...
                   hd.collage.ax(1,1).Position(2) + hd.collage.ax(1,1).Position(4) + 20, ...
                   'a', 'Units', 'pixels');
  hd.collage.lt(2) = text(hd.collage.ax(2,1).Position(1) - 0, ...
                   hd.collage.ax(2,1).Position(2) + hd.collage.ax(2,1).Position(4) + 20, ...
                   'b', 'Units', 'pixels');
  hd.collage.lt(3) = text(hd.collage.ax(3,1).Position(1) - 10, ...
                   hd.collage.ax(3,1).Position(2) + hd.collage.ax(3,1).Position(4) - 28 + 30, ...
                   'c', 'Units', 'pixels');
  hd.collage.lt(4) = text(hd.collage.ax(4,1).Position(1) - 85, ...
                   hd.collage.ax(4,1).Position(2) + hd.collage.ax(4,1).Position(4) + 30 + 15, ...
                   'd', 'Units', 'pixels');
  hd.collage.lt(5) = text(hd.collage.ax(5,1).Position(1) - 85, ...
                   hd.collage.ax(5,1).Position(2) + hd.collage.ax(5,1).Position(4) + 30, ...
                   'e', 'Units', 'pixels');
  hold off
%   [hd.collage.lt(2:5).FontSize] = deal(15);
%   [hd.collage.lt(2:5).FontWeight] = deal('bold');
  [hd.collage.lt(:).FontSize] = deal(15);
  [hd.collage.lt(:).FontWeight] = deal('bold');
  
end



function [handle,data,stim] = gainAnalysisArea(varargin)
% GAINANALYSIS and all subfunctions were written by Florian Alexander 
% Dehmelt in 2018. This is version 20180924, gainAnalysisArea20180924.m.
% Apologies for the lack of polish (and comments) at this point.
%
% All input arguments are optional:
% 
%   Resolution    - scalar, indicating your choice of resolution for the 
%                   rendering of the spherical surface as part % of the 
%                   figures. Actual data fitting is completely unaffected.
%   CamAzimuth    - scalar, horizontal angle at which sphere is shown
%   CamElevation  - scalar, vertical angle at which sphere is shown
%   Normalisation - options: nothing, or 'mean gain per hemisphere',
%                   or 'max gain per hemisphere', 
%                   or 'ratio of mean gains per hemisphere'
%
% Output arguments are:
%
%   handle  - structure containing handles to all graphics objects
%             (figures, axes, plots, decorations)
%   data    - structure containing pre-processed data, with those fields:
%   .left   - numerical array of means of OKR gain for the left eye
%   .right  - numerical array of means of OKR gain for the right eye
%   .both   - numerical array of mean of OKR gain means across both eyes
%   datafit - structure containing the fits to data, with following fields:  % removed!
%   .centre - 2x2 numerical array of azimuths (1st row) and elevation
%             angles (2nd row) of the two individial von Mises-Fisher 
%             functions fit to data (1st column = 1st function, etc.).
%             These are NOT the peaks of the actual fit (= sum of both
%             functions)! I still have to implement that feature.
%   .value  - 36x1 numerical array containing the value of the actual fit
%             (= sum of both functions) at the locations sampled by data
%   .param  - 1x9 numerical array containing all best-fit parameters
%             (offset, amplitude1, azimuth1, elevation1, concentration1, 
%             amplitude2, azimuth2, elevation2, concentration2)
%
% The recommended default is calling "hd = gainAnalysisArea20181217();".
%
% To specify camera position before analysis, call the function as follows:
%   hd = gainAnalysisArea20181217(200,'CamAzimuth',-30,'CamElevation',15);
%
% To edit camera position post-hoc, adjust and execute the following lines:
%   az = 0; el = 0; % substitute your desired camera coordinates here
%   [hd.gain.ax(1:6).CameraPosition] = ...
%     deal([cosd(az)*cosd(el),sind(az)*cosd(el),sind(el)]/19.0526);
%   figure(21)
%
% To normalise gains on each hemisphere, call the function as follows:
%  hd = gainAnalysisArea20181217('Normalisation','mean gain per hemisphere');
  

  % MANAGE VARIABLE INPUT
  p = inputParser;
  
  valid.resolution    = @(x) isnumeric(x) && numel(x)==1 && x>0;
  valid.camazimuth    = @(x) isnumeric(x) && numel(x)==1;
  valid.camelevation  = @(x) isnumeric(x) && numel(x)==1;
  valid.normalisation = @(x) ischar(x);
                        
  default.resolution    = 200;
  default.camazimuth    = -0;
  default.camelevation  = 45;
  default.normalisation  = '';
  
  addOptional(p, 'Resolution',    default.resolution,    valid.resolution);
  addOptional(p, 'CamAzimuth',    default.camazimuth,    valid.camazimuth);
  addOptional(p, 'CamElevation',  default.camelevation,  valid.camelevation);
  addOptional(p, 'Normalisation', default.normalisation, valid.normalisation);

  parse(p,varargin{:});
  
  resolution    = p.Results.Resolution;
  cam.azimuth   = p.Results.CamAzimuth;
  cam.elevation = p.Results.CamElevation;
  normalisation = p.Results.Normalisation;

  
  % GENERAL SETTINGS
  warning('off','MATLAB:interp1:UsePCHIP')
  % The line above deactivates 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)

  % Identify which stimulus entries correspond to the disk-mask stimuli.
  datarange = 2:43;
  % !!! IMPORTANT !!! These are the Excel row numbers, not stimulus types.
  
%   close all
  
  

  % IDENTIFY DATA FOLDER
  file.folder = '..\..\data\Fig5 size data\Data A Series (A-Series)\';
  

  
  % READ STIMULUS INFORMATION
%   file.folder = 'C:\Users\fdehmelt\Dropbox\ThinkPad Dropbox\MATLAB\201808 SphereAnalysis\Sample Data\D series\'; 
  file.stimulus = 'Stimulus.xlsx';
  [~,stim.tableconvention,~] = xlsread([file.folder,file.stimulus],'C1:C1');
  stim.azimuth = xlsread([file.folder,file.stimulus],['C',num2str(datarange(1)),':C',num2str(datarange(end))]);
  stim.elevation = xlsread([file.folder,file.stimulus],['D',num2str(datarange(1)),':D',num2str(datarange(end))]);
  stim.radius = xlsread([file.folder,file.stimulus],['E',num2str(datarange(1)),':E',num2str(datarange(end))]);

  % Stimulus tables may use CW convention; default code uses CCW.
  switch stim.tableconvention{:}
    case 'azimuth (CCW)'
    case 'azimuth (CW)'
      stim.azimuth = -stim.azimuth;
    otherwise
      error(['Field C1 of Stimulus.xlsx should contain either ' ...
             'string ''azimuth (CCW)'' or string ''azimuth (CW)''.'])
  end
  
  stimcentre = [stim.azimuth, stim.elevation];
  
  
  
  % READ DATA
  % Query user to manually select the main data file. Other files have
  % standard names and are expected to be in the same folder.
  [mfilefolder,~,~] = fileparts(mfilename('fullpath'));
%   [file.result,file.folder] = uigetfile([mfilefolder,'\*.mat']);
  file.result = 'result_20180228_090646.mat';
  load([file.folder,file.result],'result')
  
  % Discard unwanted trials, based on .xlsx documenting manual assessment.
  result = discardTrial(result,file.folder); %#ok<NODEF>
  
%   % 20181217: optionally, normalise gain values by mean gain across same hemisphere
%   switch normalisation
%     case 'mean gain per hemisphere'
%       meangainposazimuth = mean([result(:,:,stim.azimuth>0,:).gain]);
%       meangainnegazimuth = mean([result(:,:,stim.azimuth<0,:).gain]);
%       result(:,:,stim.azimuth>0,:).gain = result(:,:,stim.azimuth>0,:).gain / meangainposazimuth;
%       result(:,:,stim.azimuth<0,:).gain = result(:,:,stim.azimuth<0,:).gain / meangainnegazimuth;
%   end
  
  % Pool gain data in various ways (this is a significant subfunction!).
  gain = poolGainData(result);
  
  % Select the three types of data to be plotted.
  data.left  = gain.mean.FR{1};
  data.right = gain.mean.FR{2};
%   both  = gain.mean.FER;
  data.both  = mean([data.left,data.right],2);
  
  % 20181217: optionally, normalise gain values by mean gain across same hemisphere
  hemisphere1 = datarange(1:end/2) - 1; % The -1 accounts for the Excel title row.
  hemisphere2 = datarange(end/2+1:end) - 1;
  switch normalisation
    case 'mean gain per hemisphere'
      data.left(hemisphere1)  = data.left(hemisphere1)  / mean(data.left(hemisphere1));
      data.left(hemisphere2)  = data.left(hemisphere2)  / mean(data.left(hemisphere2));
      data.right(hemisphere1) = data.right(hemisphere1) / mean(data.right(hemisphere1));
      data.right(hemisphere2) = data.right(hemisphere2) / mean(data.right(hemisphere2));
      data.both(hemisphere1)  = data.both(hemisphere1)  / mean(data.both(hemisphere1));
      data.both(hemisphere2)  = data.both(hemisphere2)  / mean(data.both(hemisphere2));
    case 'max gain per hemisphere'
      data.left(hemisphere1)  = data.left(hemisphere1)  / max(data.left(hemisphere1));
      data.left(hemisphere2)  = data.left(hemisphere2)  / max(data.left(hemisphere2));
      data.right(hemisphere1) = data.right(hemisphere1) / max(data.right(hemisphere1));
      data.right(hemisphere2) = data.right(hemisphere2) / max(data.right(hemisphere2));
      data.both(hemisphere1)  = data.both(hemisphere1)  / max(data.both(hemisphere1));
      data.both(hemisphere2)  = data.both(hemisphere2)  / max(data.both(hemisphere2));
    case 'ratio of mean gains per hemisphere'
      data.left(hemisphere1)  = data.left(hemisphere1) / ...
                                (mean(data.left(hemisphere1)) / mean(data.left(hemisphere2)));
      data.left(hemisphere2)  = data.left(hemisphere2) / ...
                                (mean(data.left(hemisphere2) / mean(data.left(hemisphere1))));
      data.right(hemisphere1) = data.right(hemisphere1) / ...
                                (mean(data.right(hemisphere1)) / mean(data.right(hemisphere2)));
      data.right(hemisphere2) = data.right(hemisphere2) / ...
                                (mean(data.right(hemisphere2)) / mean(data.right(hemisphere1)));
      data.both(hemisphere1)  = data.both(hemisphere1) / ...
                                (mean(data.both(hemisphere1)) / mean(data.both(hemisphere2)));
      data.both(hemisphere2)  = data.both(hemisphere2) / ...
                                (mean(data.both(hemisphere2)) / mean(data.both(hemisphere1)));
  end
  
  
  % CHOOSE SUBSET OF DATA TO ANALYSE, AND PROCEED TO ANALYSE IT
  choicelist = {'left','both','right'};
  for choice = 1:numel(choicelist)
    
    switch choicelist{choice}
      case 'right'
        response = data.right(datarange-1);  % -1 is offset by Excel title row 
      case 'left'
        response = data.left(datarange-1);
      case 'both'
        response = data.both(datarange-1);
      otherwise
        error(['Input argument must be one of the following strings: ' ...
               '''left'', ''right'', or ''both''.'])
    end

    predictor = stimcentre;

    % compute cartesian coordinates of sphere, and its colour data values
    [x,y,z] = sphere(resolution);
%     [x,y,z] = sphere(250);
    [azimuth,elevation,radius] = car2geo(x,y,z,'CCW');
%     predictor = [azimuth(:),elevation(:),radius(:)];
%     datafit.matrix{choice} = reshape(vonmisesfisher9param(datafit.param,predictor),size(radius));

    % set 3D fish parameters
    scale = .1;
    yaw = 180;
    pitch = 0;
    roll = 0;

    % set decorative arc parameters; compute their cartesian coordinates
    equator.base = -1:.01:1;
    equator.x = [equator.base, fliplr(equator.base)];
    equator.y = [sqrt(1-equator.base.^2), -sqrt(1-equator.base.^2)];
    equator.z = 0*ones(1,numel(equator.x));

    meridian.x = equator.x;
    meridian.y = equator.z;
    meridian.z = equator.y;

    aura.x = equator.z;
    aura.y = equator.y;
    aura.z = equator.x;

    [sample.x,sample.y,sample.z] = ...
      geo2car(stimcentre(:,1),stimcentre(:,2),ones(size(stimcentre,1),1),'CCW');
    
    
    [uniquepair,firstoccurrence,~] = unique([stim.azimuth,stim.elevation],'rows');
    % MATLAB reorders unique pairs by azimuth (increasing). This does not
    % usually match the true order in which stimuli were presented! The
    % following two lines recovers this true order, and restores it.
    [~,trueorder] = sort(firstoccurrence);
    uniquepair = uniquepair(trueorder,:);
    stim.uniquecentre.azimuth   = uniquepair(:,1);
    stim.uniquecentre.elevation = uniquepair(:,2);
    
    [uniquecentre.x,uniquecentre.y,uniquecentre.z] = ...
      geo2car(uniquepair(:,1),uniquepair(:,2),ones(size(uniquepair,1),1),'CCW');

    stim.uniqueradius = unique(stim.radius);
    stim.uniquesteradian = stim.uniqueradius*(2*pi)/360;
    stim.uniquearea = stim.uniquesteradian/(2*pi);
    
    fontsize = 14;
    
    if choice == 1
      
      panel = 1;
%       hd.location.fg = figure();
%       hd.location.fg.Position = [50 300 1200 300];
%       hd.location.fg.Color = [1 1 1];

      hd.gain.ax(1) = axes;
      hd.gain.ax(1).Position = [.025 .14 .20 .8];
      hold on
      hd.gain.ob(panel,1) = surf(x,y,z);
      hd.gain.ob(panel,1).FaceAlpha = 0*.25;
      hd.gain.ob(panel,1).FaceColor = 1*[1 1 1];
      hd.gain.ob(panel,1).EdgeColor = 'none';

      hd.gain.ob(panel,2:4) = plot3dFish(scale,yaw,pitch,roll);
      hd.gain.ob(panel,5) = plot3(equator.x,equator.y,equator.z,'k','LineWidth',2);
      hd.gain.ob(panel,6) = plot3(meridian.x,meridian.y,meridian.z,'--k');
      hd.gain.ob(panel,7) = plot3(aura.x,aura.y,aura.z,':k');
      
%       colour.centre = [1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1; 0 1 1];
      colour.centre = [1 .447 .741; .85 .325 .098; .929 .694 .125; ...
                       .494 .184 .556; .466 .6740 .188; .301 .745 .933];
      colour.purple = cbrewer('seq','Purples',100);
      colour.green = cbrewer('seq','Greens',100);
      colour.centre = [colour.purple(round(end*.8),:); ...
                       colour.purple(round(end*.3),:); ...
                       colour.purple(round(end*.5),:); ...
                       colour.green(round(end*.8),:); ...
                       colour.green(round(end*.3),:); ...
                       colour.green(round(end*.5),:)];
      
      for origin = 1:numel(uniquecentre.x)
        hd.gain.ob(panel,7+origin) = scatter3(1.03*uniquecentre.x(origin), ...
                                              1.03*uniquecentre.y(origin), ...
                                              1.03*uniquecentre.z(origin), ...
                                              120,'Marker','o');
        hd.gain.ob(panel,7+origin).MarkerFaceColor = colour.centre(origin,:);
        hd.gain.ob(panel,7+origin).MarkerEdgeColor = .6*[1 1 1];
      end
      hold off
%   %     hd.gain.ob(panel,8).Color = [1 1 1];
%       hd.gain.ob(panel,1).FaceAlpha = .5;
%       hd.gain.ob(panel,1).FaceColor = hd.gain.ob(panel,1).EdgeColor;
      axis square
      
      hd.gain.ax(panel).XAxis.Visible = 'off';
      hd.gain.ax(panel).YAxis.Visible = 'off';
      hd.gain.ax(panel).ZAxis.Visible = 'off';
%       hd.gain.ax(panel).Title.String = 'stimulus centres';
      
      hd.gain.ax(panel).CameraPosition = [11.5866 -8.9228 9.2809];
      hd.gain.ax(panel).CameraViewAngle = 9.25;
      
      hd.gain.ax(4) = axes;
      hd.gain.ax(4).Position = [0 0 1 1];
      hd.gain.ax(4).Color = 'none';
      hd.gain.ax(4).Visible = 'off';
      
      hd.gain.ob(4,1) = text(.02,.93,'a','Units','normalized');
      hd.gain.ob(4,2) = text(.25,.93,'b','Units','normalized');
      [hd.gain.ob(4,1:2).FontWeight] = deal('bold');
      [hd.gain.ob(4,1:2).FontSize] = deal(14);

      
    end
    input.letter = 'e';

    panel = choice+1;
    croppeddata(:,1) = response(1:7);
    croppeddata(:,2) = response(8:14);
    croppeddata(:,3) = response(15:21);
    croppeddata(:,4) = response(22:28);
    croppeddata(:,5) = response(29:35);
    croppeddata(:,6) = response(36:42);
    hd.gain.ax(panel) = axes;
    hd.gain.ax(panel).Position = [.33+(panel-2)*.22 .26 .20 .55];
    hold on
    for origin = 1:size(croppeddata,2)
      hd.gain.ob(panel,origin) = plot(stim.uniquearea,croppeddata(:,origin));
      hd.gain.ob(panel,origin).Color = colour.centre(origin,:);
      hd.gain.ob(panel,origin).LineWidth = 2;
    end
    hold off
    
    hd.gain.ax(panel).TitleFontWeight = 'normal';
    hd.gain.ax(panel).TitleFontSizeMultiplier = 1;

    hd.gain.ax(panel).XLabel.String = 'stimulus size';
    if choice == 1
      hd.gain.ax(panel).YLabel.String = 'OKR gain';
      hd.gain.ax(panel).YLabel.Units = 'normalized';
      hd.gain.ax(2).YLabel.Position = [-.2 .5 0];
    else
      hd.gain.ax(panel).YAxis.Visible = 'off';
    end
    
    switch choicelist{choice}
      case 'left'
        hd.gain.ax(panel).Title.String = 'left eye';
      case 'right'
        hd.gain.ax(panel).Title.String = 'right eye';
      case 'both'
%         hd.gain.ax(panel).Title.String = 'mean of both eyes';
        hd.gain.ax(panel).Title.String = 'mean';
    end
    
    hd.gain.ax(panel).FontSize = fontsize;
    
    if panel > 1
      hd.gain.ax(panel).XScale = 'log';
      hd.gain.ax(panel).YScale = 'lin';
      hd.gain.ax(panel).XLim = [.9*min(hd.gain.ob(panel,1).XData), ...
                                max(hd.gain.ob(panel,1).XData)];
      switch normalisation
        case 'max gain per hemisphere'
          hd.gain.ax(panel).YLim = [0, 1];
        case 'mean gain per hemisphere'
          hd.gain.ax(panel).YLim = [0, max([data.left;data.right;data.both])];
        case 'ratio of mean gains per hemisphere'
          hd.gain.ax(panel).YLim = [0, max([data.left;data.right;data.both])];
        otherwise
          hd.gain.ax(panel).YLim = [0, .5];
      end
      hd.gain.ax(panel).TickLength = [.03 .075];
      hd.gain.ax(panel).YGrid = 'on';
      hd.gain.ax(panel).Title.Units = 'normalized';
      hd.gain.ax(panel).Title.Position = [.5 1.1 0];
    end
    
  end
  
    cam.radius = 19.0526; % MATLAB automatically reverts to this values while rotating, anyway. Might as well choose it as default.
    [cam.x,cam.y,cam.z] = geo2car(cam.azimuth,cam.elevation,cam.radius,'CCW');
    hd.gain.ax(1).CameraPosition = [cam.x cam.y cam.z];
    hd.gain.ax(1).CameraTarget = [0 0 0];
    hd.gain.ax(1).CameraUpVector = [0 0 1];
%     hd.gain.ax(panel).CameraViewAngle = 10.134;
    hd.gain.ax(1).CameraViewAngle = 9;
  
  
  % Collect output arguments, unless already done (e.g., "datafit").
%   data.left = left;
%   data.right = right;
%   data.both = both;
  handle = hd;
  
end






function [handle,data,stim] = gainAnalysisFreq(varargin)
% GAINANALYSIS and all subfunctions were written by Florian Alexander 
% Dehmelt in 2018. This is version 20180927, gainAnalysisFrequency20180924.m.
% Apologies for the lack of polish (and comments) at this point.
%
% All input arguments are optional:
% 
%   Resolution    - scalar, indicating your choice of resolution for the 
%                   rendering of the spherical surface as part % of the 
%                   figures. Actual data fitting is completely unaffected.
%   CamAzimuth    - scalar, horizontal angle at which sphere is shown
%   CamElevation  - scalar, vertical angle at which sphere is shown
%   Normalisation - options: nothing, or 'mean gain per hemisphere',
%                   or 'max gain per hemisphere', 
%                   or 'ratio of mean gains per hemisphere'
%
% Output arguments are:
%
%   handle  - structure containing handles to all graphics objects
%             (figures, axes, plots, decorations)
%   data    - structure containing pre-processed data, with those fields:
%   .left   - numerical array of means of OKR gain for the left eye
%   .right  - numerical array of means of OKR gain for the right eye
%   .both   - numerical array of mean of OKR gain means across both eyes
%   datafit - structure containing the fits to data, with following fields:  % removed!
%   .centre - 2x2 numerical array of azimuths (1st row) and elevation
%             angles (2nd row) of the two individial von Mises-Fisher 
%             functions fit to data (1st column = 1st function, etc.).
%             These are NOT the peaks of the actual fit (= sum of both
%             functions)! I still have to implement that feature.
%   .value  - 36x1 numerical array containing the value of the actual fit
%             (= sum of both functions) at the locations sampled by data
%   .param  - 1x9 numerical array containing all best-fit parameters
%             (offset, amplitude1, azimuth1, elevation1, concentration1, 
%             amplitude2, azimuth2, elevation2, concentration2)
%
% The recommended default is calling "hd = gainAnalysis20180911(200);".
%
% To specify camera position before analysis, call the function as follows:
%   hd = gainAnalysis20180911(200,'CamAzimuth',-30,'CamElevation',15);
%
% To edit camera position post-hoc, adjust and execute the following lines:
%   az = 0; el = 0; % substitute your desired camera coordinates here
%   [hd.gain.ax(1:6).CameraPosition] = ...
%     deal([cosd(az)*cosd(el),sind(az)*cosd(el),sind(el)]/19.0526);
%   figure(21)
  

  % MANAGE VARIABLE INPUT
  p = inputParser;
  
  valid.resolution    = @(x) isnumeric(x) && numel(x)==1 && x>0;
  valid.camazimuth    = @(x) isnumeric(x) && numel(x)==1;
  valid.camelevation  = @(x) isnumeric(x) && numel(x)==1;
  valid.normalisation = @(x) ischar(x);
                        
  default.resolution    = 200;
  default.camazimuth    = -0;
  default.camelevation  = 45;
  default.normalisation  = '';
  
  addOptional(p, 'Resolution',    default.resolution,    valid.resolution);
  addOptional(p, 'CamAzimuth',    default.camazimuth,    valid.camazimuth);
  addOptional(p, 'CamElevation',  default.camelevation,  valid.camelevation);
  addOptional(p, 'Normalisation', default.normalisation, valid.normalisation);

  parse(p,varargin{:});
  
  resolution    = p.Results.Resolution;
  cam.azimuth   = p.Results.CamAzimuth;
  cam.elevation = p.Results.CamElevation;
  normalisation = p.Results.Normalisation;
  
  
  % GENERAL SETTINGS
  warning('off','MATLAB:interp1:UsePCHIP')
  % The line above deactivates 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)

  % Identify which stimulus entries correspond to the disk-mask stimuli.
  datarange = 2:43;
  % !!! IMPORTANT !!! These are the Excel row numbers, not stimulus types.
  
%   close all
  
  
  % READ DATA
  % Query user to manually select the main data file. Other files have
  % standard names and are expected to be in the same folder.
  [mfilefolder,~,~] = fileparts(mfilename('fullpath'));
%   [file.result,file.folder] = uigetfile([mfilefolder,'\*.mat']);
  file.result = 'result_20180305_124441.mat';
  file.folder = '..\..\data\Fig5 frequency data\Data V Series (V-Serie)\';
  load([file.folder,file.result],'result')
  
  % Discard unwanted trials, based on .xlsx documenting manual assessment.
  result = discardTrial(result,file.folder); %#ok<NODEF>
  result = zeroTrial(result,file.folder); %#ok<NODEF>
  
  % Pool gain data in various ways (this is a significant subfunction!).
  gain = poolGainData(result);
  
  % Select the three types of data to be plotted.
  data.left  = gain.mean.FR{1};
  data.right = gain.mean.FR{2};
%   both  = gain.mean.FER;
  data.both  = mean([data.left,data.right],2);
  
  % 20181217: optionally, normalise gain values by mean gain across same hemisphere
  hemisphere1 = datarange(1:end/2) - 1; % The -1 accounts for the Excel title row.
  hemisphere2 = datarange(end/2+1:end) - 1;
  switch normalisation
    case 'mean gain per hemisphere'
      data.left(hemisphere1)  = data.left(hemisphere1)  / mean(data.left(hemisphere1));
      data.left(hemisphere2)  = data.left(hemisphere2)  / mean(data.left(hemisphere2));
      data.right(hemisphere1) = data.right(hemisphere1) / mean(data.right(hemisphere1));
      data.right(hemisphere2) = data.right(hemisphere2) / mean(data.right(hemisphere2));
      data.both(hemisphere1)  = data.both(hemisphere1)  / mean(data.both(hemisphere1));
      data.both(hemisphere2)  = data.both(hemisphere2)  / mean(data.both(hemisphere2));
    case 'max gain per hemisphere'
      data.left(hemisphere1)  = data.left(hemisphere1)  / max(data.left(hemisphere1));
      data.left(hemisphere2)  = data.left(hemisphere2)  / max(data.left(hemisphere2));
      data.right(hemisphere1) = data.right(hemisphere1) / max(data.right(hemisphere1));
      data.right(hemisphere2) = data.right(hemisphere2) / max(data.right(hemisphere2));
      data.both(hemisphere1)  = data.both(hemisphere1)  / max(data.both(hemisphere1));
      data.both(hemisphere2)  = data.both(hemisphere2)  / max(data.both(hemisphere2));
    case 'ratio of mean gains per hemisphere'
      data.left(hemisphere1)  = data.left(hemisphere1) / ...
                                (mean(data.left(hemisphere1)) / mean(data.left(hemisphere2)));
      data.left(hemisphere2)  = data.left(hemisphere2) / ...
                                (mean(data.left(hemisphere2) / mean(data.left(hemisphere1))));
      data.right(hemisphere1) = data.right(hemisphere1) / ...
                                (mean(data.right(hemisphere1)) / mean(data.right(hemisphere2)));
      data.right(hemisphere2) = data.right(hemisphere2) / ...
                                (mean(data.right(hemisphere2)) / mean(data.right(hemisphere1)));
      data.both(hemisphere1)  = data.both(hemisphere1) / ...
                                (mean(data.both(hemisphere1)) / mean(data.both(hemisphere2)));
      data.both(hemisphere2)  = data.both(hemisphere2) / ...
                                (mean(data.both(hemisphere2)) / mean(data.both(hemisphere1)));
  end

  
  % READ STIMULUS INFORMATION
  file.stimulus = 'Stimulus.xlsx';
  [~,stim.tableconvention,~] = xlsread([file.folder,file.stimulus],'C1:C1');
  stim.azimuth = xlsread([file.folder,file.stimulus],['C',num2str(datarange(1)),':C',num2str(datarange(end))]);
  stim.elevation = xlsread([file.folder,file.stimulus],['D',num2str(datarange(1)),':D',num2str(datarange(end))]);
  stim.solidangle = xlsread([file.folder,file.stimulus],['E',num2str(datarange(1)),':E',num2str(datarange(end))]);
  stim.spatialfreq = xlsread([file.folder,file.stimulus],['F',num2str(datarange(1)),':F',num2str(datarange(end))]);
  stim.tempfreq = xlsread([file.folder,file.stimulus],['G',num2str(datarange(1)),':G',num2str(datarange(end))]);

  % Stimulus tables may use CW convention; default code uses CCW.
  switch stim.tableconvention{:}
    case 'azimuth (CCW)'
    case 'azimuth (CW)'
      stim.azimuth = -stim.azimuth;
    otherwise
      error(['Field C1 of Stimulus.xlsx should contain either ' ...
             'string ''azimuth (CCW)'' or string ''azimuth (CW)''.'])
  end
  
  stimcentre = [stim.azimuth, stim.elevation];

  
  % CHOOSE SUBSET OF DATA TO ANALYSE, AND PROCEED TO ANALYSE IT
  choicelist = {'left','both','right'};
  for choice = 1:numel(choicelist)
    
    switch choicelist{choice}
      case 'right'
        response = data.right(datarange-1);  % -1 is offset by Excel title row 
      case 'left'
        response = data.left(datarange-1);
      case 'both'
        response = data.both(datarange-1);
      otherwise
        error(['Input argument must be one of the following strings: ' ...
               '''left'', ''right'', or ''both''.'])
    end

    predictor = stimcentre;

    % compute cartesian coordinates of sphere, and its colour data values
    [x,y,z] = sphere(resolution);
%     [x,y,z] = sphere(250);
    [azimuth,elevation,radius] = car2geo(x,y,z,'CCW');
%     predictor = [azimuth(:),elevation(:),radius(:)];
%     datafit.matrix{choice} = reshape(vonmisesfisher9param(datafit.param,predictor),size(radius));

    % set 3D fish parameters
    scale = .1;
    yaw = 180;
    pitch = 0;
    roll = 0;

    % set decorative arc parameters; compute their cartesian coordinates
    equator.base = -1:.01:1;
    equator.x = [equator.base, fliplr(equator.base)];
    equator.y = [sqrt(1-equator.base.^2), -sqrt(1-equator.base.^2)];
    equator.z = 0*ones(1,numel(equator.x));

    meridian.x = equator.x;
    meridian.y = equator.z;
    meridian.z = equator.y;

    aura.x = equator.z;
    aura.y = equator.y;
    aura.z = equator.x;

    [sample.x,sample.y,sample.z] = ...
      geo2car(stimcentre(:,1),stimcentre(:,2),ones(size(stimcentre,1),1),'CCW');
    
    
    [uniquepair,firstoccurrence,~] = unique([stim.azimuth,stim.elevation],'rows');
    % MATLAB reorders unique pairs by azimuth (increasing). This does not
    % usually match the true order in which stimuli were presented! The
    % following two lines recovers this true order, and restores it.
    [~,trueorder] = sort(firstoccurrence);
    uniquepair = uniquepair(trueorder,:);
    stim.uniquecentre.azimuth   = uniquepair(:,1);
    stim.uniquecentre.elevation = uniquepair(:,2);
    
    [uniquecentre.x,uniquecentre.y,uniquecentre.z] = ...
      geo2car(uniquepair(:,1),uniquepair(:,2),ones(size(uniquepair,1),1),'CCW');

    stim.uniquespatialfreq = unique(stim.spatialfreq);
    stim.uniquetempfreq = unique(stim.tempfreq);
    
    stim.uniquesolidangle = unique(stim.solidangle);
    stim.uniquesteradian = stim.uniquesolidangle*(2*pi)/360;
    stim.uniquearea = stim.uniquesteradian/(2*pi);
    
    fontsize = 14;
    
    if choice == 1
      
      panel = 1;
      hd.location.fg = figure();
      hd.location.fg.Position = [50 300 1200 300];
      hd.location.fg.Color = [1 1 1];

      hd.gain.ax(1) = axes;
      hd.gain.ax(1).Position = [.025 .14 .20 .8];
      hold on
      hd.gain.ob(panel,1) = surf(x,y,z);
      hd.gain.ob(panel,1).FaceAlpha = 0*.25;
      hd.gain.ob(panel,1).FaceColor = 1*[1 1 1];
      hd.gain.ob(panel,1).EdgeColor = 'none';

      hd.gain.ob(panel,2:4) = plot3dFish(scale,yaw,pitch,roll);
      hd.gain.ob(panel,5) = plot3(equator.x,equator.y,equator.z,'k','LineWidth',2);
      hd.gain.ob(panel,6) = plot3(meridian.x,meridian.y,meridian.z,'--k');
      hd.gain.ob(panel,7) = plot3(aura.x,aura.y,aura.z,':k');
      
%       colour.centre = [1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1; 0 1 1];
      colour.centre = [1 .447 .741; .85 .325 .098; .929 .694 .125; ...
                       .494 .184 .556; .466 .6740 .188; .301 .745 .933];
      colour.purple = cbrewer('seq','Purples',100);
      colour.green = cbrewer('seq','Greens',100);
      colour.centre = [colour.purple(round(end*.8),:); ...
                       colour.purple(round(end*.3),:); ...
                       colour.purple(round(end*.5),:); ...
                       colour.green(round(end*.8),:); ...
                       colour.green(round(end*.3),:); ...
                       colour.green(round(end*.5),:)];
      
      for origin = 1:numel(uniquecentre.x)
        hd.gain.ob(panel,7+origin) = scatter3(1.03*uniquecentre.x(origin), ...
                                              1.03*uniquecentre.y(origin), ...
                                              1.03*uniquecentre.z(origin), ...
                                              120,'Marker','o');
        hd.gain.ob(panel,7+origin).MarkerFaceColor = colour.centre(origin,:);
        hd.gain.ob(panel,7+origin).MarkerEdgeColor = .6*[1 1 1];
      end
      hold off
%   %     hd.gain.ob(panel,8).Color = [1 1 1];
%       hd.gain.ob(panel,1).FaceAlpha = .5;
%       hd.gain.ob(panel,1).FaceColor = hd.gain.ob(panel,1).EdgeColor;
      axis square
      
      hd.gain.ax(panel).XAxis.Visible = 'off';
      hd.gain.ax(panel).YAxis.Visible = 'off';
      hd.gain.ax(panel).ZAxis.Visible = 'off';
%       hd.gain.ax(panel).Title.String = 'stimulus centres';
      
      hd.gain.ax(panel).CameraPosition = [11.5866 -8.9228 9.2809];
      hd.gain.ax(panel).CameraViewAngle = 9.25;
      
      hd.gain.ax(4) = axes;
      hd.gain.ax(4).Position = [0 0 1 1];
      hd.gain.ax(4).Color = 'none';
      hd.gain.ax(4).Visible = 'off';
      
      hd.gain.ob(4,1) = text(.02,.93,'c','Units','normalized');
      hd.gain.ob(4,2) = text(.25,.93,'d','Units','normalized');
      [hd.gain.ob(4,1:2).FontWeight] = deal('bold');
      [hd.gain.ob(4,1:2).FontSize] = deal(14);

      
    end
    input.letter = 'e';

    panel = choice+1;
    croppeddata(:,1) = response(1:7);
    croppeddata(:,2) = response(8:14);
    croppeddata(:,3) = response(15:21);
    croppeddata(:,4) = response(22:28);
    croppeddata(:,5) = response(29:35);
    croppeddata(:,6) = response(36:42);
    hd.gain.ax(panel) = axes;
    hd.gain.ax(panel).Position = [.33+(panel-2)*.22 .26 .20 .55];
    hold on
    for origin = 1:size(croppeddata,2)
%       hd.gain.ob(panel,origin) = plot(stim.uniquetempfreq,croppeddata(:,origin));
      hd.gain.ob(panel,origin) = plot(stim.uniquespatialfreq,croppeddata(:,origin));
      hd.gain.ob(panel,origin).Color = colour.centre(origin,:);
      hd.gain.ob(panel,origin).LineWidth = 2;
    end
    hold off
    
    hd.gain.ax(panel).TitleFontWeight = 'normal';
    hd.gain.ax(panel).TitleFontSizeMultiplier = 1;

%     hd.gain.ax(panel).XLabel.String = 'temporal frequency';
    hd.gain.ax(panel).XLabel.String = 'spatial frequency';
    if choice == 1
      hd.gain.ax(panel).YLabel.String = 'OKR gain';
      hd.gain.ax(panel).YLabel.Units = 'normalized';
      hd.gain.ax(2).YLabel.Position = [-.2 .5 0];
    else
      hd.gain.ax(panel).YAxis.Visible = 'off';
    end
    
    switch choicelist{choice}
      case 'left'
        hd.gain.ax(panel).Title.String = 'left eye';
      case 'right'
        hd.gain.ax(panel).Title.String = 'right eye';
      case 'both'
%         hd.gain.ax(panel).Title.String = 'mean of both eyes';
        hd.gain.ax(panel).Title.String = 'mean';
    end
    
    hd.gain.ax(panel).FontSize = fontsize;
    
    if panel > 1
      hd.gain.ax(panel).XScale = 'log';
      hd.gain.ax(panel).YScale = 'lin';
      hd.gain.ax(panel).XLim = [.9*min(hd.gain.ob(panel,1).XData), ...
                                max(hd.gain.ob(panel,1).XData)];
      switch normalisation
        case 'max gain per hemisphere'
          hd.gain.ax(panel).YLim = [0, 1];
        case 'mean gain per hemisphere'
          hd.gain.ax(panel).YLim = [0, max([data.left;data.right;data.both])];
        case 'ratio of mean gains per hemisphere'
          hd.gain.ax(panel).YLim = [0, max([data.left;data.right;data.both])];
        otherwise
          hd.gain.ax(panel).YLim = [0, max([data.left;data.right;data.both])];
      end
      hd.gain.ax(panel).TickLength = [.03 .075];
      hd.gain.ax(panel).YGrid = 'on';
      hd.gain.ax(panel).Title.Units = 'normalized';
      hd.gain.ax(panel).Title.Position = [.5 1.1 0];
    end
  
  end
  
    cam.radius = 19.0526; % MATLAB automatically reverts to this values while rotating, anyway. Might as well choose it as default.
    [cam.x,cam.y,cam.z] = geo2car(cam.azimuth,cam.elevation,cam.radius,'CCW');
    hd.gain.ax(1).CameraPosition = [cam.x cam.y cam.z];
    hd.gain.ax(1).CameraTarget = [0 0 0];
    hd.gain.ax(1).CameraUpVector = [0 0 1];
%     hd.gain.ax(panel).CameraViewAngle = 10.134;
    hd.gain.ax(1).CameraViewAngle = 9;

    
  % Collect output arguments, unless already done (e.g., "datafit").
%   data.left = left;
%   data.right = right;
%   data.both = both;
  handle = hd;
  
end



function gain = poolGainData(result)
%POOLGAINDATA converts a raw result structure into plottable data by
%pooling OKR gain data in various ways, and storing those results in a new
%structure called "gain".

  % Pool data in various ways (e.g., .FR = pooled across all [F]ish and
  % all [R]epetitions, but with a separate array for each [E]ye).
  for stimtype = 1:size(result,3)

    gainpool(stimtype).FER = [result(:,:,stimtype,:).gain];

    for eye = 1:size(result,2)
      gainpool(stimtype).FR{eye} = [result(:,eye,stimtype,:).gain];
    end

    for fish = 1:size(result,1)
      gainpool(stimtype).ER{fish} = [result(fish,:,stimtype,:).gain];
    end

    for repetition = 1:size(result,4)
      gainpool(stimtype).FE{repetition} = [result(:,:,stimtype,repetition).gain];
    end

    for fish = 1:size(result,1)
      for eye = 1:size(result,2)
        gainpool(stimtype).R{fish,eye} = [result(fish,eye,stimtype,:).gain];
      end
    end
    
  end

  % Make sure that in following section, if a stimulus type is not present,
  % the array contain a NaN at that position, instead of a zero.
  nantemplate = NaN(numel(gainpool),1);

  gain.mean  .FER = nantemplate;
  gain.median.FER = nantemplate;
  gain.std   .FER = nantemplate;
  gain.sem   .FER = nantemplate;

  for eye = 1:size(result,2)
    gain.mean  .FR{eye} = nantemplate;
    gain.median.FR{eye} = nantemplate;
    gain.std   .FR{eye} = nantemplate;
    gain.sem   .FR{eye} = nantemplate;
  end

  for fish = 1:size(result,1)
    gain.mean  .ER{fish} = nantemplate;
    gain.median.ER{fish} = nantemplate;
    gain.std   .ER{fish} = nantemplate;
    gain.sem   .ER{fish} = nantemplate;
  end

  for fish = 1:size(result,1)
    for eye = 1:size(result,2)
      gain.mean  .R{fish,eye} = nantemplate;
      gain.median.R{fish,eye} = nantemplate;
      gain.std   .R{fish,eye} = nantemplate;
      gain.sem   .R{fish,eye} = nantemplate;
    end
  end

  % mean, median, stdev, sem across all pooled data, one per stimulus type
  for stimtype = 1:numel(gainpool)

    selection = gainpool(stimtype).FER;
    gain.mean  .FER(stimtype) = mean(selection);
    gain.median.FER(stimtype) = median(selection);
    gain.std   .FER(stimtype) = std(selection);
    gain.sem   .FER(stimtype) = std(selection)/sqrt(numel(selection));

    for eye = 1:size(result,2)
      selection = gainpool(stimtype).FR{eye};
      gain.mean  .FR{eye}(stimtype) = mean(selection);
      gain.median.FR{eye}(stimtype) = median(selection);
      gain.std   .FR{eye}(stimtype) = std(selection);
      gain.sem   .FR{eye}(stimtype) = std(selection)/sqrt(numel(selection));
    end

    for fish = 1:size(result,1)
      selection = gainpool(stimtype).ER{fish};
      gain.mean  .ER{fish}(stimtype) = mean(selection);
      gain.median.ER{fish}(stimtype) = median(selection);
      gain.std   .ER{fish}(stimtype) = std(selection);
      gain.sem   .ER{fish}(stimtype) = std(selection)/sqrt(numel(selection));
    end

    for fish = 1:size(result,1)
      for eye = 1:size(result,2)
        selection = gainpool(stimtype).R{fish,eye};
        gain.mean  .R{fish,eye}(stimtype) = mean(selection);
        gain.median.R{fish,eye}(stimtype) = median(selection);
        gain.std   .R{fish,eye}(stimtype) = std(selection);
        gain.sem   .R{fish,eye}(stimtype) = std(selection)/sqrt(numel(selection));
      end
    end
  end

end
  
  

function hd = fishbowlFigure(hd,cam,choicelist)
%FISHBOWLFIGURE configures an existing figure of fishbowl plots, adds decorations
%
%   HD = fishbowlFigure sets the parameters of an existing fishbowl figure,
%   including axis limits, labels, colours, axes object sizes, etc., and
%   adds additional decorations such as colorbars.
  
  cminbottomrow = min([min(min(hd.gain.ob(4,5).CData)), ...
                       min(min(hd.gain.ob(5,5).CData)), ...
                       min(min(hd.gain.ob(6,5).CData))]);
  cmaxbottomrow = max([max(max(hd.gain.ob(4,5).CData)), ...
                       max(max(hd.gain.ob(5,5).CData)), ...
                       max(max(hd.gain.ob(6,5).CData))]);
  climbottomrow = max(abs([cminbottomrow cmaxbottomrow]))*[-1 1];
  
  cmintoprow = min([min(min(hd.gain.ob(1,1).CData)), ...
                    min(min(hd.gain.ob(2,1).CData)), ...
                    min(min(hd.gain.ob(3,1).CData)), ...
                    min(min(hd.gain.ob(1,8).CData)), ...
                    min(min(hd.gain.ob(2,8).CData)), ...
                    min(min(hd.gain.ob(3,8).CData))]);
  cmaxtoprow = max([max(max(hd.gain.ob(1,1).CData)), ...
                    max(max(hd.gain.ob(2,1).CData)), ...
                    max(max(hd.gain.ob(3,1).CData)), ...
                    max(max(hd.gain.ob(1,8).CData)), ...
                    max(max(hd.gain.ob(2,8).CData)), ...
                    max(max(hd.gain.ob(3,8).CData))]);
  climtoprow = max(abs([cmintoprow cmaxtoprow]))*[-1 1];
  
  for panel = 1:3
    
    switch choicelist{panel}
      case 'left'
        string.title = 'left eye gain g_L';
      case 'right'
        string.title = 'right eye gain g_R';
      case 'both'
        string.title = 'mean = (g_L+ g_R)/2';
    end
       
    hd.gain.ax(panel).Title.String = string.title;
%     title(string.title,'Interpreter','LaTeX')
    hd.gain.ob(panel,8).MarkerFaceColor = hd.gain.ob(panel,8).MarkerEdgeColor;
    hd.gain.ob(panel,8).MarkerEdgeColor = [0 0 0];
    hd.gain.ax(panel).CLim = [0 cmaxtoprow];
%     hd.gain.ax(panel).CLim = [0 max(max(datafit.matrix{choice}))];
    hd.gain.ob(panel,1).FaceAlpha = .85;
    hd.gain.ob(panel,1).EdgeColor = 'none';
%     hd.gain.ob(panel,8).Shading = 'flat';
    shading flat
    
  end
  for panel = 4:6
   
    hd.gain.ax(panel).CLim = climbottomrow;
    hd.gain.ob(panel,5).MarkerFaceColor = hd.gain.ob(panel,5).MarkerEdgeColor;
    hd.gain.ob(panel,1).FaceColor = [1 1 1];
    hd.gain.ob(panel,1).FaceAlpha = .85;
    hd.gain.ob(panel,1).EdgeColor = 'none';

  end
  for panel = 1:6
    
    axis square
    xlabel('x')
    ylabel('y')
    zlabel('z')            
    hd.gain.ax(panel).XLim = [-1.1 1.1];
    hd.gain.ax(panel).YLim = [-1.1 1.1];
    hd.gain.ax(panel).ZLim = [-1.1 1.1];
%     hd.gain.ax(panel).CameraPosition = [12.2934, -9.0967,8.1315];
    cam.radius = 19.0526; % MATLAB automatically reverts to this values while rotating, anyway. Might as well choose it as default.
    [cam.x,cam.y,cam.z] = geo2car(cam.azimuth,cam.elevation,cam.radius,'CCW');
    hd.gain.ax(panel).CameraPosition = [cam.x cam.y cam.z];
    hd.gain.ax(panel).CameraTarget = [0 0 0];
    hd.gain.ax(panel).CameraUpVector = [0 0 1];
%     hd.gain.ax(panel).CameraViewAngle = 10.134;
    hd.gain.ax(panel).CameraViewAngle = 9;

  end
  
  hd.gain.ax(7) = axes;
  hd.gain.ax(7).Position = hd.gain.ax(1).Position - [.13 -.1 .2 .2];
  hd.gain.ax(7).CLim = hd.gain.ax(1).CLim;
  hd.gain.ob(7,1) = colorbar;
  hd.gain.ob(7,1).Label.String = 'OKR gain';
  hd.gain.ob(7,1).Label.Units = 'normalized';
  hd.gain.ob(7,1).Label.Position = [-.8 1 1] .* hd.gain.ob(7,1).Label.Position;
  
  hd.gain.ax(8) = axes;
  hd.gain.ax(8).Position = hd.gain.ax(4).Position - [.13 -.1 .2 .2];
  hd.gain.ax(8).CLim = hd.gain.ax(4).CLim;
  hd.gain.ob(8,1) = colorbar;
  hd.gain.ob(8,1).Label.String = 'fit - data';
  hd.gain.ob(8,1).Label.Units = 'normalized';
  hd.gain.ob(8,1).Label.Position = [-.8 1 1] .* hd.gain.ob(8,1).Label.Position;

  [hd.gain.ob(7,1).Label.FontSize, hd.gain.ob(7,1).FontSize, ...
   hd.gain.ob(8,1).Label.FontSize, hd.gain.ob(8,1).FontSize, ...
   hd.gain.ax.FontSize] = deal(14);
 
  [hd.gain.ax(1).Title.FontSize, ...
   hd.gain.ax(2).Title.FontSize, ...
   hd.gain.ax(3).Title.FontSize] = deal(14);
   
  [hd.gain.ax(1).Title.FontWeight, ...
   hd.gain.ax(2).Title.FontWeight, ...
   hd.gain.ax(3).Title.FontWeight] = deal('normal');

  
  [hd.gain.ax.Visible] = deal('off');
%   [hd.gain.ax.Title.Visible] = deal('on');
  set(([hd.gain.ax.Title]),'Visible','on')
  set(([hd.gain.ax.Title]),'Position',[1 1 .7] .* hd.gain.ax(1).Title.Position)
  
end



function [retained,original] = discardTrial(original,folderstring)
% DISCARDTRIAL removes data obtained from trials consider "poorly fit" etc.
%
% This function takes a "result" structure as its only input argument, and
% requires an XLSX file containing the indices of trials to be discarded.
% It then removes those bad trials, and returns a reduced structure as its
% primary output argument, "retained". For easy access, it also returns its
% original input.
%
% This function can safely be executed multiple times, as "discarded"
% trials are in fact just set to empty, so their array indices are not 
% reassigned to other trials.
%
% Written by Florian Alexander Dehmelt, Tuebingen University, in 2018.


  % Identify and load an Excel file containing the indices of all unwanted
  % trials to be discarded. This file must contain four columns of numbers
  % (from left to right: fish no., eye no., stimulus type, and repetition.
  % It may contain as many text comments as desired; these will be ignored.
  folder.badtrial = folderstring;
%   folder.badtrial = 'C:\Users\fdehmelt\Dropbox\ThinkPad Dropbox\MATLAB\201808 SphereAnalysis\Sample Data\D series\';
  file.badtrial = 'DiscardTrial.xlsx';
  badtrial = xlsread([folder.badtrial,file.badtrial]);
  
  % The "retained" structure is the same as the "original" one...
  retained = original;
  
  % ...except that unwanted trials are overwritten by empty values.
  % They thus appear the same way as non-existent trials already did.
  for trial = 1:size(badtrial,1)
    
    % Find the position of the bad trial within the structure.
    fish       = badtrial(trial,1);
    eye        = badtrial(trial,2);
    stimtype   = badtrial(trial,3) + 1; % The +1 converts index >=0 to >=1.
    repetition = badtrial(trial,4);
    
    % Overwrite all field entries with "[]".
    fieldname = fieldnames(retained);
    for field = 1:numel(fieldname)
      retained(fish,eye,stimtype,repetition).(fieldname{field}) = [];
    end
  
  end
  
end



function [altered,original] = zeroTrial(original,folderstring)
% ZEROTRIAL sets gains to zero based on visual inspection of trials.
%
% This function takes a "result" structure as its only input argument, and
% requires an XLSX file containing the indices of trials to be set to zero.
% It then adjusts those gains, and returns an altered structure as its
% primary output argument, "retained". For easy access, it also returns its
% original input.
%
% This function can safely be executed multiple times, as "altered"
% trials are in fact just set to zero, so their array indices are not 
% reassigned to other trials.
%
% Written by Florian Alexander Dehmelt, Tuebingen University, in 2018.


  % Identify and load an Excel file containing the indices of all unwanted
  % trials to be discarded. This file must contain four columns of numbers
  % (from left to right: fish no., eye no., stimulus type, and repetition.
  % It may contain as many text comments as desired; these will be ignored.
  folder.badtrial = folderstring;
%   folder.badtrial = 'C:\Users\fdehmelt\Dropbox\ThinkPad Dropbox\MATLAB\201808 SphereAnalysis\Sample Data\D series\';
  file.badtrial = 'ZeroTrial.xlsx';
  badtrial = xlsread([folder.badtrial,file.badtrial]);
  
  % The "altered" structure is the same as the "original" one...
  altered = original;
  
  % ...except that gains of visually identified trials are set to zero.
  % They thus appear the same way as non-existent trials already did.
  for trial = 1:size(badtrial,1)
    
    % Find the position of the bad trial within the structure.
    fish       = badtrial(trial,1);
    eye        = badtrial(trial,2);
    stimtype   = badtrial(trial,3) + 1; % The +1 converts index >=0 to >=1.
    repetition = badtrial(trial,4);
    
    % Overwrite the gain field entries with "0", retain all other values.
    altered(fish,eye,stimtype,repetition).gain = 0;
  
  end
  
end



function combinedvalue = vonmisesfisher9param(param,predictor)
%VONMISESFISHER9PARAM Evaluate two overlapping von Mises-Fisher distributions.
%
%   VALUE = VONMISESFISHER9PARAM computes the value of two(!) overlapping, 
%   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 = '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 = 2;
  value = NaN(size(predictor,1),numberofdistributions);
  
  for k = 1:numberofdistributions
    
    offset        = param(1);
    amplitude     = param(4*k-2);
    meanazimuth   = param(4*k-1);
    meanelevation = param(4*k);
    concentration = param(4*k+1);
    
    [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;
    
    % 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



function [x,y,z] = geo2car(azimuth,elevation,radius,convention)
 
  switch convention
    case 'CCW'
      x = radius .* cosd(elevation).*cosd(azimuth);
      y = radius .* cosd(elevation).*sind(azimuth);
      z = radius .* sind(elevation);
    case 'CW'
      x = radius .* cosd(elevation).*cosd(azimuth);
      y = radius .* cosd(elevation).*-sind(azimuth);
      z = radius .* sind(elevation);
    otherwise
      error('Convention not supported.')
  end
  
end



function [azimuth,elevation,radius] = car2geo(x,y,z,convention)
  
  azimuth = atand(y./x) + 180*((x<0).*(y>=0) - (x<0).*(y<0));
        
  switch convention
    case 'CCW'
    case 'CW'
      azimuth = -azimuth;
    otherwise
      error('Convention not supported.')
  end
  
  radius    = sqrt(x.^2+y.^2+z.^2);
  elevation = asind(z./radius);
    
end



function [lefteyehandle,righteyehandle,bodyhandle] = ...
  plot3dFish(scale,yaw,pitch,roll)

  % This corresponds to version 20170505a, a.k.a. plot3dFish20170505a.m.

  set(gcf,'Color',[1 1 1])

  resolution = 10;
  eye.colour = .1*[1 1 1];
  mainbody.colour = .7*[1 1 1];

  eye.x = scale * (cos(-pi:pi/resolution:3*pi-pi/resolution)');
  eye.y = scale * (0.7*sin(0:pi/resolution:4*pi-pi/resolution)');
  eye.z = scale * ([.7*ones(2*resolution,1); ...
                    1*ones(2*resolution,1)]);

  N = 2*resolution;

  lefteye.vertices = [eye.x, eye.y, eye.z];
  righteye.vertices = [eye.x, -eye.y, eye.z];

  lefteye.faces = [[1:2*resolution 1]; ...
                  [1+2*resolution:4*resolution 1+2*resolution]];


  for k = 1:N

    lefteye.faces = [lefteye.faces; ...
                    [mod([k-1 k],N)+1, ...
                     mod([k k-1],N)+1+N, ...
                     mod(k-1,N)+1, ...
                     NaN(1,N-4)]];

  end

  lefteye.facevertexcdata = repmat(eye.colour,[size(lefteye.faces,1) 1]);

  righteye.faces = lefteye.faces;
  righteye.facevertexcdata = lefteye.facevertexcdata;

  lefteyehandle  = patch(lefteye);
  righteyehandle = patch(righteye);
  set([lefteyehandle,righteyehandle],'FaceColor','flat','EdgeColor','none')

  mainbody.x = scale * (-.1 + [1.0*cos(-pi:pi/resolution:pi-pi/resolution), ...
                               0.4*cos(-pi:pi/resolution:pi-pi/resolution)]');
  mainbody.y = scale * ([0.7*sin(0:pi/resolution:2*pi-pi/resolution), ...
                         0.2*sin(0:pi/resolution:2*pi-pi/resolution)]');
  mainbody.z = scale * ([1.1*ones(2*resolution,1); ...
                        -7*ones(2*resolution,1)]);

  N = 2*resolution;

  body.vertices = [mainbody.x, mainbody.y, mainbody.z];

  body.faces = [[1:2*resolution 1]; ...
                [1+2*resolution:4*resolution 1+2*resolution]];

  for k = 1:N

    body.faces = [body.faces; ...
                  [mod([k-1 k],N)+1, ...
                   mod([k k-1],N)+1+N, ...
                   mod(k-1,N)+1, ...
                   NaN(1,N-4)]];

  end

  body.facevertexcdata = repmat(mainbody.colour,[size(body.faces,1) 1]);

  bodyhandle = patch(body);
  set(bodyhandle,'FaceColor','flat','EdgeColor','none')

  rotate(lefteyehandle,[1 0 0],80)
  rotate(lefteyehandle,[0 0 1],-10)
  rotate(righteyehandle,[1 0 0],-80)
  rotate(righteyehandle,[0 0 1],10)
  rotate(bodyhandle,[0 1 0],-90)
  
  rotate(lefteyehandle,[0 0 1],yaw)
  rotate(righteyehandle,[0 0 1],yaw)
  rotate(bodyhandle,[0 0 1],yaw)
  
  rotate(lefteyehandle,[0 1 0],pitch)
  rotate(righteyehandle,[0 1 0],pitch)
  rotate(bodyhandle,[0 1 0],pitch)
  
  rotate(lefteyehandle,[1 0 0],roll)
  rotate(righteyehandle,[1 0 0],roll)
  rotate(bodyhandle,[1 0 0],roll)

end



function hd = displayStimuli(typestring)
  
warning('off','MATLAB:interp1:UsePCHIP')
% The line above deactivates 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)

%   close all
%   clf
  
  switch typestring
    
    case 'disk'
      file.folder = '..\..\data\Fig3 default disk data\Data D Series (neue Experimente)\';
      datarange = 9:46;
  
    case 'area'
      file.folder = '..\..\data\Fig5 size data\Data A Series (A-Series)\';
      datarange = 2:43;
      datarange = 2:8;

    case 'freq'
      file.folder = '..\..\data\Fig5 frequency data\Data V Series (V-Serie)\';
      datarange = 2:43;
      datarange = 2;

    otherwise
      error('Argument of displayStimuli.m must be ''disk'', ''area'', or ''freq''.')
      
  end
  
  file.stimulus = 'Stimulus.xlsx';
  [~,stim.tableconvention,~] = xlsread([file.folder,file.stimulus],'C1:C1');
  stim.azimuth = xlsread([file.folder,file.stimulus],['C',num2str(datarange(1)),':C',num2str(datarange(end))]);
  stim.elevation = xlsread([file.folder,file.stimulus],['D',num2str(datarange(1)),':D',num2str(datarange(end))]);
  stim.radius = xlsread([file.folder,file.stimulus],['E',num2str(datarange(1)),':E',num2str(datarange(end))]);
  
  % compute disk boundaries
  diskcentre.geo = [stim.azimuth, stim.elevation];
  diskcentre.angle = stim.radius;
%   diskcentre.geo = [  45   20; ...
%                       30    5; ...
%                       60  -45];
%   diskcentre.angle = 40 * [ones(size(diskcentre.geo,1),1)];
  diskcentre.azim = diskcentre.geo(:,1);
  diskcentre.elev = diskcentre.geo(:,2);
  diskcentre.rads = ones(size(diskcentre.geo(:,1)));
  resolution = .001;
  numdisk = size(diskcentre.geo,1);
  for disk = 1:numdisk
    [boundary.azim{disk}, ...
     boundary.elev{disk}] = diskBoundary(diskcentre.azim(disk), ...
                            diskcentre.elev(disk), diskcentre.angle(disk), resolution);
    [boundary.x{disk}, ...
     boundary.y{disk}, ...
     boundary.z{disk}] = geo2car(boundary.azim{disk}, boundary.elev{disk}, ...
                                 ones(size(boundary.azim{disk})),'CW');
  end
  
  
  % create figure
  hd.default.fg = figure(1);
  hd.default.fg.Name = 'This is a sphere.';
  hd.default.fg.Position = [50  50 500 800];
  hd.default.fg.Color = [1 1 1];
  
  hd.default.ax(1) = axes;
  hd.default.ax(1).Position = [0 0 1 1];
  
  hold on

    % plot background sphere with shading
    hd.default.globe = displaySphere;


    % plot disk boundaries
  %   diskcolour = cbrewer('seq','Greens',numdisk*2);
    diskcolour = flipud(cbrewer('seq','Greys',numdisk*2));
  %   diskcolour = diskcolour(randperm(numdisk),:);
    croptext.xpos = [1.27 1.32 1.35 1.34 1.25 0.7 -1.58];
    croptext.ypos = [0 0 0 0 0 0 0];
    croptext.zpos = [.365 .18 .0 -.18 -.45 -.95 -.85];
    for disk = 1:numdisk
  %     plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',3,'Color',0*[1 1 1])
  %     plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',2,'Color',1*[1 1 1])
      plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',1,'Color',diskcolour(disk,:))
  %     plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',1,'Color',.2*[1 1 1])
      croptext.label(disk) = text(croptext.xpos(disk),croptext.ypos(disk),croptext.zpos(disk), ...
                                  [num2str(diskcentre.angle(disk),'%.1f'),'°'],'FontSize',12);
    end


  hold off
  
end



function globe = displaySphere
  
  [surface.x,surface.y,surface.z] = sphere(100);
  
  reference.geo = [30;20;1];  % desired azimuth, elevation, radius of maximum colour value
  reference.azim   = reference.geo(1);
  reference.elev   = reference.geo(2);
  reference.rads = ones(size(reference.geo(1)));
  [reference.x,reference.y,reference.z] = geo2car(reference.azim,reference.elev,reference.rads,'CW');
  scalarproduct = [surface.x(:),surface.y(:),surface.z(:)] * [reference.x;reference.y;reference.z];
  surface.c = reshape(scalarproduct,size(surface.x));
  surface.c = 10.^(surface.c);

  globe = surf(surface.x,surface.y,surface.z,surface.c);
  shading flat
  axis square

  greyscale = (.8:.01:.99)' * [1 1 1];  
  colormap(greyscale)
  
end



function [azimuthOut, elevationOut] = diskBoundary(azimuthIn, elevationIn, angleIn, 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;
  % This equation is used to compute border.azimuth below.

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

  intended  = -90 : border.resolution : 90;
%   limit = cosd([centre.elevation-border.angle centre.elevation+border.angle]);
%   intended = acosd(min(limit) : border.resolution : max(limit));
  % "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]);
  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;
  
  azimuthOut   = border.azimuth';
  elevationOut = border.elevation';

end



function hd = displayPattern(typestring,spatialfreq)

  %   typestring = 'freq';

  %   switch typestring
  %     
  %     case 'disk'
  %       file.folder = '..\..\data\Fig3 default disk data\Data D Series (neue Experimente)\';
  %       datarange = 9:46;
  %   
  %     case 'area'
  %       file.folder = '..\..\data\Fig5 size data\Data A Series (A-Series)\';
  %       datarange = 2:43;
  %       datarange = 2:8;
  % 
  %     case 'freq'
  %       file.folder = '..\..\data\Fig5 frequency data\Data V Series (V-Serie)\';
  %       datarange = 2:43;
  %       datarange = 2;
  % 
  %     otherwise
  %       error('Argument of displayStimuli.m must be ''disk'', ''area'', or ''freq''.')
  %       
  %   end

%   switch typestring
%     case 'area'
%       spatialfreq = 5.5 * ones(1,7);
%     case 'freq'
%       % spatialfreq = [1.5 3 5 7 9 11 20];
%       spatialfreq = logspace(0.18,1.2,7);
%     otherwise
%       error('Argument of displayStimuli.m must be ''area'' or ''freq''.')
%   end
  
  grid.colour = .3*[1 1 1];
  grid.width = 1;
  grid.height = .7;
  grid.xspace = .29;
  grid.yspace = .27 + .25;


  grid.x0 = [-grid.width-grid.xspace/2,   +grid.xspace/2, ...
             -grid.width/2*3-grid.xspace, -grid.width/2,   +grid.width/2+grid.xspace,...
             -grid.width-grid.xspace/2,   +grid.xspace/2];

  grid.y0 = [ [1 1]  * (+grid.height/2+grid.yspace), ...
             [1 1 1] * (-grid.height/2), ...
              [1 1]  * (-grid.height/2*3-grid.yspace)];


  showazimuthrange = [-180 -135];
  showazimuthfraction = abs(diff(showazimuthrange))/360;
  
  hd.pattern.fg = figure();
  hd.pattern.ax = axes();
  hold on
  dy = grid.height;
  for gridindex = 1:7
    numbar(gridindex) = ceil(grid.width*360*spatialfreq(gridindex)*showazimuthfraction);
    dx = grid.width / (360*spatialfreq(gridindex)*showazimuthfraction) / 2;
    for barindex = 1:numbar(gridindex)
      x0 = grid.x0(gridindex) + (barindex-1) * dx * 2;
      xmax = grid.x0(gridindex) + grid.width;
      y0 = grid.y0(gridindex);
      darkbar(gridindex,barindex) = patch([x0*[1 1],min(x0+dx,xmax)*[1 1]], ...
                                          [y0,y0+dy,y0+dy,y0], grid.colour);
    end
    outline(gridindex) = plot(grid.x0(gridindex)*[1 1 1 1 1] + grid.width*[0 0 1 1 0], ...
                              grid.y0(gridindex)*[1 1 1 1 1] + grid.height*[0 1 1 0 0], ...
                              'Color', grid.colour);
    textlabel(gridindex) = text(grid.x0(gridindex)+grid.width/2, ...
                                grid.y0(gridindex)+grid.height/30, ...
                                [num2str(spatialfreq(gridindex),'%.2f'),' cpd']);
    textlabel(gridindex).HorizontalAlignment = 'center';
    textlabel(gridindex).VerticalAlignment = 'top';
    textlabel(gridindex).FontSize = 12;
  end
  hold off

  hd.pattern.ax.PlotBoxAspectRatio = [grid.width, grid.height, grid.height];
  [hd.pattern.ax.XAxis.Visible, hd.pattern.ax.YAxis.Visible] = deal('off');
  hd.pattern.ax.XLim = [-1 1] * 2.2;
%   hd.pattern.ax.XLim = [0 .25] * 2.2;
  hd.pattern.ax.YLim = [-1 1] * 1.8;
  
end