function varargout = ising(varargin)
% ISING M-file for ising.fig
% 2D Ising Model
% Copyright 2002-2003 The MathWorks, Inc.
% Last Modified by GUIDE v2.5 21-Aug-2006 23:18:49

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ising_OpeningFcn, ...
                   'gui_OutputFcn',  @ising_OutputFcn, ...
                   '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 ising is made visible.
function ising_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   command line arguments to ising (see VARARGIN)

global temp_cur temp_new temperature stop reset scaling;

temp_cur = 2.269;
temp_new = 2.269;
temperature = temp_cur;
stop = 1;
reset = 1;
scaling = 10;

set(handles.text_current, 'String', num2str(temp_cur));
set(handles.text_newtemp, 'String', num2str(temp_cur));
set(handles.popupmenu,'Value',3)

% display image.  
handles.axes_ising = image([1;1],'Tag','image');
%The reason we multiply by 5 and add 2 is just to pick colors.  It has no meaning.  The two different colors represent the two different spins.
axis off  % Since the axis do not carry meaning, we will turn them off. 
set(handles.text_step,'String', '0');

disp('Welcome to 2D Ising Simulation');

% --- Outputs from this function are returned to the command line.
function varargout = ising_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)



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

% Update handles structure
guidata(hObject, handles);

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

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

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

global temperature stop reset scaling;

stop = 5;
reset = 5;

disp (strcat('RUN: Running Ising Model Simulation with temperature=',num2str(temperature))) % shows up on your command window
n=100; % number of rows in the lattice
o=100; % number of columns in the lattice

%Generate a random initial configuration

    for i=1:n  %for every row
        for j=1:o  % for every column
            r1(i)=rand*n;  % pick a random row coordinate
            r12(j)= rand*o;  % pick a random column coordinate

            if r1(i)<(n/2)   % Assign a value of 1 or -1 to each coordinate. The probability of either choice is 1/2.   
                if r12(j)<(o/2)
                    A(i,j)=-1;  % The lattice of points is called A.
                end
            else
                A(i,j)=1;
            end
            if r1(i)>(n/2)
                if r12(j)>(o/2)
                    A(i,j)=1;
                end
            else
                A(i,j)=-1;
            end
        end
    end
  % end of generating a random initial configuration
    
    % Assign A1 to be the same as A. Now we can make a change in A.  
    % If we don't decide to keep the change, we can revert to A1 
    % by assigning A to be same as A1.
    A1=A;  
   count = 0;
    while ((stop > 0) && (reset > 0))  % 10,000 images will be shown
        m = 0;
        count = count+1;
        set(handles.text_step,'String', num2str(count));
        while ((stop > 0) && (reset > 0) && (m < 1000))% 1000 steps between each image 
            m = m+1;

            % Pick a random point to flip


            r2=ceil(rand*(n)); % pick a random horizontal position
            r21= ceil(rand*(o)); % pick a random vertical position



            % Calculate energy around A(r2, r21)



            if r2 ~= 1 & r2 ~= n & r21~=1 & r21~=o  %  exclude boundary points for now

                Elatticepoint1 = -A(r2,r21)*( A(r2+1, r21)+A(r2-1, r21)+A(r2, r21+1)+A(r2, r21-1)  ); % energy of the current lattice 

            else  % each boundary point case is considered
                if r2 == 1 
                    AQ = -A(r2,r21)*A(n, r21);
                else
                    AQ = -A(r2,r21)*A(r2-1, r21);
                end


                if r2 == n
                    AW = -A(r2,r21)*A(1, r21);
                else
                    AW = -A(r2,r21)*A(r2+1, r21);
                end


                if r21==1
                    AE = -A(r2,r21)*A(r2, o);
                else
                    AE = -A(r2,r21)*A(r2, r21-1);
                end


                if  r21==o
                    AR = -A(r2,r21)*A(r2, 1);
                else
                    AR = -A(r2,r21)*A(r2, r21+1);
                end

                Elatticepoint1= AQ+ AW + AE+ AR; % energy of the current lattice 


            end

            % Cause change
            if A(r2,r21) ==1
                A(r2,r21)=-1;
            else A(r2,r21)=1;
            end
            A2=A;  % Assign A2 to be the perturbed configuration.  If we decide to keep the change, we will let A1 be the same as A2.

            % Calculate energy around A(r2, r21).  This is done exactly the
            % same way as the previous energy calculation.


            if r2 ~= 1 & r2 ~= n & r21~=1 & r21~=o

                Elatticepoint2 = -A(r2,r21)*( A(r2+1, r21)+A(r2-1, r21)+A(r2, r21+1)+A(r2, r21-1)  );% energy of the lattice with a suggested change

            else
                if r2 == 1
                    AQ = -A(r2,r21)*A(n, r21);
                else
                    AQ = -A(r2,r21)*A(r2-1, r21);
                end


                if r2 == n
                    AW = -A(r2,r21)*A(1, r21);
                else
                    AW = -A(r2,r21)*A(r2+1, r21);
                end


                if r21==1
                    AE = -A(r2,r21)*A(r2, o);
                else
                    AE = -A(r2,r21)*A(r2, r21-1);
                end


                if  r21==o
                    AR = -A(r2,r21)*A(r2, 1);
                else
                    AR = -A(r2,r21)*A(r2, r21+1);
                end

                Elatticepoint2= AQ+ AW + AE+ AR; % energy of the lattice with a suggested change


            end

            % Decide if new configuration should be kept.  New configuration is
            % kept if rand < P.  P is the ratio of probability of occurence of perturbed
            % setup over occurence of current setup.  Probability of a setup is
            % calculated by Gibbs expression.

            r4=rand;
            changeE = Elatticepoint2-Elatticepoint1; % calculate difference in energy of the lattice due to the suggested change
            P = exp(-changeE/temperature);  % calculate probability of accepting a higher energy change

            if changeE < 0
                A1=A2;  % if the suggested change leads to lower energy of system, we will definitely keep the change
            end

            if r4 < P
                A1=A2;  % if the suggested change leads to a higher energy system, we will keep the change if the random number is less than the probability we just calculated.
            end
            A=A1;
        end
        h1 = image(5*(A1+2),'Tag','image'); % display image.  
        handles.axes_ising = h1;
        %The reason we multiply by 5 and add 2 is just to pick colors.  It has no meaning.  The two different colors represent the two different spins.
        axis off  % Since the axis do not carry meaning, we will turn them off.  
        pause(0.01*scaling) % pause 0.01 seconds before displaying each image
    end


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

global stop;
stop = -1;

disp('STOP: Simulation stopped, click Run to start over');

% --- Executes on slider movement.
function slider_temp_Callback(hObject, eventdata, handles)
% hObject    handle to slider_temp (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,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

global temp_new;

temp_new = get(hObject, 'Value');
set(handles.text_newtemp, 'String', num2str(temp_new));

set(handles.text_newtemplbl, 'ForegroundColor', 'red')


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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

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

global temp_new temp_cur reset temperature;

if temp_cur == temp_new reset = 1;
    disp('RESET: New temperature is the same as current, nothing to reset')
else

    disp(strcat('RESET: Temperature changed from t', num2str(temp_cur),...
            ' to t=', num2str(temp_new)));
    temp_cur = temp_new;

    set(handles.text_current, 'String', num2str(temp_cur));
    set(handles.text_newtemplbl, 'ForegroundColor', 'blue');

    % This parameter can be changed to observe paramagnetism and ferromagnetism.  
    % The threshold value is 2.269.
    temperature = temp_cur;
    reset = -1;

end

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

global stop;

stop = -1;
pause(0.2);
close all;
disp('EXIT: Exited Ising')

% --- Executes on selection change in popupmenu.
function popupmenu_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popupmenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu

global scaling;

value = get(hObject,'Value');

if value == 1 
    scaling = 25;
else if value == 2
        scaling = 20;
    else if value == 3
            scaling = 10;
        else if value == 4
                scaling = 5;
            else if value == 5
                    scaling = 1;
                end
            end
        end
    end
end
                   
% --- Executes during object creation, after setting all properties.
function popupmenu_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu 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