function varargout = chaotic_scattering_gui(varargin)
%CHAOTIC_SCATTERING_GUI M-file for chaotic_scattering_gui.fig
%      CHAOTIC_SCATTERING_GUI, by itself, creates a new CHAOTIC_SCATTERING_GUI or raises the existing
%      singleton*.
%
%      H = CHAOTIC_SCATTERING_GUI returns the handle to a new CHAOTIC_SCATTERING_GUI or the handle to
%      the existing singleton*.
%
%      CHAOTIC_SCATTERING_GUI('Property','Value',...) creates a new CHAOTIC_SCATTERING_GUI using the
%      given property value pairs. Unrecognized properties are passed via
%      varargin to chaotic_scattering_gui_OpeningFcn.  This calling syntax produces a
%      warning when there is an existing singleton*.
%
%      CHAOTIC_SCATTERING_GUI('CALLBACK') and CHAOTIC_SCATTERING_GUI('CALLBACK',hObject,...) call the
%      local function named CALLBACK in CHAOTIC_SCATTERING_GUI.M with the given input
%      arguments.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help chaotic_scattering_gui

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @chaotic_scattering_gui_OpeningFcn, ...
                   'gui_OutputFcn',  @chaotic_scattering_gui_OutputFcn, ...
                   'gui_LayoutFcn',  @chaotic_scattering_gui_LayoutFcn, ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
   gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before chaotic_scattering_gui is made visible.
function chaotic_scattering_gui_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   unrecognized PropertyName/PropertyValue pairs from the
%            command line (see VARARGIN)

% Choose default command line output for chaotic_scattering_gui
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes chaotic_scattering_gui wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% Init graphics
initGraphics();

% --- Outputs from this function are returned to the command line.
function varargout = chaotic_scattering_gui_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in togglebutton3.
function togglebutton3_Callback(hObject, eventdata, handles)
% hObject    handle to togglebutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of togglebutton3

% Computing fractal dimension

% Show input dialog
prompt={'Number of points for the curve','Number of trials for a given epsilon'};
name = '';
defaultans = {'100','1000'};
answer = inputdlg(prompt,name,[1 30],defaultans);
drawnow;
numPoints = str2num(answer{1});
numTrials = str2num(answer{2});

% Create progress bar
h = waitbar(0,'Computing path-length array...','CreateCancelBtn',...
            'setappdata(gcbf,''canceling'',1)');
setappdata(h,'canceling',0);

% Change bar color
set(findobj(h,'type','patch'), ...
'edgecolor','b','facecolor','b');

% Set perturbation array
epsilon = logspace(-2,-8,numPoints);

% Number of uncertain initial conditions
numUncertain = 0;

% Cancel flag for waitbar
cancel_waitbar = 0;

% Loop over different epsilon perturbation
for i=1:numPoints
   % Generating points
   pos = -2+4*rand(numTrials,1);
   % Update progress bar
   waitbar(i/numPoints, h, sprintf('%d %%', floor(100*i/numPoints)));
   % Check for Cancel button press
   if getappdata(h,'canceling')
     cancel_waitbar = 1;
     break;
   end
   % Compute first initial condition
   numUncertain(i) = 0;
   % Loop over all random positions
   for j=1:numTrials
      % Compute first initial condition
      [theta1(j), sum_path1(j) numBounces1(j)] = compute_path(pos(j),'nodraw');
      % Compute same initial condition with perturbation
      [theta2(j), sum_path2(j) numBounces2(j)] = compute_path(pos(j)+epsilon(i),'nodraw');
      if (numBounces1(j) ~= numBounces2(j))
        numUncertain(i)=numUncertain(i)+1;
      end
   end
   functionFraction(i) = numUncertain(i)/numTrials;
end

% Delete progress bar
delete(h);

% Return if cancel button pressed
if cancel_waitbar
  return;
end

% Remove Nan values with Log10 array for functionFraction
indexNotInf = ~isinf(log10(functionFraction));
log10_epsilon = log10(epsilon(indexNotInf));
log10_function = log10(functionFraction(indexNotInf));

if ((range(log10_epsilon) == 0) || (range(log10_function) == 0) ...
   || (numel(log10_epsilon) < 2) || (numel(log10_function) < 2))
  errordlg('Not enough points : try with a higher number of points or trials','Error');
  return;
end

% Plot log10_epsilon vs log10_function
figure(3);
h1 = plot(log10_epsilon,log10_function);
xlim([min(log10_epsilon) max(log10_epsilon)]);
ylim([min(log10_function) max(log10_function)]);

% Save by printing into file
dlmwrite('UncertaintyFunction.dat',[epsilon;functionFraction]','delimiter',' ');

% 15-digit floating point
format long e;

% Linear regression
beta = polyfit(log10_epsilon,log10_function,1);
if (isinf(beta(1)) || isinf(beta(2)))
  errordlg('Regression could not be done','Error');
  return;
end

% Computing standard error on slope coefficient
stats=regstats(log10_function,log10_epsilon,'linear',{'beta','covb'});
error_coefficients=sqrt(diag(stats.covb));
error_slope=error_coefficients(2);

% Plot linear regresion
hold on;
log10_yFit = polyval(beta,log10_epsilon);
h2 = plot(log10_epsilon,log10_yFit,'r');
title('Computing uncertainty coefficient');
xlabel('$log_{10}(\epsilon)$','interpreter','latex');
ylabel('$log_{10}f(\epsilon)$','interpreter','latex');

% legends
legend([h1 h2],'Fraction function','Linear regression');
text(min(log10_epsilon)+0.1*(max(log10_epsilon)-min(log10_epsilon)),...
     min(log10_function)+0.9*(max(log10_function-min(log10_function))),sprintf('\\alpha = %f \\pm %f',beta(1),error_slope));

% Release button3
set(handles.togglebutton3,'Value',0); 


function edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double


% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function edit4_Callback(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit4 as text
%        str2double(get(hObject,'String')) returns contents of edit4 as a double


% --- Executes during object creation, after setting all properties.
function edit4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in togglebutton2.
function togglebutton2_Callback(hObject, eventdata, handles)
% hObject    handle to togglebutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of togglebutton2

% Get input ordinate
input_ordinate_start = str2double(get(handles.edit3,'String'));
input_ordinate_final = str2double(get(handles.edit4,'String'));

% Warning if input_ordinate not in range
if (input_ordinate_start >= 2.25 || input_ordinate_start <= -2.25 || input_ordinate_final >= 2.25 || input_ordinate_final <= -2.25)
  errordlg('Input must be < 2.25 and > -2.25','Error');
  return;
end

% Warning if input_ordinate_start > input_ordinate_final
if (input_ordinate_start >= input_ordinate_final)
  errordlg('Start input must be < Final input');
  return;
end

% Show input dialog
prompt={'Number of points'};
name = '';
defaultans = {'1000'};
answer = inputdlg(prompt,name,[1 30],defaultans);
drawnow;
numPoints = str2num(answer{1});

% Create progress bar
h = waitbar(0,'Computing path-length array...','CreateCancelBtn',...
            'setappdata(gcbf,''canceling'',1)');
setappdata(h,'canceling',0);

% Change bar color
set(findobj(h,'type','patch'), ...
'edgecolor','b','facecolor','b');

% Main loop
for i=1:numPoints
  % Check for Cancel button press
  if getappdata(h,'canceling')
    break;
  end
  waitbar(i/numPoints, h, sprintf('%d %%', floor(100*i/numPoints)));
  input_ordinate(i)=input_ordinate_start+(input_ordinate_final-input_ordinate_start)*(i-1)/numPoints;
  % Compute and draw path
  [theta(i), sum_path(i) numBounces(i)]= compute_path(input_ordinate(i),'nodraw');
end

% Delete progress bar
delete(h);

% Plot theta versus ordinate input %%
figure(1);
plot(input_ordinate,theta);
xlim([input_ordinate(1) input_ordinate(numel(input_ordinate))]);
ylim([0 2*pi]);

title('Theta versus Input ordinate');
xlabel('Input ordinate');
ylabel('Output theta in radian');

% Plot length versus ordinate input %%
figure(2);
plot(input_ordinate,sum_path);
xlim([input_ordinate(1) input_ordinate(numel(input_ordinate))]);
ylim([min(sum_path) max(sum_path)]);
title('Length path versus Input ordinate');
xlabel('Input ordinate');
ylabel('Output path length');

% Release button1
set(handles.togglebutton2,'Value',0); 

function edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double


% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double


% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in togglebutton1.
function togglebutton1_Callback(hObject, eventdata, handles)
% hObject    handle to togglebutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of togglebutton1

% Get input ordinate
input_ordinate = str2double(get(handles.edit1,'String'));

% Warning if input_ordinate not in range
if (input_ordinate >= 2.25 || input_ordinate <= -2.25)
  errordlg('Input must be < 2.25 and > -2.25','Error');
  return;
end

% Init graphics
cla(handles.figure1);
initGraphics();

% Compute and draw path
[output_angle sum_path numBounces] = compute_path(input_ordinate,'draw');

% Display output angle
set(handles.edit2,'string',num2str(output_angle));

% Release button1
set(handles.togglebutton1,'Value',0); 

% --- Creates and returns a handle to the GUI figure. 
function h1 = chaotic_scattering_gui_LayoutFcn(policy)
% policy - create a new figure or use a singleton. 'new' or 'reuse'.

persistent hsingleton;
if strcmpi(policy, 'reuse') & ishandle(hsingleton)
    h1 = hsingleton;
    return;
end

appdata = [];
appdata.GUIDEOptions = struct(...
    'active_h', [], ...
    'taginfo', struct(...
    'figure', 2, ...
    'axes', 2, ...
    'uipanel', 5, ...
    'edit', 9, ...
    'text', 9, ...
    'togglebutton', 6), ...
    'override', 1, ...
    'release', 13, ...
    'resize', 'none', ...
    'accessibility', 'callback', ...
    'mfile', 1, ...
    'callbacks', 1, ...
    'singleton', 1, ...
    'syscolorfig', 1, ...
    'blocking', 0);

appdata.lastValidTag = 'figure1';
appdata.GUIDELayoutEditor = [];
appdata.initTags = struct(...
    'handle', [], ...
    'tag', 'figure1');


% Force position for main window
pos_screen = get(0,'ScreenSize');
pos_window_x = pos_screen(3)/4;
pos_window_y = pos_screen(4)/4;

h1 = figure(...
'Color',[0.701960784313725 0.701960784313725 0.701960784313725],...
'Colormap',[0 0 0.5625;0 0 0.625;0 0 0.6875;0 0 0.75;0 0 0.8125;0 0 0.875;0 0 0.9375;0 0 1;0 0.0625 1;0 0.125 1;0 0.1875 1;0 0.25 1;0 0.3125 1;0 0.375 1;0 0.4375 1;0 0.5 1;0 0.5625 1;0 0.625 1;0 0.6875 1;0 0.75 1;0 0.8125 1;0 0.875 1;0 0.9375 1;0 1 1;0.0625 1 1;0.125 1 0.9375;0.1875 1 0.875;0.25 1 0.8125;0.3125 1 0.75;0.375 1 0.6875;0.4375 1 0.625;0.5 1 0.5625;0.5625 1 0.5;0.625 1 0.4375;0.6875 1 0.375;0.75 1 0.3125;0.8125 1 0.25;0.875 1 0.1875;0.9375 1 0.125;1 1 0.0625;1 1 0;1 0.9375 0;1 0.875 0;1 0.8125 0;1 0.75 0;1 0.6875 0;1 0.625 0;1 0.5625 0;1 0.5 0;1 0.4375 0;1 0.375 0;1 0.3125 0;1 0.25 0;1 0.1875 0;1 0.125 0;1 0.0625 0;1 0 0;0.9375 0 0;0.875 0 0;0.8125 0 0;0.75 0 0;0.6875 0 0;0.625 0 0;0.5625 0 0],...
'IntegerHandle','off',...
'InvertHardcopy',get(0,'defaultfigureInvertHardcopy'),...
'MenuBar','none',...
'Name','Chaotic scattering in Gaspard-Rice''s system',...
'NumberTitle','off',...
'PaperPosition',[0.0294117647058824 0.227272727272727 0.941176470588235 0.545454545454545],...
'PaperSize',[1 1],...
'Position',[pos_window_x pos_window_y 951 649],...
'Resize','off',...
'HandleVisibility','callback',...
'UserData',[],...
'Tag','figure1',...
'Visible','on',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'axes1';

h2 = axes(...
'Parent',h1,...
'Units','pixels',...
'FontUnits','normalized',...
'Position',[350 53 551 551],...
'CameraPosition',[0.5 0 41.1201920231798],...
'CameraPositionMode',get(0,'defaultaxesCameraPositionMode'),...
'Color',get(0,'defaultaxesColor'),...
'ColorOrder',get(0,'defaultaxesColorOrder'),...
'DataAspectRatio',get(0,'defaultaxesDataAspectRatio'),...
'DataAspectRatioMode','manual',...
'FontSize',0.0310945273631841,...
'LooseInset',[168.588720291759 76.3148937857383 123.199449443978 52.0328821266397],...
'XColor',get(0,'defaultaxesXColor'),...
'XLim',[0 1],...
'XLimMode','manual',...
'YColor',get(0,'defaultaxesYColor'),...
'YLim',get(0,'defaultaxesYLim'),...
'YLimMode','manual',...
'ZColor',get(0,'defaultaxesZColor'),...
'Tag','axes1',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.SerializedAnnotationV7 = struct(...
    'LegendInformation', struct(...
    'IconDisplayStyle', 'on'));

h3 = get(h2,'title');

set(h3,...
'Parent',h2,...
'Units','data',...
'FontUnits','points',...
'BackgroundColor','none',...
'Color',[0 0 0],...
'DisplayName',blanks(0),...
'EdgeColor','none',...
'DVIMode','auto',...
'FontAngle','normal',...
'FontName','Helvetica',...
'FontSize',10,...
'FontWeight','normal',...
'HorizontalAlignment','center',...
'LineStyle','-',...
'LineWidth',0.5,...
'Margin',2,...
'Position',[0.485480943738658 4.13793103448276 1.00024509545847],...
'Rotation',0,...
'String',blanks(0),...
'Interpreter','tex',...
'VerticalAlignment','bottom',...
'ButtonDownFcn',[],...
'CreateFcn', {@local_CreateFcn, [], appdata} ,...
'DeleteFcn',[],...
'BusyAction','queue',...
'HandleVisibility','off',...
'HelpTopicKey',blanks(0),...
'HitTest','on',...
'Interruptible','on',...
'SelectionHighlight','on',...
'Serializable','on',...
'Tag',blanks(0),...
'UserData',[],...
'Visible','on',...
'XLimInclude','on',...
'YLimInclude','on',...
'ZLimInclude','on',...
'CLimInclude','on',...
'ALimInclude','on',...
'IncludeRenderer','on',...
'Clipping','off');

appdata = [];
appdata.SerializedAnnotationV7 = struct(...
    'LegendInformation', struct(...
    'IconDisplayStyle', 'on'));

h4 = get(h2,'xlabel');

set(h4,...
'Parent',h2,...
'Units','data',...
'FontUnits','points',...
'BackgroundColor','none',...
'Color',[0 0 0],...
'DisplayName',blanks(0),...
'EdgeColor','none',...
'DVIMode','auto',...
'FontAngle','normal',...
'FontName','Helvetica',...
'FontSize',10,...
'FontWeight','normal',...
'HorizontalAlignment','center',...
'LineStyle','-',...
'LineWidth',0.5,...
'Margin',2,...
'Position',[0.485480943738658 -4.38475499092559 1.00024509545847],...
'Rotation',0,...
'String',blanks(0),...
'Interpreter','tex',...
'VerticalAlignment','cap',...
'ButtonDownFcn',[],...
'CreateFcn', {@local_CreateFcn, [], appdata} ,...
'DeleteFcn',[],...
'BusyAction','queue',...
'HandleVisibility','off',...
'HelpTopicKey',blanks(0),...
'HitTest','on',...
'Interruptible','on',...
'SelectionHighlight','on',...
'Serializable','on',...
'Tag',blanks(0),...
'UserData',[],...
'Visible','on',...
'XLimInclude','on',...
'YLimInclude','on',...
'ZLimInclude','on',...
'CLimInclude','on',...
'ALimInclude','on',...
'IncludeRenderer','on',...
'Clipping','off');

appdata = [];
appdata.SerializedAnnotationV7 = struct(...
    'LegendInformation', struct(...
    'IconDisplayStyle', 'on'));

h5 = get(h2,'ylabel');

set(h5,...
'Parent',h2,...
'Units','data',...
'FontUnits','points',...
'BackgroundColor','none',...
'Color',[0 0 0],...
'DisplayName',blanks(0),...
'EdgeColor','none',...
'DVIMode','auto',...
'FontAngle','normal',...
'FontName','Helvetica',...
'FontSize',10,...
'FontWeight','normal',...
'HorizontalAlignment','center',...
'LineStyle','-',...
'LineWidth',0.5,...
'Margin',2,...
'Position',[-0.37114337568058 -0.0145190562613435 1.00024509545847],...
'Rotation',90,...
'String',blanks(0),...
'Interpreter','tex',...
'VerticalAlignment','bottom',...
'ButtonDownFcn',[],...
'CreateFcn', {@local_CreateFcn, [], appdata} ,...
'DeleteFcn',[],...
'BusyAction','queue',...
'HandleVisibility','off',...
'HelpTopicKey',blanks(0),...
'HitTest','on',...
'Interruptible','on',...
'SelectionHighlight','on',...
'Serializable','on',...
'Tag',blanks(0),...
'UserData',[],...
'Visible','on',...
'XLimInclude','on',...
'YLimInclude','on',...
'ZLimInclude','on',...
'CLimInclude','on',...
'ALimInclude','on',...
'IncludeRenderer','on',...
'Clipping','off');

appdata = [];
appdata.SerializedAnnotationV7 = struct(...
    'LegendInformation', struct(...
    'IconDisplayStyle', 'on'));

h6 = get(h2,'zlabel');

set(h6,...
'Parent',h2,...
'Units','data',...
'FontUnits','points',...
'BackgroundColor','none',...
'Color',[0 0 0],...
'DisplayName',blanks(0),...
'EdgeColor','none',...
'DVIMode','auto',...
'FontAngle','normal',...
'FontName','Helvetica',...
'FontSize',10,...
'FontWeight','normal',...
'HorizontalAlignment','right',...
'LineStyle','-',...
'LineWidth',0.5,...
'Margin',2,...
'Position',[-8.57441016333938 4.71869328493648 1.00024509545847],...
'Rotation',0,...
'String',blanks(0),...
'Interpreter','tex',...
'VerticalAlignment','middle',...
'ButtonDownFcn',[],...
'CreateFcn', {@local_CreateFcn, [], appdata} ,...
'DeleteFcn',[],...
'BusyAction','queue',...
'HandleVisibility','off',...
'HelpTopicKey',blanks(0),...
'HitTest','on',...
'Interruptible','on',...
'SelectionHighlight','on',...
'Serializable','on',...
'Tag',blanks(0),...
'UserData',[],...
'Visible','off',...
'XLimInclude','on',...
'YLimInclude','on',...
'ZLimInclude','on',...
'CLimInclude','on',...
'ALimInclude','on',...
'IncludeRenderer','on',...
'Clipping','off');

appdata = [];
appdata.lastValidTag = 'uipanel1';

h7 = uipanel(...
'Parent',h1,...
'Units','pixels',...
'Title','Draw one path',...
'Tag','uipanel1',...
'Clipping','on',...
'Position',[50 398 251 238],...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'edit2';

h8 = uicontrol(...
'Parent',h7,...
'Callback',@(hObject,eventdata)chaotic_scattering_gui('edit2_Callback',hObject,eventdata,guidata(hObject)),...
'Position',[61 65 129 36],...
'String','0.0',...
'Style','edit',...
'CreateFcn', {@local_CreateFcn, @(hObject,eventdata)chaotic_scattering_gui('edit2_CreateFcn',hObject,eventdata,guidata(hObject)), appdata} ,...
'Tag','edit2');

appdata = [];
appdata.lastValidTag = 'text1';

h9 = uicontrol(...
'Parent',h7,...
'CData',[],...
'Position',[42 178 158 23],...
'String','input ordinate',...
'Style','text',...
'UserData',[],...
'Tag','text1',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'edit1';

h10 = uicontrol(...
'Parent',h7,...
'Callback',@(hObject,eventdata)chaotic_scattering_gui('edit1_Callback',hObject,eventdata,guidata(hObject)),...
'Position',[62 140 128 36],...
'String','0.0',...
'Style','edit',...
'CreateFcn', {@local_CreateFcn, @(hObject,eventdata)chaotic_scattering_gui('edit1_CreateFcn',hObject,eventdata,guidata(hObject)), appdata} ,...
'Tag','edit1');

appdata = [];
appdata.lastValidTag = 'text2';

h11 = uicontrol(...
'Parent',h7,...
'Position',[29 101 203 23],...
'String','output angle (radian)',...
'Style','text',...
'Tag','text2',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'togglebutton1';

h12 = uicontrol(...
'Parent',h7,...
'Callback',@(hObject,eventdata)chaotic_scattering_gui('togglebutton1_Callback',hObject,eventdata,guidata(hObject)),...
'Position',[62 16 129 35],...
'String','start',...
'Style','togglebutton',...
'Tag','togglebutton1',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'uipanel2';

h13 = uipanel(...
'Parent',h1,...
'Units','pixels',...
'Title','Generate chaotic graph',...
'Tag','uipanel2',...
'Clipping','on',...
'Position',[50 126 251 250],...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'edit3';

h14 = uicontrol(...
'Parent',h13,...
'Callback',@(hObject,eventdata)chaotic_scattering_gui('edit3_Callback',hObject,eventdata,guidata(hObject)),...
'Position',[59 144 131 37],...
'String','0.0',...
'Style','edit',...
'CreateFcn', {@local_CreateFcn, @(hObject,eventdata)chaotic_scattering_gui('edit3_CreateFcn',hObject,eventdata,guidata(hObject)), appdata} ,...
'Tag','edit3');

appdata = [];
appdata.lastValidTag = 'text3';

h15 = uicontrol(...
'Parent',h13,...
'Position',[49 188 151 19],...
'String','start value',...
'Style','text',...
'Tag','text3',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'text4';

h16 = uicontrol(...
'Parent',h13,...
'Position',[42 112 158 22],...
'String','end value',...
'Style','text',...
'Tag','text4',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'edit4';

h17 = uicontrol(...
'Parent',h13,...
'Callback',@(hObject,eventdata)chaotic_scattering_gui('edit4_Callback',hObject,eventdata,guidata(hObject)),...
'Position',[60 72 131 37],...
'String','0.0',...
'Style','edit',...
'CreateFcn', {@local_CreateFcn, @(hObject,eventdata)chaotic_scattering_gui('edit4_CreateFcn',hObject,eventdata,guidata(hObject)), appdata} ,...
'Tag','edit4');

appdata = [];
appdata.lastValidTag = 'togglebutton2';

h18 = uicontrol(...
'Parent',h13,...
'Callback',@(hObject,eventdata)chaotic_scattering_gui('togglebutton2_Callback',hObject,eventdata,guidata(hObject)),...
'Position',[61 22 130 36],...
'String','generate',...
'Style','togglebutton',...
'Tag','togglebutton2',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'uipanel3';

h19 = uipanel(...
'Parent',h1,...
'Units','pixels',...
'Title','Compute fractal dimension',...
'Tag','uipanel3',...
'Clipping','on',...
'Position',[50 25 251 86],...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

appdata = [];
appdata.lastValidTag = 'togglebutton3';

h20 = uicontrol(...
'Parent',h19,...
'Callback',@(hObject,eventdata)chaotic_scattering_gui('togglebutton3_Callback',hObject,eventdata,guidata(hObject)),...
'Position',[60 15 129 36],...
'String','compute',...
'Style','togglebutton',...
'Tag','togglebutton3',...
'CreateFcn', {@local_CreateFcn, blanks(0), appdata} );

hsingleton = h1;

% --- Set application data first then calling the CreateFcn. 
function local_CreateFcn(hObject, eventdata, createfcn, appdata)

if ~isempty(appdata)
   names = fieldnames(appdata);
   for i=1:length(names)
       name = char(names(i));
       setappdata(hObject, name, getfield(appdata,name));
   end
end

if ~isempty(createfcn)
   if isa(createfcn,'function_handle')
       createfcn(hObject, eventdata);
   else
       eval(createfcn);
   end
end


% --- Handles default GUIDE GUI creation and callback dispatch
function varargout = gui_mainfcn(gui_State, varargin)

gui_StateFields =  {'gui_Name'
    'gui_Singleton'
    'gui_OpeningFcn'
    'gui_OutputFcn'
    'gui_LayoutFcn'
    'gui_Callback'};
gui_Mfile = '';
for i=1:length(gui_StateFields)
    if ~isfield(gui_State, gui_StateFields{i})
        error(message('MATLAB:guide:StateFieldNotFound', gui_StateFields{ i }, gui_Mfile));
    elseif isequal(gui_StateFields{i}, 'gui_Name')
        gui_Mfile = [gui_State.(gui_StateFields{i}), '.m'];
    end
end

numargin = length(varargin);

if numargin == 0
    % CHAOTIC_SCATTERING_GUI
    % create the GUI only if we are not in the process of loading it
    % already
    gui_Create = true;
elseif local_isInvokeActiveXCallback(gui_State, varargin{:})
    % CHAOTIC_SCATTERING_GUI(ACTIVEX,...)
    vin{1} = gui_State.gui_Name;
    vin{2} = [get(varargin{1}.Peer, 'Tag'), '_', varargin{end}];
    vin{3} = varargin{1};
    vin{4} = varargin{end-1};
    vin{5} = guidata(varargin{1}.Peer);
    feval(vin{:});
    return;
elseif local_isInvokeHGCallback(gui_State, varargin{:})
    % CHAOTIC_SCATTERING_GUI('CALLBACK',hObject,eventData,handles,...)
    gui_Create = false;
else
    % CHAOTIC_SCATTERING_GUI(...)
    % create the GUI and hand varargin to the openingfcn
    gui_Create = true;
end

if ~gui_Create
    % In design time, we need to mark all components possibly created in
    % the coming callback evaluation as non-serializable. This way, they
    % will not be brought into GUIDE and not be saved in the figure file
    % when running/saving the GUI from GUIDE.
    designEval = false;
    if (numargin>1 && ishghandle(varargin{2}))
        fig = varargin{2};
        while ~isempty(fig) && ~ishghandle(fig,'figure')
            fig = get(fig,'parent');
        end
        
        designEval = isappdata(0,'CreatingGUIDEFigure') || isprop(fig,'__GUIDEFigure');
    end
        
    if designEval
        beforeChildren = findall(fig);
    end
    
    % evaluate the callback now
    varargin{1} = gui_State.gui_Callback;
    if nargout
        [varargout{1:nargout}] = feval(varargin{:});
    else       
        feval(varargin{:});
    end
    
    % Set serializable of objects created in the above callback to off in
    % design time. Need to check whether figure handle is still valid in
    % case the figure is deleted during the callback dispatching.
    if designEval && ishghandle(fig)
        set(setdiff(findall(fig),beforeChildren), 'Serializable','off');
    end
else
    if gui_State.gui_Singleton
        gui_SingletonOpt = 'reuse';
    else
        gui_SingletonOpt = 'new';
    end

    % Check user passing 'visible' P/V pair first so that its value can be
    % used by oepnfig to prevent flickering
    gui_Visible = 'auto';
    gui_VisibleInput = '';
    for index=1:2:length(varargin)
        if length(varargin) == index || ~ischar(varargin{index})
            break;
        end

        % Recognize 'visible' P/V pair
        len1 = min(length('visible'),length(varargin{index}));
        len2 = min(length('off'),length(varargin{index+1}));
        if ischar(varargin{index+1}) && strncmpi(varargin{index},'visible',len1) && len2 > 1
            if strncmpi(varargin{index+1},'off',len2)
                gui_Visible = 'invisible';
                gui_VisibleInput = 'off';
            elseif strncmpi(varargin{index+1},'on',len2)
                gui_Visible = 'visible';
                gui_VisibleInput = 'on';
            end
        end
    end
    
    % Open fig file with stored settings.  Note: This executes all component
    % specific CreateFunctions with an empty HANDLES structure.

    
    % Do feval on layout code in m-file if it exists
    gui_Exported = ~isempty(gui_State.gui_LayoutFcn);
    % this application data is used to indicate the running mode of a GUIDE
    % GUI to distinguish it from the design mode of the GUI in GUIDE. it is
    % only used by actxproxy at this time.   
    setappdata(0,genvarname(['OpenGuiWhenRunning_', gui_State.gui_Name]),1);
    if gui_Exported
        gui_hFigure = feval(gui_State.gui_LayoutFcn, gui_SingletonOpt);

        % make figure invisible here so that the visibility of figure is
        % consistent in OpeningFcn in the exported GUI case
        if isempty(gui_VisibleInput)
            gui_VisibleInput = get(gui_hFigure,'Visible');
        end
        set(gui_hFigure,'Visible','off')

        % openfig (called by local_openfig below) does this for guis without
        % the LayoutFcn. Be sure to do it here so guis show up on screen.
        movegui(gui_hFigure,'onscreen');
    else
        gui_hFigure = local_openfig(gui_State.gui_Name, gui_SingletonOpt, gui_Visible);
        % If the figure has InGUIInitialization it was not completely created
        % on the last pass.  Delete this handle and try again.
        if isappdata(gui_hFigure, 'InGUIInitialization')
            delete(gui_hFigure);
            gui_hFigure = local_openfig(gui_State.gui_Name, gui_SingletonOpt, gui_Visible);
        end
    end
    if isappdata(0, genvarname(['OpenGuiWhenRunning_', gui_State.gui_Name]))
        rmappdata(0,genvarname(['OpenGuiWhenRunning_', gui_State.gui_Name]));
    end

    % Set flag to indicate starting GUI initialization
    setappdata(gui_hFigure,'InGUIInitialization',1);

    % Fetch GUIDE Application options
    gui_Options = getappdata(gui_hFigure,'GUIDEOptions');
    % Singleton setting in the GUI M-file takes priority if different
    gui_Options.singleton = gui_State.gui_Singleton;

    if ~isappdata(gui_hFigure,'GUIOnScreen')
        % Adjust background color
        if gui_Options.syscolorfig
            set(gui_hFigure,'Color', get(0,'DefaultUicontrolBackgroundColor'));
        end

        % Generate HANDLES structure and store with GUIDATA. If there is
        % user set GUI data already, keep that also.
        data = guidata(gui_hFigure);
        handles = guihandles(gui_hFigure);
        if ~isempty(handles)
            if isempty(data)
                data = handles;
            else
                names = fieldnames(handles);
                for k=1:length(names)
                    data.(char(names(k)))=handles.(char(names(k)));
                end
            end
        end
        guidata(gui_hFigure, data);
    end

    % Apply input P/V pairs other than 'visible'
    for index=1:2:length(varargin)
        if length(varargin) == index || ~ischar(varargin{index})
            break;
        end

        len1 = min(length('visible'),length(varargin{index}));
        if ~strncmpi(varargin{index},'visible',len1)
            try set(gui_hFigure, varargin{index}, varargin{index+1}), catch break, end
        end
    end

    % If handle visibility is set to 'callback', turn it on until finished
    % with OpeningFcn
    gui_HandleVisibility = get(gui_hFigure,'HandleVisibility');
    if strcmp(gui_HandleVisibility, 'callback')
        set(gui_hFigure,'HandleVisibility', 'on');
    end

    feval(gui_State.gui_OpeningFcn, gui_hFigure, [], guidata(gui_hFigure), varargin{:});

    if isscalar(gui_hFigure) && ishghandle(gui_hFigure)
        % Handle the default callbacks of predefined toolbar tools in this
        % GUI, if any
        guidemfile('restoreToolbarToolPredefinedCallback',gui_hFigure); 
        
        % Update handle visibility
        set(gui_hFigure,'HandleVisibility', gui_HandleVisibility);

        % Call openfig again to pick up the saved visibility or apply the
        % one passed in from the P/V pairs
        if ~gui_Exported
            gui_hFigure = local_openfig(gui_State.gui_Name, 'reuse',gui_Visible);
        elseif ~isempty(gui_VisibleInput)
            set(gui_hFigure,'Visible',gui_VisibleInput);
        end
        if strcmpi(get(gui_hFigure, 'Visible'), 'on')
            figure(gui_hFigure);
            
            if gui_Options.singleton
                setappdata(gui_hFigure,'GUIOnScreen', 1);
            end
        end

        % Done with GUI initialization
        if isappdata(gui_hFigure,'InGUIInitialization')
            rmappdata(gui_hFigure,'InGUIInitialization');
        end

        % If handle visibility is set to 'callback', turn it on until
        % finished with OutputFcn
        gui_HandleVisibility = get(gui_hFigure,'HandleVisibility');
        if strcmp(gui_HandleVisibility, 'callback')
            set(gui_hFigure,'HandleVisibility', 'on');
        end
        gui_Handles = guidata(gui_hFigure);
    else
        gui_Handles = [];
    end

    if nargout
        [varargout{1:nargout}] = feval(gui_State.gui_OutputFcn, gui_hFigure, [], gui_Handles);
    else
        feval(gui_State.gui_OutputFcn, gui_hFigure, [], gui_Handles);
    end

    if isscalar(gui_hFigure) && ishghandle(gui_hFigure)
        set(gui_hFigure,'HandleVisibility', gui_HandleVisibility);
    end
end

function gui_hFigure = local_openfig(name, singleton, visible)

% openfig with three arguments was new from R13. Try to call that first, if
% failed, try the old openfig.
if nargin('openfig') == 2
    % OPENFIG did not accept 3rd input argument until R13,
    % toggle default figure visible to prevent the figure
    % from showing up too soon.
    gui_OldDefaultVisible = get(0,'defaultFigureVisible');
    set(0,'defaultFigureVisible','off');
    gui_hFigure = openfig(name, singleton);
    set(0,'defaultFigureVisible',gui_OldDefaultVisible);
else
    gui_hFigure = openfig(name, singleton, visible);  
    %workaround for CreateFcn not called to create ActiveX
    if feature('HGUsingMATLABClasses')
        peers=findobj(findall(allchild(gui_hFigure)),'type','uicontrol','style','text');    
        for i=1:length(peers)
            if isappdata(peers(i),'Control')
                actxproxy(peers(i));
            end            
        end
    end
end

function result = local_isInvokeActiveXCallback(gui_State, varargin)

try
    result = ispc && iscom(varargin{1}) ...
             && isequal(varargin{1},gcbo);
catch
    result = false;
end

function result = local_isInvokeHGCallback(gui_State, varargin)

try
    fhandle = functions(gui_State.gui_Callback);
    result = ~isempty(findstr(gui_State.gui_Name,fhandle.file)) || ...
             (ischar(varargin{1}) ...
             && isequal(ishghandle(varargin{2}), 1) ...
             && (~isempty(strfind(varargin{1},[get(varargin{2}, 'Tag'), '_'])) || ...
                ~isempty(strfind(varargin{1}, '_CreateFcn'))) );
catch
    result = false;
end

function [ h ] = circles(x,y,r,varargin)
% h = circles(x,y,r,varargin) plots circles of radius r at points x and y. 
% x, y, and r can be scalars or N-D arrays.  
% 
% Chad Greene, March 2014. Updated August 2014. 
% University of Texas Institute for Geophysics. 
% 
%% Syntax 
%  circles(x,y,r)
%  circles(...,'points',numberOfPoints)
%  circles(...,'rotation',degreesRotation)
%  circles(...,'ColorProperty',ColorValue)
%  circles(...,'LineProperty',LineValue)
%  h = circles(...)
% 
%% Description
% 
% circles(x,y,r) plots circle(s) of radius or radii r centered at points given by 
% x and y.  Inputs x, y, and r may be any combination of scalar,
% vector, or 2D matrix, but dimensions of all nonscalar inputs must agree. 
% 
% circles(...,'points',numberOfPoints) allows specification of how many points to use 
% for the outline of each circle. Default value is 1000, but this may be
% increased to increase plotting resolution.  Or you may specify a small
% number (e.g. 4 to plot a square, 5 to plot a pentagon, etc.). 
% 
% circles(...,'rotation',degreesRotation) rotates the shape by a given
% degreesRotation, which can be a scalar or a matrix. This is useless for
% circles, but may be desired for polygons with a discernible number of corner points. 
% 
% circles(...,'ColorProperty',ColorValue) allows declaration of
% 'facecolor' or 'facealpha'
% as name-value pairs. Try declaring any fill property as name-value pairs. 
%
% circles(...,'LineProperty',LineValue) allows declaration of 'edgecolor', 
% 'linewidth', etc.
%
% h = circles(...) returns the handle(s) h of the plotted object(s). 
% 
% 
%% EXAMPLES: 
%
% Example 1: 
% circles(5,10,3)
% 
% % Example 2: 
% x = 2:7;
% y = [5,15,12,25,3,18]; 
% r = [3 4 5 5 7 3]; 
% figure
% circles(x,y,r)
% 
% % Example 3: 
% figure
% circles(1:10,5,2)
% 
% % Example 4: 
% figure
% circles(5,15,1:5,'facecolor','none')
% 
% % Example 5: 
% figure 
% circles(5,10,3,'facecolor','green')
% 
% % Example 6: 
% figure
% h = circles(5,10,3,'edgecolor',[.5 .2 .9])
% 
% % Example 7: 
% lat = repmat((10:-1:1)',1,10); 
% lon = repmat(1:10,10,1); 
% r = .4; 
% figure
% h1 = circles(lon,lat,r,'linewidth',4,'edgecolor','m','facecolor',[.6 .4 .8]);
% hold on;
% h2 = circles(1:.5:10,((1:.5:10).^2)/10,.12,'edgecolor','k','facecolor','none');
% axis equal 
% 
% % Example 8: Circles have corners
% This script approximates circles with 1000 points. If all those points
% are too complex for your Pentium-II, you can reduce the number of points
% used to make each circle.  If 1000 points is not high enough resolution,
% you can increase the number of points.  Or if you'd like to draw
% triangles or squares, or pentagons, you can significantly reduce the
% number of points. Let's try drawing a stop sign: 
% 
% figure
% h = circles(1,1,10,'points',8,'color','red'); 
% axis equal
% % and we see that our stop sign needs to be rotated a little bit, so we'll
% % delete the one we drew and try again: 
% delete(h)
% h = circles(1,1,10,'points',8,'color','red','rot',45/2); 
% text(1,1,'STOP','fontname','helvetica CY',...
%     'horizontalalignment','center','fontsize',140,...
%     'color','w','fontweight','bold')
% 
% figure
% circles([1 3 5],2,1,'points',4,'rot',[0 45 35])
% 
%
% TIPS: 
% 1. Include the name-value pair 'facecolor','none' to draw outlines
% (non-filled) circles. 
% 
% 2. Follow the circles command with axis equal to fix distorted circles. 
%
% See also: fill, patch, and scatter. 

%% Check inputs: 

assert(isnumeric(x),'Input x must be numeric.') 
assert(isnumeric(y),'Input y must be numeric.') 
assert(isnumeric(r),'Input r must be numeric.') 

if ~isscalar(x) && ~isscalar(y)
    assert(numel(x)==numel(y),'If neither x nor y is a scalar, their dimensions must match.')
end
if ~isscalar(x) && ~isscalar(r)
    assert(numel(x)==numel(r),'If neither x nor r is a scalar, their dimensions must match.')
end
if ~isscalar(r) && ~isscalar(y)
    assert(numel(r)==numel(y),'If neither y nor r is a scalar, their dimensions must match.')
end

%% Parse inputs: 

% Define number of points per circle: 
tmp = strcmpi(varargin,'points')|strcmpi(varargin,'NOP')|strcmpi(varargin,'corners')|...
    strncmpi(varargin,'vert',4); 
if any(tmp)
    NOP = varargin{find(tmp)+1}; 
    tmp(find(tmp)+1)=1; 
    varargin = varargin(~tmp); 
else
    NOP = 1000; % 1000 points on periphery by default 
end

% Define rotation
tmp = strncmpi(varargin,'rot',3);
if any(tmp)
    rotation = varargin{find(tmp)+1}; 
    assert(isnumeric(rotation)==1,'Rotation must be numeric.')
    rotation = rotation*pi/180; % converts to radians
    tmp(find(tmp)+1)=1; 
    varargin = varargin(~tmp); 
else
    rotation = 0; % no rotation by default.
end

% Be forgiving if the user enters "color" instead of "facecolor"
tmp = strcmpi(varargin,'color');
if any(tmp)
    varargin{tmp} = 'facecolor'; 
end

%% Begin operations:

% Make inputs column vectors: 
x = x(:); 
y = y(:);
r = r(:); 
rotation = rotation(:); 

% Determine how many circles to plot: 
numcircles = max([length(x) length(y) length(r) length(rotation)]); 

% Create redundant arrays to make the plotting loop easy: 
if length(x)<numcircles
    x(1:numcircles) = x; 
end

if length(y)<numcircles
    y(1:numcircles) = y; 
end

if length(r)<numcircles
    r(1:numcircles) = r; 
end

if length(rotation)<numcircles
    rotation(1:numcircles) = rotation; 
end

% Define an independent variable for drawing circle(s):
t = 2*pi/NOP*(1:NOP); 

% Query original hold state:
holdState = ishold; 
hold on; 

% Preallocate object handle: 
h = NaN(size(x)); 

% Plot circles singly: 
for n = 1:numcircles
    h(n) = fill(x(n)+r(n).*cos(t+rotation(n)), y(n)+r(n).*sin(t+rotation(n)),'',varargin{:});
end

% Return to original hold state: 
if ~holdState
    hold off
end

% Delete object handles if not requested by user: 
if nargout==0 
    clear h 
end

%%%%%%%%%%% compute_path function for computing output angle %%%%%%%%%%%
function [theta, sum_path, numBounces] = compute_path(input, drawing)

%% 15-digit floating point
format long e;

%% Initial number of reflexions   %%
%% set to 0 for the current input %%
numBounces=0;

%% Initial slope %%
a=0;

%% Initial coordinates %%
current=input;
x_begin=-4;y_begin=input;
x_final=-4;y_final=input;

%% Current last coordinates %%
t=x_begin;
yt=y_begin;

%% Initial length of the particle path %%
sum_path=0;    

%% Cooordinates of 3 disks : %%
%% Disk 1 on top             %%
%% Disk 2 on right           %%
%% Disk 3 on bottom          %%
d=2.5;
xc(1)=-d*realsqrt(3)/6; yc(1)=d/2;
xc(2)=d*realsqrt(3)/3; yc(2)=0;           
xc(3)=-d*realsqrt(3)/6; yc(3)=-d/2;

%% Compute the 3 discriminants %%
%% for the first reflexion     %%               
for j=1:3
q=-2*xc(j)+2*a*(current-yc(j));
k=((current-yc(j))^2)-1+xc(j)^2 ;
m=1+a^2;
delta(j)=q^2-(4*k*m);
end

%% Get the first disk collided %%
if (delta(1)>0 && delta(2)>0)
 i=1;
else if (delta(3)>0 && delta(2)>0)
 i=3;       
else if (delta(2)>0 && delta(1)>=0)
 i=2;
else if (delta(2)>0 && delta(3)>=0)
 i=2;
else if (delta(2)>0)
 i=2;
else if (delta(1)>0)
 i=1;
else if (delta(3)>0)
 i=3;
end
end
end
end  
end
end
end

%% Secondary loop : until particle is not collided with a disk %%
while ( (delta(1)>0 && delta(2)<=0 && delta(3)<=0) || (delta(2)>0 && delta(1)<0 && delta(3)<0) || (delta(3)>0 && delta(2)<=0 && delta(1)<=0) || (delta(1)>=0 && delta(2)>0) || (delta(3)>=0 && delta(2)>0) || (delta(1)>0 && delta(3)>0)) 

%% Break si current = 0 %%
if (current == 0)

 %% Drawing path between (x_begin,y_begin) and (x_final,y_final)
 if (strcmp(drawing,'draw'))
   vector1=[-4,xc(2)-1];
   vector2=[0,0];
   line(vector1,vector2);
   pause(0.5);
   %% For last segment 
   t=-4;
   yt=0;
 end
 %% Length of path %%
 sum_path=2*(xc(2)-1)+4;                
 break;
end

%% Break if current = 1.25 || -1.25 %%
if (current == 1.25 || current == -1.25)

 %% Drawing path between (x_begin,y_begin) and (x_final,y_final)
 if (strcmp(drawing,'draw'))
   vector1=[-4,xc(1)-1];
   vector2=[current,current];
   line(vector1,vector2);
   pause(0.5);
   %% For last segment 
   t=-4;
   yt=current;
 end
 %% Length of path %%
 sum_path=2*(xc(1)-1)+4;                
 break;
end

 %% Compute of intersection point for the l-th reflexion %%
 q=-2*xc(i)+2*a*(current-yc(i));
 k=((current-yc(i))^2)-1+xc(i)^2 ;
 m=1+a^2;
 delta(i)=q^2-(4*k*m);

 % Secondary equation : "a" and "b" %%
 a_l=m; b_l=q;

 % 2 solutions of secondary equation %%
 x1_l=(-b_l+realsqrt(delta(i)))/(2*a_l);
 x2_l=(-b_l-realsqrt(delta(i)))/(2*a_l);

 y1_l=a*x1_l+current;
 y2_l=a*x2_l+current;

 %% Save x_begin and y_begin %% 
 %% from x_final and y_final %%       
 x_begin=x_final;
 y_begin=y_final;

 %% Take the solution of minimum distance %%
 if (numBounces>0)
 if ((x_begin-x1_l)^2+(y_begin-y1_l)^2)>=((x_begin-x2_l)^2+(y_begin-y2_l)^2)  
 x_final=x2_l;
 y_final=y2_l;
 else
 x_final=x1_l;
 y_final=y1_l;
 end
 else
 %% Case for the first input %%
 x_final=min([x1_l,x2_l]);
 y_final=current;
 end

 %% Drawing path between (x_begin,y_begin) and (x_final,y_final)
 if (strcmp(drawing,'draw'))
 vector1=[x_begin,x_final];
 vector2=[y_begin,y_final];
 line(vector1,vector2);
 pause(0.5);
 end

 %% Sum the length of path %%
 distance=realsqrt((x_final-x_begin)^2+(y_final-x_begin)^2);
 sum_path=sum_path+distance;                

 %% Increment the number of reflexion %%
 numBounces=numBounces+1;

 %% Compute the new straight line equation %%
 p=(y_final-yc(i))/(x_final-xc(i));  
 a=(2*p-a+a*p^2)/(1-p^2+2*p*a);    
 current=y_final-a*x_final;

 % Compute the 3 new discriminants %%
 for k=1:3
 f(k)=-2*xc(k)+2*a*(current-yc(k));
 g(k)=((current-yc(k))^2)-1+xc(k)^2 ;
 h(k)=1+a^2;   
 delta(k)=f(k)^2-(4*g(k)*h(k));
 end

 %% Searching for the indice of new disk collided %%
 %% Delta position %%
 pos_delta=0.00000001;
 %% Adaptative step deplacement of straight    %%
 %% line for a high value of slope :           %%
 %% minimum of 10000 points on smaller distance %%
 %% between disks (disk 1 and disk 3) = 0.5    %%
 step=0.5/(10000*abs(a));
 % Affine parameter %%
 t=x_final;
 %% Ordinate start %%
 yt=y_final;
 %% Boolean for intersection %%
 intersection='false';

 %% Special case where there is only one delta > 0 %%
 if ((i==1) && (delta(1)>0) && (delta(2)<0) && (delta(3)<0))
  x_shift=x_final+pos_delta;
  y_shift=a*x_shift+current;
 if (((x_shift-xc(i))^2+(y_shift-yc(i))^2)<1)
 delta_line=-step; 
 else
 delta_line=step;
 end

 while ((yt<4) && (yt>-4) && (t<4) && (t>-4))   
 t=t+delta_line;
 yt=a*t+current;
 end

 if strcmp(intersection,'false')
  if ((((yt>4) && (t<4)) || ((yt>0) && (t>4))) && (delta_line>0))
  direction='topright';
 break; 
 else if ((((yt>4) && (t>-4)) || ((yt>0) && (t<-4))) && (delta_line<0))
 direction='topleft';
 break;
 else if ((((yt<-4) && (t<4)) || ((yt<0) && (t>4)))  && (delta_line>0))
 direction='bottomright';
 break;
 else if ((((yt<-4) && (t>-4)) || ((yt<0) && (t<-4))) && (delta_line<0))
 direction='bottomleft';
 break;
 end
 end
 end
 end
 end  

 %% Special case where there is only one delta > 0 %%
 else if ((i==2) && (delta(2)>0) && (delta(1)<0) && (delta(3)<0))
  x_shift=x_final+pos_delta;
  y_shift=a*x_shift+current;
 if (((x_shift-xc(i))^2+(y_shift-yc(i))^2)<1)
 delta_line=-step; 
 else
 delta_line=step;
 end

 while ((yt<4) && (yt>-4) && (t<4) && (t>-4))   
 t=t+delta_line;
 yt=a*t+current;
 end

 if strcmp(intersection,'false')
  if ((((yt>4) && (t<4)) || ((yt>0) && (t>4))) && (delta_line>0))
  direction='topright';
 break; 
 else if ((((yt>4) && (t>-4)) || ((yt>0) && (t<-4))) && (delta_line<0))
 direction='topleft';
 break;
 else if ((((yt<-4) && (t<4)) || ((yt<0) && (t>4)))  && (delta_line>0))
 direction='bottomright';
 break;
 else if ((((yt<-4) && (t>-4)) || ((yt<0) && (t<-4))) && (delta_line<0))
 direction='bottomleft';
 break;
 end
 end
 end
 end
 end  

 %% Special case where there is only one delta > 0 %%
 else if ((i==3) && (delta(3)>0) && (delta(1)<0) && (delta(2)<0))
  x_shift=x_final+pos_delta;
  y_shift=a*x_shift+current;
  if (((x_shift-xc(i))^2+(y_shift-yc(i))^2)<1)
  delta_line=-step; 
  else
  delta_line=step;
  end

 while ((yt<4) && (yt>-4))   
 t=t+delta_line;
 yt=a*t+current;
 end

 if strcmp(intersection,'false')
  if ((((yt>4) && (t<4)) || ((yt>0) && (t>4))) && (delta_line>0))
  direction='topright';
 break; 
 else if ((((yt>4) && (t>-4)) || ((yt>0) && (t<-4))) && (delta_line<0))
 direction='topleft';
 break;
 else if ((((yt<-4) && (t<4)) || ((yt<0) && (t>4)))  && (delta_line>0))
 direction='bottomright';
 break;
 else if ((((yt<-4) && (t>-4)) || ((yt<0) && (t<-4))) && (delta_line<0))
 direction='bottomleft';
 break;
 end
 end
 end
 end
 end  

 %% General case with i=1 (top disk) %%
 else if (((i==1) && (delta(1)>=0) && (delta(2)>=0)) || ((i==1) && (delta(1)>=0) && (delta(3)>=0)))
  x_shift=x_final+pos_delta;
  y_shift=a*x_shift+current;
 if (((x_shift-xc(i))^2+(y_shift-yc(i))^2)<1)
 delta_line=-step; 
 else
 delta_line=step;
 end

 while ((yt<4) && (yt>-4) && (t<4) && (t>-4) && strcmp(intersection,'false'))   
 t=t+delta_line;
 yt=a*t+current;

 if (((t-xc(1))^2+(yt-yc(1))^2)<1)
  i=1; 
  intersection='true';
 else if (((t-xc(2))^2+(yt-yc(2))^2)<1)
 i=2;
 intersection='true';
 else if (((t-xc(3))^2+(yt-yc(3))^2)<1)
 i=3;
 intersection='true';
 end
 end
 end  
 end

 if strcmp(intersection,'false')
  if ((((yt>4) && (t<4)) || ((yt>0) && (t>4))) && (delta_line>0))
  direction='topright';
 break; 
 else if ((((yt>4) && (t>-4)) || ((yt>0) && (t<-4))) && (delta_line<0))
 direction='topleft';
 break;
 else if ((((yt<-4) && (t<4)) || ((yt<0) && (t>4)))  && (delta_line>0))
 direction='bottomright';
 break;
 else if ((((yt<-4) && (t>-4)) || ((yt<0) && (t<-4))) && (delta_line<0))
 direction='bottomleft';
 break;
 end
 end
 end
 end
 end  

 %% General case with i=2 (right disk) %%
 else if (((i==2) && (delta(2)>=0) && (delta(1)>=0)) || ((i==2) && (delta(2)>=0) && (delta(3)>=0)))
  x_shift=x_final+pos_delta;
  y_shift=a*x_shift+current;
 if (((x_shift-xc(i))^2+(y_shift-yc(i))^2)<1)
 delta_line=-step; 
 else
 delta_line=step;
 end

 while ((yt<4) && (yt>-4) && (t<4) && (t>-4) && strcmp(intersection,'false'))   
 t=t+delta_line;
 yt=a*t+current;

 if (((t-xc(1))^2+(yt-yc(1))^2)<1)
  i=1; 
  intersection='true';
 else if (((t-xc(2))^2+(yt-yc(2))^2)<1)
 i=2;
 intersection='true';
 else if (((t-xc(3))^2+(yt-yc(3))^2)<1)
 i=3;
 intersection='true';
 end
 end
 end  
 end

 if strcmp(intersection,'false')
  if ((((yt>4) && (t<4)) || ((yt>0) && (t>4))) && (delta_line>0))
  direction='topright';
 break; 
 else if ((((yt>4) && (t>-4)) || ((yt>0) && (t<-4))) && (delta_line<0))
 direction='topleft';
 break;
 else if ((((yt<-4) && (t<4)) || ((yt<0) && (t>4)))  && (delta_line>0))
 direction='bottomright';
 break;
 else if ((((yt<-4) && (t>-4)) || ((yt<0) && (t<-4))) && (delta_line<0))
 direction='bottomleft';
 break;
 end
 end
 end
 end
 end  

 %% General case with i=3 (bottom disk) %%
 else if (((i==3) && (delta(3)>=0) && (delta(1)>=0)) || ((i==3) && (delta(3)>=0) && (delta(2)>=0)))
  x_shift=x_final+pos_delta;
  y_shift=a*x_shift+current;
 if (((x_shift-xc(i))^2+(y_shift-yc(i))^2)<1)
 delta_line=-step; 
 else
 delta_line=step;
 end

 while ((yt<4) && (yt>-4) && (t<4) && (t>-4) && strcmp(intersection,'false'))   
 t=t+delta_line;
 yt=a*t+current;

 if (((t-xc(1))^2+(yt-yc(1))^2)<1)
  i=1; 
  intersection='true';
 else if (((t-xc(2))^2+(yt-yc(2))^2)<1)
 i=2;
 intersection='true';
 else if (((t-xc(3))^2+(yt-yc(3))^2)<1)
 i=3;
 intersection='true';
 end
 end
 end  
 end

 if strcmp(intersection,'false')
  if ((((yt>4) && (t<4)) || ((yt>0) && (t>4))) && (delta_line>0))
  direction='topright';
  break; 
 else if ((((yt>4) && (t>-4)) || ((yt>0) && (t<-4))) && (delta_line<0))
 direction='topleft';
 break;
 else if ((((yt<-4) && (t<4)) || ((yt<0) && (t>4)))  && (delta_line>0))
 direction='bottomright';
 break;
 else if ((((yt<-4) && (t>-4)) || ((yt<0) && (t<-4))) && (delta_line<0))
 direction='bottomleft';
 break;
 end
 end
 end
 end
 end  
 end
 end
 end
 end
 end
 end
 end

%% 4 possible directions for output %%
if (a==0)
 theta=pi;
 else if (strcmp(direction,'topright'))
 theta=atan(a);
 else if (strcmp(direction,'topleft'))
 theta=atan(a)+pi;
 else if (strcmp(direction,'bottomright'))
 theta=atan(a)+2*pi;
 else if (strcmp(direction,'bottomleft'))
 theta=atan(a)+pi;
 end
 end
 end
 end
 end

 %% Drawing last segment 
 if (strcmp(drawing,'draw'))
 vector1=[x_final,t];
 vector2=[y_final,yt];
 line(vector1,vector2);
 pause(0.5);
 end

function initGraphics()

% Plot axes
box on;
x(1:8/0.05+1)=0;
y=(-4:0.05:4);
plot(x,y,':k','LineWidth',1.5);
hold on;
x=(-4:0.05:4);
y(1:8/0.05+1)=0;
plot(x,y,':k','LineWidth',1.5);
xlim([-4 4]);
ylim([-4 4]);

% Center coordinates of 3 circles
d=2.5;
xc(1)=-d*realsqrt(3)/6;
yc(1)=d/2;
xc(2)=d*realsqrt(3)/3;
yc(2)=0;
xc(3)=-d*realsqrt(3)/6;
yc(3)=-d/2;

% Plot circles and fill them
radius=1;
grey = hex2rgb('#B3B3BC'); 
for i=1:3
circles(xc(i),yc(i),radius,'color',grey);
end

function [ rgb ] = hex2rgb(hex,range)
% hex2rgb converts hex color values to rgb arrays on the range 0 to 1. 
% 
% 
% * * * * * * * * * * * * * * * * * * * * 
% SYNTAX:
% rgb = hex2rgb(hex) returns rgb color values in an n x 3 array. Values are
%                    scaled from 0 to 1 by default. 
%                    
% rgb = hex2rgb(hex,256) returns RGB values scaled from 0 to 255. 
% 
% 
% * * * * * * * * * * * * * * * * * * * * 
% EXAMPLES: 
% 
% myrgbvalue = hex2rgb('#334D66')
%    = 0.2000    0.3020    0.4000
% 
% 
% myrgbvalue = hex2rgb('334D66')  % <-the # sign is optional 
%    = 0.2000    0.3020    0.4000
% 
%
% myRGBvalue = hex2rgb('#334D66',256)
%    = 51    77   102
% 
% 
% myhexvalues = ['#334D66';'#8099B3';'#CC9933';'#3333E6'];
% myrgbvalues = hex2rgb(myhexvalues)
%    =   0.2000    0.3020    0.4000
%        0.5020    0.6000    0.7020
%        0.8000    0.6000    0.2000
%        0.2000    0.2000    0.9020
% 
% 
% myhexvalues = ['#334D66';'#8099B3';'#CC9933';'#3333E6'];
% myRGBvalues = hex2rgb(myhexvalues,256)
%    =   51    77   102
%       128   153   179
%       204   153    51
%        51    51   230
% 
% HexValsAsACharacterArray = {'#334D66';'#8099B3';'#CC9933';'#3333E6'}; 
% rgbvals = hex2rgb(HexValsAsACharacterArray)
% 
% * * * * * * * * * * * * * * * * * * * * 
% Chad A. Greene, April 2014
%
% Updated August 2014: Functionality remains exactly the same, but it's a
% little more efficient and more robust. Thanks to Stephen Cobeldick for
% the improvement tips. In this update, the documentation now shows that
% the range may be set to 256. This is more intuitive than the previous
% style, which scaled values from 0 to 255 with range set to 255.  Now you
% can enter 256 or 255 for the range, and the answer will be the same--rgb
% values scaled from 0 to 255. Function now also accepts character arrays
% as input. 
% 
% * * * * * * * * * * * * * * * * * * * * 
% See also rgb2hex, dec2hex, hex2num, and ColorSpec. 
% 

%% Input checks:

assert(nargin>0&nargin<3,'hex2rgb function must have one or two inputs.') 

if nargin==2
assert(isscalar(range)==1,'Range must be a scalar, either "1" to scale from 0 to 1 or "256" to scale from 0 to 255.')
end

%% Tweak inputs if necessary: 

if iscell(hex)
assert(isvector(hex)==1,'Unexpected dimensions of input hex values.')

% In case cell array elements are separated by a comma instead of a
% semicolon, reshape hex:
if isrow(hex)
hex = hex'; 
end

% If input is cell, convert to matrix: 
hex = cell2mat(hex);
end

if strcmpi(hex(1,1),'#')
hex(:,1) = [];
end

if nargin == 1
range = 1; 
end

%% Convert from hex to rgb: 

switch range
case 1
rgb = reshape(sscanf(hex.','%2x'),3,[]).'/255;

case {255,256}
rgb = reshape(sscanf(hex.','%2x'),3,[]).';

otherwise
error('Range must be either "1" to scale from 0 to 1 or "256" to scale from 0 to 255.')
end