%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, main program (control)
% dice_experiment_WD(1,6, 10, 1000, 1, 0,0,"tmp")
function dice_experiment_WD(
  N_dots_min,          % minimal nb of dots ("Augen") of dice
  N_dots_max,          % maximal nb of dots ("Augen") of dice
  N_dice,              % nb of dice to sum up
  N_exp,               % nb of repetitions of experiment
  N_animate,           % nb of frames for animation (if set to 1: no animation)
  ref_curve_en,        % plot Gauss function of same mean and variance value as reference
  occ_prob_en,         % y-axis frequ. of occurrence, or prob.
  svg_en,
  webdemo_en,                   % if 0, then test-mode in Octave/Matlab commandline, no webserver running
  filename                      % output file for webdemo
)


  landscape = "-S930,350";
  portrait  = "-S640,480";
  square    = "-S480,480";
  tall  = "-S320,480";

  output_format = portrait;

  if (webdemo_en==1)
    fh = figure(1);
    set(fh, "visible", "off");
  end

  set(0, "defaultlinelinewidth", 2);

  if (svg_en)
    set(0, "defaultaxesfontsize", 6);
    set(0, "defaultaxesfontname", "Arial");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Arial");
    const_lw = 1;
  else
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Arial");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Arial");
    const_lw = 2;
  end


  % RGB color definitions, http://www.rapidtables.com/web/color/RGB_Color.htm; from dark to light
  gray =    [[0,0,0]; [32,32,32]; [64,64,64]; [96,96,96]; [128,128,128]; [160,160,160]; [192,192,192]; [224,224,224]; [255,255,255]]/255;
  violett = [[51,0,25]; [102,0,51]; [153,0,76]; [204,0,102]; [255,0,127]; [255,51,153]; [255,102,178]; [255,153,204]; [255,204,229]]/255;
  magenta = [[51,0,51]; [102,0,102]; [153,0,153]; [204,0,204]; [255,0,255]; [255,51,255]; [255,102,255]; [255,153,255]; [255,204,255]]/255;
  lilac   = [[25,0,51]; [51,0,102]; [76,0,153]; [102,0,204]; [127,0,255]; [153,51,255]; [178,102,255]; [204,153,255]; [229,204,255]]/255;
  blue    = [[0,0,51]; [0,0,102]; [0,0,153]; [0,0,204]; [0,0,255]; [51,51,255]; [102,102,255]; [153,153,255]; [204,204,255]]/255;
  seablue = [[0,25,51]; [0,51,102]; [0,76,153]; [0,102,204]; [0,128,255]; [51,153,255]; [102,178,255]; [153,204,255]; [204,229,255]]/255;
  cyan    = [[0,51,51]; [0,102,102]; [0,153,153]; [0,204,204]; [0,255,255]; [51,255,255]; [102,255,255]; [153,255,255]; [204,255,255]]/255;
  seagreen = [[0,51,25]; [0,102,51]; [0,153,76]; [0,204,102]; [0,255,128]; [51,255,153]; [102,255,178]; [153,255,204]; [204,255,229]]/255;
  green   = [[0,51,0]; [0,102,0]; [0,153,0]; [0,204,0]; [0,255,0]; [51,255,51]; [102,255,102]; [153,255,153]; [204,255,204]]/255;
  lightgreen = [[25,51,0]; [51,102,0]; [76,153,0]; [102,204,0]; [128,255,0]; [153,255,51]; [178,255,102]; [204,255,153]; [229,255,204]]/255;
  yellow  = [[51,51,0]; [102,102,0]; [153,153,0]; [204,204,0]; [255,255,0]; [255,255,51]; [255,255,102]; [255,255,153]; [255,255,204]]/255;
  amber   = [[51,25,0]; [102,51,0]; [153,76,0]; [204,102,0]; [255,128,0]; [255,153,51]; [255,178,102]; [255,204,153]; [255,229,204]]/255;
  red     = [[51,0,0]; [102,0,0]; [153,0,0]; [204,0,0]; [255,0,0]; [255,51,51]; [255,102,102]; [255,153,153]; [255,204,204]]/255;


  N_outcomes=N_dots_max-N_dots_min+1;

  N_exp_frame=N_exp/N_animate;


  xmax=N_dice*N_dots_max;
  x_acc=zeros(1,xmax);

%% animation loop
for (frame_nb=1:N_animate)

  clf;
  hold on

  title([ "nach ", num2str(N_exp_frame*frame_nb), " Wuerfen"]);

  xlabel('Ausgang des Experiments x_i');

  if (occ_prob_en==0)
    ylabel('Auftrittshaeufigkeit fuer X=x_i');
  else
    ylabel('Auftrittswahrscheinlichkeit P[X=x_i]');
  end



  [x,x_acc,xmin,xmax]=experiment_sum_dice(x_acc, N_dots_min, N_dots_max, N_dice, N_exp_frame);
  [mean_val,variance]=compute_mean_variance(x,x_acc./(N_exp_frame*frame_nb));


  if (occ_prob_en==1) % prob. rather than frequ. of occ.
    x_acc_plot = x_acc./(N_exp_frame*frame_nb);
  else
    x_acc_plot = x_acc;
  end


  xlim([xmin-0.2,xmax+0.2]);

  ymax=max(x_acc_plot)*1.2;
  ylim([0,ymax]);
  plot(x, x_acc_plot ,"linewidth", const_lw+1, "o", "color", "b");
%  bar(x, x_acc_plot ,"linewidth", const_lw+1, "color", "b");  % maybe later...



  % reference curve
  if (ref_curve_en==1)
    xx=linspace(x(1),x(end),100);
    if (N_dice==1)
      uniform_ref = (N_exp_frame*frame_nb)/N_outcomes;
      if (occ_prob_en==1)
        uniform_ref = uniform_ref ./ (N_exp_frame*frame_nb);
      end
      plot(xx, ones(1,length(xx))*uniform_ref,"linewidth", const_lw+1, "-", "color", "m");
    else
      gauss_ref = gauss(xx-mean_val,variance);
      if (occ_prob_en==0)
        gauss_ref = gauss_ref.*(N_exp_frame*frame_nb);
      end
      plot(xx, gauss_ref ,"linewidth", const_lw+1, "-", "color", "m");
    end
  end
  legend ("Wuerfelergebnisse", "Referenzkurve");


  hold off
  
  if (N_animate>1)
    newFilename = [filename '_' num2str(frame_nb-1)];
  else
    newFilename = filename;
  end

  print(newFilename, '-dpng',output_format);
end % for, animation loop

end

% END, main program (control)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




function [x,x_acc,xmin,xmax]=experiment_sum_dice(x_acc, N_dots_min, N_dots_max, N_dice, N_exp);

  xmin=N_dice*N_dots_min;
  xmax=N_dice*N_dots_max;
  N_dots_range=N_dots_max-N_dots_min+1;

  for i=1:N_exp
    acc=0;
    for k=1:N_dice
      rv=ceil(rand(1,1)*N_dots_range)+N_dots_min-1;
      acc=acc+rv;
    end
    x_acc(acc)=x_acc(acc)+1;
  end

  x=(1:1:xmax);
end



function [mean_val,variance]=compute_mean_variance(x,x_acc)
  N=length(x);

  m=0;
  for i=1:N
    m=m+x(i)*x_acc(i);
  end

  mean_val=m;

  v=0;
  for i=1:N
    v=v+(x(i)-m)^2*x_acc(i);
  end

  variance=v;

end








function y = gauss(x,sn2)
  y = 1/sqrt(2*pi*sn2)*exp(-1/2/sn2.*x.^2);
end


function y = Q(x)
  y = 0.5*(1-erf(x/sqrt(2)));
end


function y = Qinv(x)
  y = sqrt(2)*erfinv(1-2*x);
end


function y = PHI(x)
  y = Q(-x);
end




