clear all
close all

grating = displayPattern(360/22);
velocity = displayCosine(12.5);
[leftdisk,leftparam] = displayDiskStimuli('left');
[rightdisk,rightparam] = displayDiskStimuli('right');
singledisk = displayDiskStimuli('single');
[diskhandle,diskdata,fitdata] = gainAnalysisDisk('Resolution',50);
variant = 3; % 1 = regular, 2 = rotated, 3 = combined/corrected     %%%%% !!!! %%%%%

hd.collage.fg = figure();
hd.collage.fg.Name = 'Sphere manuscript (Suppl Fig); Yoking dependence on stimulus location';
hd.collage.fg.Color = [1 1 1];
hd.collage.fg.Position = [0 40 1070 960];

% Populate the combined figure from previously generated partial figures...
hd.collage.ax(1,1) = copyobj(grating.pattern.ax(1), hd.collage.fg);
hd.collage.ax(1,2) = copyobj(velocity.cosine.ax(1), hd.collage.fg);
hd.collage.ax(1,3) = copyobj(leftdisk.default.ax(1), hd.collage.fg);
hd.collage.ax(1,4) = copyobj(rightdisk.default.ax(1), hd.collage.fg);
hd.collage.ax(1,5) = copyobj(singledisk.default.ax(1), hd.collage.fg);
% hd.collage.ax(2,1) = copyobj(diskdata.gain.ax(7), hd.collage.fg);
hd.collage.ax(2,1) = axes;
for subpanel = 2:4
  hd.collage.ax(2,subpanel) = copyobj(diskhandle.gain(variant).ax(subpanel-1), hd.collage.fg);
end
for subpanel = 5:7
  hd.collage.ax(2,subpanel) = copyobj(hd.collage.ax(2,subpanel-3), hd.collage.fg);
end

% ...get rid of the partial figures once pillaged...
% close(1:7)

% ...then adjust the size, position and other properties of those objects.


[hd.collage.ax(1,1:5).Units, hd.collage.ax(2,1:7).Units] = deal('pixels');
hd.collage.ax(1,1).Position = [100 850 180  70];
hd.collage.ax(1,2).Position = [130 730 140  55];
hd.collage.ax(1,3).Position = [350 730 200 200];
hd.collage.ax(1,4).Position = hd.collage.ax(1,3).Position + [250 0 0 0];
hd.collage.ax(1,5).Position = hd.collage.ax(1,3).Position + [520 0 0 0];
hd.collage.ax(2,1).Position = [ 50 240  50 160];
hd.collage.ax(2,2).Position = [ 90 280 350 350];
for subpanel = 3:4
  hd.collage.ax(2,subpanel).Position = hd.collage.ax(2,subpanel-1).Position + [320 0 0 0];
end
for subpanel = 5:7
  hd.collage.ax(2,subpanel).Position = hd.collage.ax(2,subpanel-3).Position - [0 280 0 0];
end

darkgreyscale = (.8:.01:.99)' * [1 1 1];  
lightgreyscale = (.9:.01:.99)' * [1 1 1];  
[hd.collage.ax(1,3:4).Colormap] = deal(lightgreyscale);
[hd.collage.ax(1,5).Colormap] = deal(darkgreyscale);


  
% % Shared properties of all axis objects anywhere:
[hd.collage.ax(1,1:5).FontSize, ...
 hd.collage.ax(2,1:7).FontSize] = deal(15);
[hd.collage.ax(1,1:5).LabelFontSizeMultiplier, ...
 hd.collage.ax(1,1:5).TitleFontSizeMultiplier, ...
 hd.collage.ax(2,1:7).LabelFontSizeMultiplier, ...
 hd.collage.ax(2,1:7).TitleFontSizeMultiplier] = deal(1);


% Create content of first subpanel of panel b (colorbar):
hd.collage.ax(2,1).CLim = hd.collage.ax(2,2).CLim;
hd.collage.ob(2,1,1) = colorbar;
hd.collage.ob(2,1,1).Label.String = 'OKR gain';
hd.collage.ob(2,1,1).Label.Units = 'normalized';
hd.collage.ob(2,1,1).AxisLocation = 'in';
hd.collage.ax(2,1).Visible = 'off';
hd.collage.ob(2,1,1).Units = 'pixels';
hd.collage.ob(2,1,1).Position(3) = 9;
hd.collage.ob(2,1,1).FontSize = hd.collage.ax(2,1).FontSize;
hd.collage.ob(2,1,1).Label.Position = [-.8 1 1] .* hd.collage.ob(2,1,1).Label.Position;


% Add and format titles to panel a:
for panel = 1
  hd.collage.ax(panel,1).Title.String = 'stimulus pattern';
  hd.collage.ax(panel,3).Title.String = '----------------------- stimulus crop -----------------------';
  hd.collage.ax(panel,3).Title.Interpreter = 'Latex';
%   hd.collage.ax(panel,4).Title.String = {'stimulus crop','(right hemisphere)'};
  hd.collage.ax(panel,2).XLabel.String = 'time (sec)';
  hd.collage.ax(panel,2).YLabel.String = {'velocity','(deg/sec)'};
  hd.collage.ax(panel,3).XLabel.String = 'left hemisphere';
  hd.collage.ax(panel,4).XLabel.String = 'right hemisphere';
  [hd.collage.ax(panel,3).XLabel.Visible, ...
   hd.collage.ax(panel,4).XLabel.Visible] = deal('on');
  hd.collage.ax(panel,5).Title.String = 'stimulus';
  for subpanel = 1:5
    hd.collage.ax(panel,subpanel).Title.FontWeight = 'normal';
    [hd.collage.ax(panel,subpanel).Title.Units, ...
     hd.collage.ax(panel,subpanel).XLabel.Units] = deal('pixels');
  end
%   hd.collage.ax(panel,1).Title.Position(2) = hd.collage.ax(panel,1).Position(4) - 14;
  hd.collage.ax(panel,1).XLabel.Position(2) = hd.collage.ax(panel,1).XLabel.Position(2) + 25;
%   hd.collage.ax(panel,2).XLabel.Position(2) = hd.collage.ax(panel,2).XLabel.Position(2) - 65;
  hd.collage.ax(panel,2).XLabel.Position(2) = hd.collage.ax(panel,2).XLabel.Position(2) + 25;
  hd.collage.ax(panel,3).Title.Position(1) = hd.collage.ax(1,3).Title.Position(1) + 126;
  for subpanel = 3:4
    hd.collage.ax(panel,subpanel).Title.Position(2) = hd.collage.ax(panel,subpanel).Position(4) - 2;
    hd.collage.ax(panel,subpanel).XLabel.Position(2) = hd.collage.ax(panel,subpanel).XLabel.Position(2) + 23;
  end
%   hd.collage.ax(panel,5).Title.Position(2) = hd.collage.ax(panel,4).Position(4) + 2;
end
% axes(hd.collage.ax(1,3))
% leftline = text(hd.collage.ax(1,3).Title.Position(1) - 50, ...
%                 hd.collage.ax(1,3).Title.Position(2) + 20, ...
%                 '-----------------------','Interpreter','Latex','HorizontalAlignment','right','Units','pixels')


% Force visibility of axis ticks on bottom left (sine) plot of panel a:
% hd.collage.ax(1,2).XLim = [-13 13];
hd.collage.ax(1,2).YLim = [-13 13];
  

% Format titles on top row of panels b to d:
for subpanel = 2:4
  hd.collage.ax(2,subpanel).Title.Units = 'pixels';
  hd.collage.ax(2,subpanel).Title.Position(2) = hd.collage.ax(2,subpanel).Title.Position(2) - 25;
end
% Correct for odd treatment of indices (or their absence):
hd.collage.ax(2,3).Title.Position(2) = hd.collage.ax(2,3).Title.Position(2) + 9;


% Remove titles from bottom row of panels b to d:
for subpanel = 5:7
  hd.collage.ax(2,subpanel).Title.Visible = 'off';
end


% % Make fish on panels b to d more visible by adjusting sphere alpha:
% [hd.collage.ax(

% Adjust camera perspective on the lower row of panels b to d:
[hd.collage.ax(2,5:7).CameraPosition] = deal([19.5 0 0]);


% Display panel lettering:
hd.collage.ax(4,1) = axes;
hd.collage.ax(4,1).Color = 'none';
hd.collage.ax(4,1).Position = [0 0 1 1];
hd.collage.ax(4,1).Units = 'pixels';
hd.collage.ax(4,1).Visible = 'off';
hold on
hd.collage.lt(1) = text(hd.collage.ax(1,1).Position(1) - 85 , ...
                        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,2).Position(1) - 50, ...
                        hd.collage.ax(2,2).Position(2) + hd.collage.ax(2,2).Position(4) - 10, ...
                        'b', 'Units', 'pixels');
hd.collage.lt(3) = text(hd.collage.ax(2,3).Position(1) + 35, ...
                        hd.collage.ax(2,3).Position(2) + hd.collage.ax(2,3).Position(4) - 10, ...
                        'c', 'Units', 'pixels');
hd.collage.lt(4) = text(hd.collage.ax(2,4).Position(1) + 35, ...
                        hd.collage.ax(2,4).Position(2) + hd.collage.ax(2,4).Position(4) - 10, ...
                        'd', 'Units', 'pixels');
hold off
[hd.collage.lt(:).FontSize] = deal(15);
[hd.collage.lt(:).FontWeight] = deal('bold');

if variant == 3
  plotAsymmetry(diskdata(variant))
  plotYoking(diskdata(variant))
end
mercatorPanel20190520

delete(hd.collage.ax(1,1:2)) % remove panels about basic stimulus pattern
close(1:8)
for subpanel = 3:4
  hd.collage.ax(1,subpanel).Position = hd.collage.ax(1,subpanel).Position - [200 0 0 0];
end
hd.collage.ax(1,5).Position = hd.collage.ax(1,5).Position - [100 0 0 0];




% [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);




function plotYoking(diskdata)
  
  col.barface = [.2 .8 .4];
  col.neutralbarface = .8*[1 1 1];
  col.baredge = 'none';
  col.errorbar = .3*[1 1 1];

  hd.yoke.fg = figure();
  hd.yoke.fg.Color = [1 1 1];
  hd.yoke.fg.Position = [0 100 1900 390];

  hd.yoke.ax(1,1) = axes();
  % hd.yoke.ax(1,2) = axes();

  [hd.yoke.ax(:,:).Units] = deal('pixels');

  hd.yoke.ax(1,1).Position = [100 70 1780 300];

  y5a = diskdata.yoking.regular;
  y5b = diskdata.yoking.rotated;
  y5c = diskdata.yoking.corrected;
  y5d = diskdata.yoking.black;
  % e5a = mscdisk.botheyes.stdgain;
  % e5b = upside.botheyes.stdgain;
  e5a = NaN(size(y5a));
  e5b = NaN(size(y5b));
  e5c = NaN(size(y5c));
  e5d = NaN(size(y5d));

%   pvalue.asymmetry = signrank(y5a,y5b);

  x5a = (1:numel(y5a))' - .24;
  x5b = (1:numel(y5b))' - .08;
  x5c = (1:numel(y5b))' + .08;
  x5d = (1:numel(y5c))' + .24;


  axes(hd.yoke.ax(1,1))
  hold on
  hd.yoke.barregular = bar(x5a,y5a);
  hd.yoke.barcontrol = bar(x5b,y5b);
  hd.yoke.barcorrect = bar(x5c,y5c);
  hd.yoke.barblack = bar(x5d,y5d);
  hd.yoke.errorregular = errorbar(x5a,y5a,e5a,'LineStyle','none');
  hd.yoke.errorcontrol = errorbar(x5b,y5b,e5b,'LineStyle','none');
  hd.yoke.errorcorrect = errorbar(x5c,y5c,e5c,'LineStyle','none');
  hd.yoke.errorblack = errorbar(x5d,y5d,e5d,'LineStyle','none');
  hd.yoke.explanation = text(29,.351,{['regular (light grey), rotated (dark grey), ', ...
                                      'corrected (green), reflections screened (black)']}, ...
                             'HorizontalAlignment','left');

  hold off

%   [hd.yoke.ax(1,:).XTick] = deal([]);
  hd.yoke.ax(1,1).XTick = mean([x5b,x5c],2);
  
  oldlabel = hd.yoke.ax(1,1).XTickLabel;
  for entry = 1:numel(oldlabel)
    if entry < 8
      newlabel{entry} = ['H',num2str(str2num(oldlabel{entry}))];
    else
      newlabel{entry} = ['D',num2str(str2num(oldlabel{entry})-7)];
    end
  end
  hd.yoke.ax(1,1).XTickLabel = newlabel;
%   hd.yoke.ax(1,1).XTickLabel{1} = '';
  

  hd.yoke.barregular.FaceColor = col.neutralbarface;
  hd.yoke.barcontrol.FaceColor = col.neutralbarface*3/4;
  hd.yoke.barcorrect.FaceColor = col.barface;
  hd.yoke.barblack.FaceColor = col.neutralbarface*1/4;

  [hd.yoke.barregular.BarWidth, ...
   hd.yoke.barcontrol.BarWidth, ...
   hd.yoke.barcorrect.BarWidth, ...
   hd.yoke.barblack.BarWidth] = deal(.16);

  [hd.yoke.errorregular.Color, ...
   hd.yoke.errorcontrol.Color, ...
   hd.yoke.errorcorrect.Color, ...
   hd.yoke.errorblack.Color] = deal(col.errorbar);

  [hd.yoke.errorregular.LineWidth, ...
   hd.yoke.errorcontrol.LineWidth, ...
   hd.yoke.errorcorrect.LineWidth, ...
   hd.yoke.errorblack.LineWidth] = deal(1);

  hd.yoke.ax(1,1).XLim = [.5 46];
  % hd.yoke.ax(1,1).YLim = hd.yoke.ax(1,2).YLim;
  hd.yoke.ax(1,1).YLabel.String = 'yoking index  (g_L-g_R)/(g_L+g_R)';
  for panel = 1:1
    hd.yoke.ax(1,panel).YLabel.Units = 'pixels';
    hd.yoke.ax(1,panel).YLabel.Position(1) = hd.yoke.ax(1,panel).YLabel.Position(1) - 20;
  end
  hd.yoke.ax(1,1).XLabel.String = 'stimulus types';
  % hd.yoke.ax(1,2).XLabel.String = 'stimulus types';
  [hd.yoke.explanation.FontSize, ...
   hd.yoke.pvalue.FontSize] = deal(13);
  for panel = 1:1
    hd.yoke.ax(1,panel).XLabel.FontSize = 13;
    hd.yoke.ax(1,panel).FontSize = 13;
    hd.yoke.ax(1,panel).LabelFontSizeMultiplier = 1;
    hd.yoke.ax(1,panel).YGrid = 'on';
    hd.yoke.ax(1,panel).TickLength = [0 0];
  end
  
  hd.yoke.ax(1,1).YMinorGrid = 'on';
  
end



function plotAsymmetry(diskdata)
  
  col.barface = [.2 .8 .4];
  col.neutralbarface = .8*[1 1 1];
  col.baredge = 'none';
  col.errorbar = .3*[1 1 1];

  hd.asymindex.fg = figure();
  hd.asymindex.fg.Color = [1 1 1];
  hd.asymindex.fg.Position = [0 100 1900 390];

  hd.asymindex.ax(1,1) = axes();
  % hd.asymindex.ax(1,2) = axes();

  [hd.asymindex.ax(:,:).Units] = deal('pixels');

  hd.asymindex.ax(1,1).Position = [100 70 1780 300];

  y5a = diskdata.asymmetry.A1;
  y5b = diskdata.asymmetry.A2;
  % e5a = mscdisk.botheyes.stdgain;
  % e5b = upside.botheyes.stdgain;
  e5a = NaN(size(y5a));
  e5b = NaN(size(y5b));

  pvalue.asymmetry = signrank(y5a,y5b);

  x5a = (1:numel(y5a))' - .2;
  x5b = (1:numel(y5b))' + .2;


  axes(hd.asymindex.ax(1,1))
  hold on
  hd.asymindex.barregular = bar(x5a,y5a);
  hd.asymindex.barcontrol = bar(x5b,y5b);
  hd.asymindex.errorregular = errorbar(x5a,y5a,e5a,'LineStyle','none');
  hd.asymindex.errorcontrol = errorbar(x5b,y5b,e5b,'LineStyle','none');
  hd.asymindex.explanation = text(24,.351,{'biological asymmetry shown in green'},'HorizontalAlignment','center');
%   hd.asymindex.explanation = text(24,.351,{'biological asymmetry shown in green', ...
%     '(error bars show s.e.m.)'},'HorizontalAlignment','center');
%   hd.asymindex.pvalue = text(36,.351,{['p = ',num2str(pvalue.asymmetry)], ...
%     '(Wilcoxon signed-rank test)'},'HorizontalAlignment','center');

  hold off

%   [hd.asymindex.ax(1,:).XTick] = deal([]);
  hd.asymindex.ax(1,1).XTick = mean([x5a,x5b],2);
  
  oldlabel = hd.asymindex.ax(1,1).XTickLabel;
  for entry = 1:numel(oldlabel)
    if entry < 8
      newlabel{entry} = ['H',num2str(str2num(oldlabel{entry}))];
    else
      newlabel{entry} = ['D',num2str(str2num(oldlabel{entry})-7)];
    end
  end
  hd.asymindex.ax(1,1).XTickLabel = newlabel;

  hd.asymindex.barregular.FaceColor = col.neutralbarface;

  hd.asymindex.barcontrol.FaceColor = col.barface;

  [hd.asymindex.barregular.BarWidth, ...
   hd.asymindex.barcontrol.BarWidth] = deal(.35);

  [hd.asymindex.errorregular.Color, ...
   hd.asymindex.errorcontrol.Color] = deal(col.errorbar);

  [hd.asymindex.errorregular.LineWidth, ...
   hd.asymindex.errorcontrol.LineWidth] = deal(1);


  hd.asymindex.ax(1,1).XLim = [.5 46];
  % hd.asymindex.ax(1,1).YLim = hd.asymindex.ax(1,2).YLim;
  hd.asymindex.ax(1,1).YLabel.String = 'asymmetries A1 (grey), A2+A3 (green)';
%   hd.asymindex.ax(1,1).YLabel.String = 'mean OKR gain  (pooled across eyes)';
  % hd.asymindex.ax(1,2).YLabel.String = 'mean OKR gain  (right eye)';
  for panel = 1:1
    hd.asymindex.ax(1,panel).YLabel.Units = 'pixels';
    hd.asymindex.ax(1,panel).YLabel.Position(1) = hd.asymindex.ax(1,panel).YLabel.Position(1) - 20;
  end
  hd.asymindex.ax(1,1).XLabel.String = 'stimulus types';
  % hd.asymindex.ax(1,2).XLabel.String = 'stimulus types';
  [hd.asymindex.explanation.FontSize, ...
   hd.asymindex.pvalue.FontSize] = deal(13);
  for panel = 1:1
    hd.asymindex.ax(1,panel).XLabel.FontSize = 13;
    hd.asymindex.ax(1,panel).FontSize = 13;
    hd.asymindex.ax(1,panel).LabelFontSizeMultiplier = 1;
    hd.asymindex.ax(1,panel).YGrid = 'on';
    hd.asymindex.ax(1,panel).TickLength = [0 0];
  end

  
end




function hd = displayPattern(barwidth)

  spatialfreq = 180/barwidth;
  
  grid.colour = .3*[1 1 1];
  grid.width = 360;
  grid.height = 180;
  grid.x0 = -180;
  grid.y0 = -90;

  hd.pattern.fg = figure();
  hd.pattern.ax = axes();
  hold on
  dy = grid.height;
  numbar = ceil(grid.width*spatialfreq);
  dx = grid.width / spatialfreq / 2;
  for barindex = 1:numbar
    x0 = grid.x0 + (barindex-1) * dx * 2;
    xmax = grid.x0 + grid.width;
    y0 = grid.y0;
    darkbar(barindex) = patch([x0*[1 1],min(x0+dx,xmax)*[1 1]], ...
                              [y0,y0+dy,y0+dy,y0], grid.colour);
  end
  [darkbar(:).EdgeColor] = deal('none');
  outline = plot(grid.x0*[1 1 1 1 1] + grid.width*[0 0 1 1 0], ...
                 grid.y0*[1 1 1 1 1] + grid.height*[0 1 1 0 0], ...
                 'Color', grid.colour);                            
  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 = [-180 180];
  hd.pattern.ax.YLim = [-90 90];
  hd.pattern.ax.XTick = [-180 180];
  hd.pattern.ax.YTick = [-90 0 90];
  
  hd.pattern.ax.XLabel.String = 'azimuth';
  hd.pattern.ax.YLabel.String = 'elevation';
  
end



function hd = displayCosine(maxvelocity)

  time = 0:.2:25;
  velocity = sin(time*2*pi/10) * maxvelocity;
  
  hd.cosine.fg = figure();
  hd.cosine.ax = axes();
  hold on
  hd.cosine.ob(1) = plot(time([1 end]),[0 0]);
  hd.cosine.ob(2) = plot(time,velocity);
  hold off
  
  hd.cosine.ob(1).LineStyle = ':';
  hd.cosine.ob(1).Color = .0*[1 1 1];
  hd.cosine.ob(2).Color = .2*[1 1 1];
  hd.cosine.ax.Box = 'off';
  hd.cosine.ax.XGrid = 'on';
%   hd.cosine.ax.XMinorGrid = 'on';
%   hd.cosine.ax.YGrid = 'on';
%   hd.cosine.ax.XAxisLocation = 'origin';
  hd.cosine.ax.XTick = [0 25];
%   hd.cosine.ax.XTick = 25;
%   hd.cosine.ax.XTickLabel = hd.cosine.ax.XTick;
  hd.cosine.ax.YTick = [-1 1] * maxvelocity;

end




function [handle,data,datafit] = gainAnalysisDisk(varargin) %#ok<FNDEF>
% GAINANALYSIS and all subfunctions were written by Florian Alexander 
% Dehmelt in 2018. This is version 20180911, a.k.a. gainAnalysis20180911.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:
%   .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(end).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);
%   valid.rotation      = @(x) ischar(x);
                        
  default.resolution    = 200;
  default.camazimuth    = 0;
  default.camelevation  = 45;
  default.normalisation = '';
%   default.rotation      = 'none';
  
  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);
%   addOptional(p, 'Rotation',      default.rotation,      valid.rotation);

  parse(p,varargin{:});
  
  resolution    = p.Results.Resolution;
  cam.azimuth   = p.Results.CamAzimuth;
  cam.elevation = p.Results.CamElevation;
  normalisation = p.Results.Normalisation;
%   rotation      = p.Results.Rotation;
  
  
  % 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 = 9:46;
%   datarange = 9:44;
  % !!! IMPORTANT !!! These are the Excel row numbers, not stimulus types.
    
  
  % READ DATA
  % Query user to manually select the main data regularset. Other files have
  % standard names and are expected to be in the same folder.
  [mfilefolder,~,~] = fileparts(mfilename('fullpath'));
%   [regularset.result,regularset.folder] = uigetfile([mfilefolder,'\*.mat']);
  regularset.file =  'result_20180110_125537.mat';
  regularset.folder = ['..\..\data\Fig3 default disk data\Data D Series (neue Experimente)\'];
  rotatedset.file = 'result_20181001_114504.mat';
  rotatedset.folder = ['..\..\data\Fig3 default disk data\Control Data Rotation (arena flip)\'];
  blackset.file = 'result_20181026_111904.mat';
  blackset.folder = ['..\..\data\Fig3FigSupp4 asymmetry data\Control Data Blocked (one side black)\'];
  
  regularset.structure = load([regularset.folder,regularset.file],'result');
  rotatedset.structure = load([rotatedset.folder,rotatedset.file],'result');
  blackset.structure = load([blackset.folder,blackset.file],'result');
  regularset.result = regularset.structure.result;
  rotatedset.result = rotatedset.structure.result;
  blackset.result = blackset.structure.result;
  
  % Discard unwanted trials, based on .xlsx documenting manual assessment.
  regularset.result = discardTrial(regularset.result,regularset.folder);
  rotatedset.result = discardTrial(rotatedset.result,rotatedset.folder);
  blackset.result = discardTrial(blackset.result,blackset.folder);
  
  % Pool gain data in various ways (this is a significant subfunction!).
  regularset.gain = poolGainData(regularset.result);
  rotatedset.gain = poolGainData(rotatedset.result);
  blackset.gain = poolGainData(blackset.result);
  
  
  % READ STIMULUS INFORMATION
  regularset.stimcentre = readStimulus(regularset.folder,datarange);
  rotatedset.stimcentre = readStimulus(rotatedset.folder,datarange);
  blackset.stimcentre = readStimulus(blackset.folder,datarange);
  
%   reordering = [1;2;3;5;4;7;6]; % hemisphere stimuli, U vs. D and L vs. R flipped; now add disks:
%   for regularindex = 1:size(regularset.stimcentre,1)
%     % Find rotated-arena stimulus appearing at same position relative to fish:
%     rotatedindex = find(sum(regularset.stimcentre(regularindex,:)==rotatedset.stimcentre,2)==2);
%     reordering = [reordering; rotatedindex]; %#ok<AGROW>
%   end
%    
  % Correction above is unnecessary if stimulus coordinates were already in
  % fish-centred coordinates. In this case, do not reorder data:
  reordering = 1:45;
    
  % CHOOSE SUBSET OF DATA TO ANALYSE, AND PROCEED TO ANALYSE IT
  variantlist = {'regular','rotated','corrected'};
  for variant = 1:numel(variantlist)
    
    switch variantlist{variant}
      case 'regular'   
        data(variant).left  = regularset.gain.mean.FR{1};
        data(variant).right = regularset.gain.mean.FR{2};
%         for selection = datarange-1
%           gLreg = regularset.gain.mean.FR{1}(selection);
%           gRreg = regularset.gain.mean.FR{2}(selection);
%           data(variant).yoking.regular(selection) = (gLreg-gRreg)/(gLreg+gRreg);
%         end
      case 'rotated'
        data(variant).left  = rotatedset.gain.mean.FR{1}(reordering);
        data(variant).right = rotatedset.gain.mean.FR{1}(reordering);
%         for selection = datarange-1
%           gLrot = rotatedset.gain.mean.FR{1}(reordering(selection));
%           gRrot = rotatedset.gain.mean.FR{2}(reordering(selection));
%           data(variant).yoking.rotated(selection) = (gLrot-gRrot)/(gLrot+gRrot);
%         end
      case 'corrected'
%         data(variant).left  = (regularset.gain.mean.FR{1} + rotatedset.gain.mean.FR{1}) / 2;
%         data(variant).right = (regularset.gain.mean.FR{2} + rotatedset.gain.mean.FR{1}) / 2;
        data(variant).left  = (regularset.gain.mean.FR{1} + rotatedset.gain.mean.FR{1}(reordering)) / 2;
        data(variant).right = (regularset.gain.mean.FR{2} + rotatedset.gain.mean.FR{1}(reordering)) / 2;
        data(variant).unstimblack.left  = blackset.gain.mean.FR{1};
        data(variant).unstimblack.right = blackset.gain.mean.FR{2};
                
%         for selection = datarange-1
        for selection = 1:(max(datarange)-1)
          gLreg = regularset.gain.mean.FR{1}(selection);
          gRreg = regularset.gain.mean.FR{2}(selection);
          gLrot = rotatedset.gain.mean.FR{1}(reordering(selection));
          gRrot = rotatedset.gain.mean.FR{2}(reordering(selection));
          gLblk = blackset.gain.mean.FR{1}(selection);
          gRblk = blackset.gain.mean.FR{2}(selection);
          data(variant).asymmetry.A1(selection) = ((gLreg-gRreg)/(gLreg+gRreg) - (gLrot-gRrot)/(gLrot+gRrot)) / 2;
          data(variant).asymmetry.A2(selection) = ((gLreg-gRreg)/(gLreg+gRreg) + (gLrot-gRrot)/(gLrot+gRrot)) / 2;
          data(variant).yoking.regular(selection) = (gLreg-gRreg)/(gLreg+gRreg);
          data(variant).yoking.rotated(selection) = (gLrot-gRrot)/(gLrot+gRrot);
          data(variant).yoking.corrected(selection) = (mean([gLreg,gLrot])-mean([gRreg,gRrot])) / ...
                                                      (mean([gLreg,gLrot])+mean([gRreg,gRrot]));
          data(variant).yoking.black(selection) = (gLblk-gRblk)/(gLblk+gRblk);

        end
    end
    data(variant).both  = [data(variant).left(datarange(1:end/2)-1);data(variant).right(datarange(end/2+1:end)-1)]; 
%     data(variant).both  = mean([data(variant).left,data(variant).right],2);  % important change! no more averaging 

    % Create figure(s).
    hd.gain(variant).fg = figure();
    hd.gain(variant).fg.Name = 'OKR gain';
    hd.gain(variant).fg.Color = [1 1 1];
    hd.gain(variant).fg.Position = [50 50 1000 650];
    hd.gain(variant).fg.Colormap = flipud(cbrewer('div','RdBu',100));
  %   hd.gain(variant).fg.Colormap = cbrewer('div','PRGn',100);

    choicelist = {'left','both','right'};
    for choice = 1:numel(choicelist)

      switch choicelist{choice}
        case 'right'
          response = data(variant).right(datarange-1);  % -1 is offset by Excel title row 
        case 'left'
          response = data(variant).left(datarange-1);
        case 'both'
          response = data(variant).both;
        otherwise
          error(['Input argument must be one of the following strings: ' ...
                 '''left'', ''right'', or ''both''.'])
      end
      
      [x,y,z] = sphere(resolution);
      [datafit(choice).matrix,datafit(choice).param] = ...
        vonMisesFisherFit(regularset.stimcentre,response,x,y,z);

      % 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(regularset.stimcentre(:,1), regularset.stimcentre(:,2), ...
                ones(size(regularset.stimcentre,1),1),'CCW');

      

      % plot all of the above onto the correct panel of the existing figure
      panel = choice;

      figure(hd.gain(variant).fg.Number)
      hd.gain(variant).ax(panel) = axes;
      hd.gain(variant).ax(panel).Position = [.13 .45 .35 .45] + (panel-1)*[.25 0 0 0];
      hold on
      hd.gain(variant).ob(panel,1) = surf(x,y,z,datafit(panel).matrix);
      shading interp
      colour.eye = 0*[1 1 1];
      colour.body = 1*[1 1 1];
      hd.gain(variant).ob(panel,2:4) = plot3dFish(scale,yaw,pitch,roll,colour);
      hd.gain(variant).ob(panel,5) = plot3(equator.x,equator.y,equator.z,'k','LineWidth',2);
      hd.gain(variant).ob(panel,6) = plot3(meridian.x,meridian.y,meridian.z,'--k');
      hd.gain(variant).ob(panel,7) = plot3(aura.x,aura.y,aura.z,':k');
      hd.gain(variant).ob(panel,8) = scatter3(1.03*sample.x,1.03*sample.y,1.03*sample.z, ...
                                180*ones(size(sample.x)),response,'Marker','o');
      hold off

    end

    % Apply full formatting to the figure(s). Subfunction included below.
    hd = fishbowlFigure(hd,cam,choicelist);

%     % Compute yoking/lateralisation index per individual stimulus location.
%     convention = 'CCW';
%     switch convention
%       case 'CCW'
%         stim.isleft  = stim.azimuth > 0;
%         stim.isright = stim.azimuth < 0;
%       case 'CW'
%         stim.isleft  = stim.azimuth < 0;
%         stim.isright = stim.azimuth > 0;
%       otherwise
%         error('Sign convention must either be ''CCW'' or ''CW''.')
%     end
  
  end
     
  % Collect output arguments, unless already done (e.g., "datafit").
  handle = hd;
  
end



function stimcentre = readStimulus(datafolder,datarange)

  stimulus = 'Stimulus.xlsx';
  [~,tableconvention,~] = xlsread([datafolder,stimulus],'C1:C1');
  azimuth = xlsread([datafolder,stimulus],['C',num2str(datarange(1)),':C',num2str(datarange(end))]);
  elevation = xlsread([datafolder,stimulus],['D',num2str(datarange(1)),':D',num2str(datarange(end))]);

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

  stimcentre = [azimuth,elevation];

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 [fitresult,fitparam] = vonMisesFisherFit(stimcentre,response,x,y,z)

  predictor = stimcentre;

  % Initialise fit parameters.
  offset = 0;
  amplitude(1) = 1;
  amplitude(2) = 1;
  meanazimuth(1)   =  70;
  meanelevation(1) =  20; % azimuth, elevation
  meanazimuth(2)   = -70;
  meanelevation(2) =  20; % azimuth, elevation
  concentration(1) = 2.5;
  concentration(2) = 2.5;

  initialparam = [offset,        ...
                  amplitude(1),     ...
                  meanazimuth(1),   ...
                  meanelevation(1), ...
                  concentration(1), ...
                  amplitude(2),     ...
                  meanazimuth(2),   ...
                  meanelevation(2), ...
                  concentration(2)];

  % Fit.
  option.MaxIter = 1e4;
  datafit.param = nlinfit(predictor,response,@vonmisesfisher9param,initialparam,option);
  datafit.value = vonmisesfisher9param(datafit.param,predictor);

  % Identify centres.
  azimuth(1)   = datafit.param(3);
  elevation(1) = datafit.param(4);
  azimuth(2)   = datafit.param(7);
  elevation(2) = datafit.param(8);
  for k = 1:2

%       azimuth(k)          = datafit.param(3); % This was incorrect - same value read twice.
%       elevation(k)        = datafit.param(4);
    uniqueazimuth(k)    = mod(azimuth(k)+180,360) - 180;
    uniqueelevation(k)  = mod(elevation(k)+90,180) - 90;
    datafit.centre(:,k) = [uniqueazimuth(k), uniqueelevation(k)];

  end

  % compute geographic coordinates of sphere, and its colour data values
  [azimuth,elevation,radius] = car2geo(x,y,z,'CCW');
  predictor = [azimuth(:),elevation(:),radius(:)];
  fitresult = reshape(vonmisesfisher9param(datafit.param,predictor),size(radius));
  fitparam = datafit.param;

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.
    
  cmin = min([min(min(hd.gain(end).ob(1,1).CData)), ...
              min(min(hd.gain(end).ob(2,1).CData)), ...
              min(min(hd.gain(end).ob(3,1).CData)), ...
              min(min(hd.gain(end).ob(1,8).CData)), ...
              min(min(hd.gain(end).ob(2,8).CData)), ...
              min(min(hd.gain(end).ob(3,8).CData))]);
  cmax = max([max(max(hd.gain(end).ob(1,1).CData)), ...
              max(max(hd.gain(end).ob(2,1).CData)), ...
              max(max(hd.gain(end).ob(3,1).CData)), ...
              max(max(hd.gain(end).ob(1,8).CData)), ...
              max(max(hd.gain(end).ob(2,8).CData)), ...
              max(max(hd.gain(end).ob(3,8).CData))]);
  clim = max(abs([cmin cmax]))*[-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 = 'merge stimulated';
%         string.title = 'mean = (g_L+ g_R)/2';
    end
       
    hd.gain(end).ax(panel).Title.String = string.title;
%     title(string.title,'Interpreter','LaTeX')
    hd.gain(end).ob(panel,8).MarkerFaceColor = hd.gain(end).ob(panel,8).MarkerEdgeColor;
    hd.gain(end).ob(panel,8).MarkerEdgeColor = [0 0 0];
    hd.gain(end).ax(panel).CLim = [0 cmax];
%     hd.gain(end).ax(panel).CLim = [0 max(max(datafit(choice).matrix))];
%     hd.gain(end).ob(panel,1).FaceAlpha = .85;
    hd.gain(end).ob(panel,1).FaceAlpha = .75;
    hd.gain(end).ob(panel,1).EdgeColor = 'none';
%     hd.gain(end).ob(panel,8).Shading = 'flat';
    shading interp
    
    axis square
    xlabel('x')
    ylabel('y')
    zlabel('z')            
    hd.gain(end).ax(panel).XLim = [-1.1 1.1];
    hd.gain(end).ax(panel).YLim = [-1.1 1.1];
    hd.gain(end).ax(panel).ZLim = [-1.1 1.1];
%     hd.gain(end).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(end).ax(panel).CameraPosition = [cam.x cam.y cam.z]; % FD20190507: this hides 3rd fish!!!
    hd.gain(end).ax(panel).CameraTarget = [0 0 0];
    hd.gain(end).ax(panel).CameraUpVector = [0 0 1];
%     hd.gain(end).ax(panel).CameraViewAngle = 10.134;
    hd.gain(end).ax(panel).CameraViewAngle = 9;

  end
  
  hd.gain(end).ax(4) = axes;
  hd.gain(end).ax(4).Position = hd.gain(end).ax(1).Position - [.13 -.1 .2 .2];
  hd.gain(end).ax(4).CLim = hd.gain(end).ax(1).CLim;
  hd.gain(end).ob(4,1) = colorbar;
  hd.gain(end).ob(4,1).Label.String = 'OKR gain';
  hd.gain(end).ob(4,1).Label.Units = 'normalized';
  hd.gain(end).ob(4,1).Label.Position = [-.8 1 1] .* hd.gain(end).ob(4,1).Label.Position;
  
  [hd.gain(end).ob(4,1).Label.FontSize, hd.gain(end).ob(4,1).FontSize, ...
   hd.gain(end).ax.FontSize] = deal(14);
 
  [hd.gain(end).ax(1).Title.FontSize, ...
   hd.gain(end).ax(2).Title.FontSize, ...
   hd.gain(end).ax(3).Title.FontSize] = deal(14);
   
  [hd.gain(end).ax(1).Title.FontWeight, ...
   hd.gain(end).ax(2).Title.FontWeight, ...
   hd.gain(end).ax(3).Title.FontWeight] = deal('normal');

  
  [hd.gain(end).ax.Visible] = deal('off');
%   [hd.gain(end).ax.Title.Visible] = deal('on');
  set(([hd.gain(end).ax.Title]),'Visible','on')
  set(([hd.gain(end).ax.Title]),'Position',[1 1 .7] .* hd.gain(end).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\';
  regularset.badtrial = 'DiscardTrial.xlsx';
  badtrial = xlsread([folder.badtrial,regularset.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\';
  regularset.badtrial = 'ZeroTrial.xlsx';
  badtrial = xlsread([folder.badtrial,regularset.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,colour)

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

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

  resolution = 10;
%   colour.eye = .1*[1 1 1];
%   colour.body = .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(colour.eye,[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(colour.body,[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,stim] = displayDiskStimuli(choice)
  
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
  
  regularset.folder = '..\..\data\Fig3 default disk data\Data D Series (neue Experimente)\';
  datarange = 9:46;
  
%   regularset.folder = '..\..\data\Fig5 size data\Data A Series (A-Series)\';
%   datarange = 2:43;
%   datarange = 2:8;

%   regularset.folder = '..\..\data\Fig5 frequency data\Data V Series (V-Serie)\';
%   datarange = 2:43;
  
  regularset.stimulus = 'Stimulus.xlsx';
  [~,stim.tableconvention,~] = xlsread([regularset.folder,regularset.stimulus],'C1:C1');
  stim.azimuth = xlsread([regularset.folder,regularset.stimulus],['C',num2str(datarange(1)),':C',num2str(datarange(end))]);
  stim.elevation = xlsread([regularset.folder,regularset.stimulus],['D',num2str(datarange(1)),':D',num2str(datarange(end))]);
%   stim.radius = xlsread([regularset.folder,regularset.stimulus],['E',num2str(datarange(1)),':E',num2str(datarange(end))]);
% stim.azimuth = stim.azimuth(1:2:end);
% stim.elevation = stim.elevation(1:2:end);
  stim.radius = 40*ones(size(stim.elevation));

  
% %   [surface.]backgroundSphere
%   % compute background sphere with shading
%   [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+1);
% %   surface.c = 2.^(surface.c+1);
  
  
  % 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 = .2;
  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
  [diskcentre.x,diskcentre.y,diskcentre.z] = ...
    geo2car(diskcentre.azim,diskcentre.elev,diskcentre.rads,'CW');
  [textcentre.x,textcentre.y,textcentre.z] = ...
    geo2car(diskcentre.azim,.91*diskcentre.elev,1.10*diskcentre.rads,'CW');  
    
  % create figure
  hd.default.fg = figure();
  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;
%   hd.default.globe = surf(surface.x,surface.y,surface.z,surface.c);
%   shading flat
%   axis square
%   greyscale = (.8:.01:.99)' * [1 1 1];  
%   colormap(greyscale)
  fish.scale = .05;
  fish.yaw = 180;
  fish.pitch = 0;
  fish.roll = 0;
  colour.eye = 0*[1 1 1];
  colour.body = .5*[1 1 1];
  plot3dFish(fish.scale,fish.yaw,fish.pitch,fish.roll,colour);
  hd.default.globe.FaceAlpha = .5;
  
  
  cropdisk = 8;
   switch choice
    case 'left'
      disklist = 1 : numdisk/2;
    case 'right'
      disklist = numdisk/2+1 : numdisk;
    case 'single'
      disklist = cropdisk;
  end
  
  % plot disk boundaries
%   diskcolour = cbrewer('seq','Greens',numdisk*2);
  cbrewergreen = cbrewer('seq','Greens',numdisk);
  diskcolour = repmat(cbrewergreen(round(numdisk/1.4),:),[numdisk,1]);
%   diskcolour = flipud(cbrewer('seq','Greys',numdisk*2));
%   diskcolour = diskcolour(randperm(numdisk),:);
  for disk = disklist
    if sum(disk==cropdisk)
      plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',2,'Color',0*[1 1 1])
      disktext(disk) = text(textcentre.x(disk),textcentre.y(disk),textcentre.z(disk), ...
        ['',num2str(disk)], ...
        'FontSize', 13, 'HorizontalAlignment','center','VerticalAlignment','middle');
    else
        plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',1,'Color',.6*[1 1 1])
      if sum(disk==disklist)
      disktext(disk) = text(textcentre.x(disk),textcentre.y(disk),textcentre.z(disk), ...
        ['',num2str(disk)], ...
        'FontSize', 13, 'HorizontalAlignment','center','VerticalAlignment','middle');
      end
%     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])
    end
  end
  
  hold off
  
  switch choice
    
    case 'single'

      disktext(cropdisk).Visible = 'off';
      
      barwidth = 360/22/2;
    %   gratingazimuth = -180:barwidth:(180-barwidth);
      gratingazimuth = -180:barwidth:-80;
      for bar = 1:2:numel(gratingazimuth)-1
    %     if bar < numel(gratingazimuth)
          azimlimit = gratingazimuth([bar bar+1]);
    %     else
    %       azimlimit = gratingazimuth([bar 1]);
    %     end
        displaySphereBar(azimlimit(1),azimlimit(2),boundary.azim{cropdisk},boundary.elev{cropdisk})
      end
      
    case {'left','right'}

      disktext(disklist(ceil(end/2))).Position = ...
        disktext(disklist(ceil(end/2))).Position - [0 0 .16];
      
  end
    

%   hd.default.ax(1).CameraPosition = [2.3 17 2.2];

  switch choice
    case {'left','single'}
      hd.default.ax(1).CameraPosition = [0 20 0];
    case 'right'
      hd.default.ax(1).CameraPosition = [0 -20 0];
  end
  hd.default.ax(1).CameraViewAngle = 6;
  [hd.default.ax(1).XAxis.Visible, ...
   hd.default.ax(1).YAxis.Visible, ...
   hd.default.ax(1).ZAxis.Visible] = deal('off');
  
  
end



function globe = displaySphere
  
  [surface.x,surface.y,surface.z] = sphere(300);
  
  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 interp
  axis square

  greyscale = (.8:.01:.99)' * [1 1 1];  
%   greyscale = (.9:.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+1e-9 1-1e-9];
  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)-168;
  
  azimuthOut   = border.azimuth';
  elevationOut = border.elevation';

end



function hd = displaySphereBar(azimlimit1,azimlimit2,boundaryazim,boundaryelev)

  resolution = .01;
  if azimlimit2 < azimlimit1
  error('sdasdsa')
  end
  
  crop.index = logical((azimlimit1 <= boundaryazim) .* (boundaryazim <= azimlimit2));
  
  if sum(crop.index) > 0
    
    crop.azim = boundaryazim(crop.index)';
    crop.elev = boundaryelev(crop.index)';

    crop.azim1 = crop.azim(diff(crop.azim) > 0);
    crop.azim2 = crop.azim(diff(crop.azim) < 0);
    crop.elev1 = crop.elev(diff(crop.azim) > 0);
    crop.elev2 = crop.elev(diff(crop.azim) < 0);

    sequence = @(a,b) a:(sign(b-a)*resolution):b;
    edge.elev1 = sequence(crop.elev1(1), crop.elev2(end));
    edge.elev2 = sequence(crop.elev2(1), crop.elev1(end));
%     edge.elev1 = sequence(crop.elev1(1), crop.elev2(end));
%     edge.elev2 = sequence(crop.elev2(1), crop.elev1(end));
    edge.azim1 = azimlimit1*ones(size(edge.elev1));
    edge.azim2 = azimlimit2*ones(size(edge.elev2));
    
    azim = [edge.azim1,fliplr(crop.azim2),edge.azim2,fliplr(crop.azim1)];
    elev = [edge.elev1,fliplr(crop.elev2),edge.elev2,fliplr(crop.elev1)];
    rads = 1.009*ones(size(elev));

    [bar.x,bar.y,bar.z] = geo2car(azim,elev,rads,'CW');
    hold on, p = patch(bar.x,bar.y,bar.z,'k'); hold off
    % p.FaceColor = [.2 .8 .4];
    p.FaceColor = .2*[1 1 1];
    p.EdgeColor = 'none';
%     hold on, plot3(bar.x,bar.y,bar.z,'r'); hold off

  end

end