classdef HW13_Solution_2DIsing < handle
    %HW13_Solution_2DIsing 2D Ising model simulation
    %   Configurations are updated using MCMC
    
    properties
        B = 0;              %Magnetic Field
        J = -.1;              %Spin Interaction Energy
        Beta = 1;           %1/KbT
        UpdateInterval = 100;%How often to show updated spin state
        NumJumps = 10000;  %Number of MCMC jumps
        NSize= 100;         %Number of spins in each row and column
        Spins;              %This is Spin State Matrix
        
        FigH = 0;
        ShowCorFit = 1;     %Show fitting of Autocorrelation to exp model
        
        GuiGo=1;
    end
    
    methods
        function runChain(obj,InGUI)
            
            if nargin <2
                InGUI=0;
            end
            
            %Start with random orientation
            obj.Spins=double(randn(obj.NSize,obj.NSize)>0);
            obj.Spins(obj.Spins==0)=-1;
            
            if InGUI
                obj.showSpinsGUI();
            else
                obj.showSpins();
            end
            
            for nn=1:obj.NumJumps
                
                %Stop if GUI stop button is pressed
                if ~obj.GuiGo; break; end
                
                %Make test jump
                
                %Pick site and flip spin
                IDX=ceil(obj.NSize*rand);
                IDY=ceil(obj.NSize*rand);
                SpinsTest=obj.Spins;
                SpinsTest(IDX,IDY)=-obj.Spins(IDX,IDY);
                
                %Energy from nearest neighbor interactions
                E_Spin=obj.J*(sum(sum(obj.Spins(1:end-1,:).*obj.Spins(2:end,:)))+... %Y neighbor
                    sum(obj.Spins(1,:).*obj.Spins(end,:))+... %Y periodic BC
                    sum(sum(obj.Spins(:,1:end-1).*obj.Spins(:,2:end)))+... %X neighbor
                    sum(obj.Spins(:,1).*obj.Spins(:,end))); %X periodic BC
                
                E_SpinTest=obj.J*(sum(sum(SpinsTest(1:end-1,:).*SpinsTest(2:end,:)))+... %Y neighbor
                    sum(SpinsTest(1,:).*SpinsTest(end,:))+... %Y periodic BC
                    sum(sum(SpinsTest(:,1:end-1).*SpinsTest(:,2:end)))+... %X neighbor
                    sum(SpinsTest(:,1).*SpinsTest(:,end))); %X periodic BC
                
                %Energy from Interactions with feild
                E_Field=-obj.B*sum(obj.Spins(:));
                E_FieldTest=-obj.B*sum(SpinsTest(:));
                
                %Total Energy (Hamiltonian)
                H=E_Spin+E_Field;
                HTest=E_SpinTest+E_FieldTest;
                
                %Test acceptance probability for Metropolis-Hastings
                Alpha=exp(-obj.Beta*(HTest-H));
                
                %Metropolis-Hastings rule
                if rand()< (Alpha)
                    obj.Spins(IDX,IDY)=-obj.Spins(IDX,IDY);
                end
                
                %Show every ShowIter interations
                if nn/obj.UpdateInterval==round(nn/obj.UpdateInterval)
                    if InGUI
                        obj.showSpinsGUI();
                    else
                    FE=obj.showSpins(1);
                    if ~FE;break;end
                    end
                end
                
            end
        end
        
        function L=findCorLength(obj)
            
            AutoCor=real(ifft2(fft2(obj.Spins).*conj(fft2(obj.Spins))));
            %Shift zero to center
            ACS=fftshift(AutoCor);
            
            %Make X and Y array
            N=size(AutoCor,1);
            [X Y]=meshgrid((-N/2:N/2-1),(-N/2:N/2-1));
            %Make R array
            R=sqrt(X.^2+Y.^2);
            
            %Do radial averaging of Auto Correlation
            RbinSize=1;
            
            for rr=0:obj.NSize/2
                Rmask=(R<(rr+RbinSize))&(R>=rr);
                Rave(rr+1)=mean(ACS(Rmask));
            end
            
            Rvec=(0:obj.NSize/2);
            
            %Fit curve to exponential model
            X0=[Rave(1)  Rave(find(Rave<Rave(1)/exp(1)))];
            Optim=optimset();
            Xfound=fminsearch(@obj.fitmodel,X0,Optim,Rave,Rvec);
            L=abs(Xfound(2));
            
            if obj.ShowCorFit
                figure
                plot(Rvec,Rave)
                xlabel('distance (Pixels)')
                ylabel('correlation')
                
                hold on
                [~,Model]=obj.fitmodel(Xfound,Rave,Rvec);
                plot(Rvec,Model)
                legend('AutoCorrelation','Fit')
                
            end
            
        end
        
        function [MSE,Model]=fitmodel(obj,X,Rave,Rvec)
            Model=X(1)*exp(-Rvec/abs(X(2)));
            Data = Rave;
            MSE=mean( (Model-Data).^2);
        end
        
        function [FigExists]=showSpins(obj,InLoop)
            
            if nargin <2
                InLoop=0;
            end
            
            
            
            if InLoop
                if (ishandle(obj.FigH) && strcmp(get(obj.FigH, 'type'), 'figure'))
                    figure(obj.FigH)
                    imshow(obj.Spins,[],'InitialMagnification',400)
                    drawnow();
                    FigExists=1;
                else
                    FigExists=0;
                end
            else
                if ~(ishandle(obj.FigH) && strcmp(get(obj.FigH, 'type'), 'figure'))
                    obj.FigH=figure;
                end
                figure(obj.FigH)
                imshow(obj.Spins,[],'InitialMagnification',400)
                drawnow();
                FigExists=1;
            end
        end
        
        function showSpinsGUI(obj)
            imshow(obj.Spins,[],'InitialMagnification',400)
            drawnow();
        end
        
        
        function gui(obj)
            
            %Create figure and position on the screen
            %Position is [left bottom width height]
            GUIFigH=figure('Position',[100 100 600 400])
            GUIFigH.CloseRequestFcn=@myclose;
            
            %Create axis for our plot
            AxH=axes('Units','Pixels')
            %Set axis position
            AxH.Position = [150 50 400 300]
            
            %Setup our edit boxes
            EditTxtJ = uicontrol('Style','edit',...
                'Position',[70 350 50 20],...
                'String','J');
            uicontrol('Style','text',...
                'Position',[10 350 60 20],...
                'String','J');
            
            
            EditTxtB = uicontrol('Style','edit',...
                'Position',[70 300 50 20],...
                'String','B');
            uicontrol('Style','text',...
                'Position',[10 300 60 20],...
                'String','B');
            
            EditTxtNumJumps = uicontrol('Style','edit',...
                'Position',[70 250 50 20],...
                'String','NumJumps');
            uicontrol('Style','text',...
                'Position',[10 250 60 20],...
                'String','NumJumps');
            
            
            %Setup pushbutton
            PushButton = uicontrol('Style','pushbutton',...
                'Position',[20 200 75 20],...
                'String','Start/Reset');
            PushButton.Callback=@start;
            
            %Setup pushbutton
            PushButton = uicontrol('Style','pushbutton',...
                'Position',[20 150 75 20],...
                'String','Stop');
            PushButton.Callback=@stop;
            
            
            
            %Put object properties (A and B) into edit boxes
            params2gui()
            
            
            %This function is called when figure is closed
            function myclose(~,~) %(~,~) is required to catch two inputs
                gui2params();
                stop();  
                delete(GUIFigH)
            end
            
            %Update object properties with GUI values.
            function gui2params()
                obj.B=str2double(EditTxtB.String);
                obj.J=str2double(EditTxtJ.String);
                obj.NumJumps=str2double(EditTxtNumJumps.String);
            end
            
            %Update GUI with objecte properties.
            function params2gui()
                EditTxtNumJumps.String=num2str(obj.NumJumps)
                EditTxtB.String=num2str(obj.B)
                EditTxtJ.String=num2str(obj.J)
            end
            
            %This function is called when we hit the 'start' pushbutton
            function start(~,~) %(~,~) is required to catch two inputs
                gui2params();
                obj.GuiGo=1;
                obj.runChain(1);
            end
            
            %This function is called when we hit the 'stop' pushbutton
            function stop(~,~) %(~,~) is required to catch two inputs
                obj.GuiGo=0;
            end
            
                      
        end
    end
    
end

