function isingP1
% modelo de ISING with periodic BC

% N^2 = num de espines
h1 = figure;
hp10 = uipanel('parent',h1,'pos',[0 0 0.15 1]);
hp11 = uipanel('parent',h1,'pos',[0.15 0.5 0.85 0.5]);
hp12 = uipanel('parent',h1,'pos',[0.15 0 0.85 0.5]);

ha11 = axes('parent', hp11,'pos',[0.1 0.1 0.87 0.87]);
ylabel('energía','parent',ha11)
ha12 = axes('parent',hp12,'pos', [0.1 0.1 0.87 0.87]);
ylabel('momento','parent',ha12)
h= figure;
set(h1,'pos',get(h,'pos')+[100 100 0 0])
set(h,'closerequestfcn',@my_closefcn)
hp1 = uipanel('parent',h,'pos',[0 0 0.2 1]);
hp2 = uipanel('parent',h,'pos',[0.2 0 0.8 1]);
ha = axes('parent', hp2);
colormap(gray(2))
axis square
N = 10;
ss = [];
H = 0 ; %campo
Hmax = 3;
mu = 1;
J = 1;
T = 1;  % T en unidades de energía
ngo = 10;
%timecalc = 0.001;
refre = 0.5;
d = [];
s = [];
ener = [];
ener0 = 0;
mom = [];
mom0 = 0;
inicializa



% controls ***********************
uicontrol('style','togglebutton','parent',hp1,...
    'pos',[10,5,40,25],...
    'max',1,'min',0,'value',1,...
    'string','go',...
    'callback',@avan)
    function avan(src,~)
        st = get(src,'value');
        switch st
            case 0
                %start(tim)
                set(src,'string','stop')
                set(hepaso,'enable','off')
                set(herefre,'enable','off')
                set(hen,'enable','off')    
                go
                start(tim2)
            case 1
                %stop(tim)
                set(src,'string','go')
                set(hepaso,'enable','on')
                set(herefre,'enable','on')
                set(hen,'enable','on')
                goTAG = 0;
                stop(tim2)
        end
    end
uicontrol('style','pushbutton','parent',hp1,...
    'pos',[60,5,40,25],...
    'string','show',...
    'callback',@para)
    function para(~,~)
        va = strcmpi(get(tim2,'running'),'on');
        if va, stop(tim2), end
        plot(ener,'parent',ha11)
        plot(mom,'parent',ha12)
        ylabel('energía','parent',ha11)
        ylabel('momento','parent',ha12)
        set(ht10e,'string',['Eav=' num2str(mean(ener))])
        set(ht10e1,'string',['std=' num2str(std(ener))])
        set(ht10m,'string',['Mav=' num2str(mean(mom))])
        set(ht10m1,'string',['std=' num2str(std(mom))])
        if va,  start(tim2), end
    end
    
% reinicializa
uicontrol('style','pushbutton','parent',hp1,...
    'pos',[70,100,40,25],...
    'string','ini',...
    'callback',@bini)
    function bini(~,~)
        inicializa
    end

uicontrol('style','pushbutton','parent',hp1,...
    'pos',[20,100,40,25],...
    'string','cl E',...
    'callback',@bclE)
    function bclE(~,~)
        ener = [];
        mom = [];
    end



% cambia Temperatura
uicontrol('style','edit','parent',hp1,...
    'pos',[50 130 60 20],...
    'string',num2str(T),...
    'callback',@eT)
    function eT(src,~)
        T = str2num(get(src,'string'));
    end
uicontrol('style','text','parent',hp1,...
    'pos',[10 130 40 20],...
    'string','T')

% cambia Número de espines
hen = uicontrol('style','edit','parent',hp1,...
    'pos',[50 155 60 20],...
    'string',num2str(N),...    
    'callback',@eN);
    function eN(src,~)
        N = str2num(get(src,'string'));
        inicializa
    end
uicontrol('style','text','parent',hp1,...
    'pos',[10 155 40 20],...
    'string','N(>1)')

% cambia interacción
uicontrol('style','edit','parent',hp1,...
    'pos',[50 180 60 20],...
    'string',num2str(J),...
    'callback',@eJ)
    function eJ(src,~)
        J = str2num(get(src,'string'));
    end
uicontrol('style','text','parent',hp1,...
    'pos',[10 180 40 20],...
    'string','J')

%cambia Hmax
uicontrol('style','edit','parent',hp1,...
    'pos',[50 360 40 20],...
    'string',num2str(Hmax),...
    'callback',@eHmax)
    function eHmax(src,~)
        Hmax = str2num(get(src,'string'));
        set(hsl,'max',Hmax);
        set(hsl,'min',-Hmax);
    end
uicontrol('style','text','parent',hp1,...
    'pos',[50 380 40 20],...
    'string','Hmax')
% control de H
uicontrol('style','pushbutton','parent',hp1,...
    'pos',[50 280 40 20],...
    'string','ZERO',...
    'callback',@bzero);
    function bzero(~,~)
        H = 0;
        set(hsl,'value',H)
        set(ht,'string',['H=' num2str(H)])
    end
hsl = uicontrol('style','slider','parent',hp1,...
    'pos',[20 210 20 200],...
    'max',Hmax,'min',-Hmax,...
    'value',H,...
    'callback',@slH);
    function slH(src,~)
        H = get(src,'value');
        set(ht,'string',['H=' num2str(H)]) 
        ener0 = energia;
        mom0 = momento;
    end
ht = uicontrol('style','text','parent',hp1,...
    'pos',[40 300 80 40],...
    'string',['H=' num2str(H)]);

% refresco de pintaS
herefre = uicontrol('style','edit','parent',hp1,...
    'pos',[50 35 60 20],...
    'string',num2str(refre),...
    'callback',@eRefre);
    function eRefre(src,~)
        refre = str2num(get(src,'string'));
        set(tim2,'period',refre);
    end
uicontrol('style','text','parent',hp1,...
    'pos',[10 35 40 20],...
    'string','refre')
% paso de tiempo
hepaso = uicontrol('style','edit','parent',hp1,...
    'pos',[50 60 60 20],...
    'string',num2str(ngo));
    
uicontrol('style','text','parent',hp1,...
    'pos',[10 60 40 20],...
    'string','ciclos')
htiempo = uicontrol('style','text','parent',hp1,...
    'pos',[10 80 100 20],...
    'string','');
ht10e = uicontrol('style','text','parent',hp10,...
    'pos',[10 330 85 20],...
    'string','ee');
ht10e1 = uicontrol('style','text','parent',hp10,...
    'pos',[10 300 85 20],...
    'string','ee');
ht10m = uicontrol('style','text','parent',hp10,...
    'pos',[10 130 85 20],...
    'string','mm');
ht10m1 = uicontrol('style','text','parent',hp10,...
    'pos',[10 100 85 20],...
    'string','mm');


 %tim = timer('TimerFcn',@timm,...
 %   'period',timecalc,...
 %   'executionmode','fixedDelay');         
 %   function timm(~,~)
 %       tic
 %       [den,dmo] = flipS;
 %       ener0 = ener0+den;
        %disp([num2str(ener0) '   ' num2str(energia)])
 %       mom0 = mom0+dmo;
        %disp([num2str(mom0) '   ' num2str(momento)])
        %disp('   ')
 %       ener = [ener ener0]; 
 %       mom = [mom mom0];
 %       toc
 %   end

tim2 = timer('TimerFcn',@timm2,...
    'period',refre,...
    'executionmode','fixedDelay'); 
    function timm2(~,~)
        pintaS
        go
    end


% uicontrol('style','pushbutton','parent',hp1,...
%     'pos',[10,440,40,25],...
%     'string','special',...
%     'callback',@espe)
%     function espe(~,~)
%         xs = xcorr2(s);
%         figure(111)
%         image(xs)
%         colormap(gray(max(xs(:))))
%     end

%***********************************************



% funciones *********************************
    function go
        ngo = str2num(get(hepaso,'string'));
        tic
        for kgo = 1:ngo
            [den,dmo] = flipS;
            ener0 = ener0+den;
            %disp([num2str(ener0) '   ' num2str(energia)])
            mom0 = mom0+dmo;
            %disp([num2str(mom0) '   ' num2str(momento)])
            %disp('   ')
            ener = [ener ener0]; 
            mom = [mom mom0];
            
        end
        tfin = toc;
        set(htiempo,'string',['dt = ' num2str(tfin)])
    end  
            

    function pintaS
        image((s+1)/2+1,'parent',ha)
        axis off
        drawnow
    end

    function en = energia
        % calcula la energía total de una configuración
        a1 = circshift(s,1);
        a2 = circshift(s,-1);
        a3 = circshift(s,[0 1]);
        a4 = circshift(s,[0 -1]);

        aa = s.*(a1+a2+a3+a4);
        ss =  sum(aa(:))/2;
        e0 = - H*mu*sum(s(:));
        ei = - J*ss;
        en = e0 + ei;
    end
    function mo = momento
        mo = mu*sum(s(:));
    end

    function [den,dmo] = flipS
        % cambia aleatoriamente un solo espín y calcula el cambio de
        % energía
        ind = randi(N^2);
        switch d(ind)
            case 0
                v = ind+[1 -1 N -N];
            case 1
                v = ind+[1 N-1 N -N];
            case 2
                v = ind+[-1 -N+1 N -N];
            case 3
                v = ind+[1 -1 N  N*(N-1)];
            case 4
                v = ind+[1 -1 -N*(N-1) -N];
            case 13
                v = ind+[1 N-1 N  N*(N-1)];
            case 14
                v = ind+[1 N-1 -N*(N-1) -N];
            case 23
                v = ind+[-1 -N+1 N  N*(N-1)];
            case 24
                v = ind+[-1 -N+1  -N*(N-1) -N];
        end
        den = 2*s(ind)*( H*mu + J*sum(s(v)) );
        if den <= 0
            dmo = -2*mu*s(ind); s(ind) = -s(ind);
        else
            x = rand;
            if T==0, c=0; else  c = exp(-den/T); end
            if c >= x, dmo=-2*mu*s(ind); s(ind) = -s(ind); else den=0; dmo=0; end
        end

    end


    function inicializa
        d = zeros(N);
        d(1,:) = 1;
        d(N,:) = 2;
        d(:,1) = 3;
        d(:,N) = 4;
        d(1,1) = 13;
        d(1,N) = 14;
        d(N,1) = 23;
        d(N,N) = 24;
        s = 2*randi(2,N)-3; %matriz de espines
        %s = ones(N);
        pintaS
        ener = [];
        mom = [];
        ener0 = energia;
        disp(ener0)
        ener = [ener ener0];
        mom0 = momento;
        disp(mom0)
        mom = [mom mom0];
    end

   function my_closefcn(~,~)
        %delete(tim)
        delete(tim2)
        if ishandle(h1), delete(h1), end
        delete(gcf)
     end

end
