% Ising Monte Carlo simulation: analyse 
% Version 08.03.2007
% 
% Reads and analyses the data produced by the function ising.
% Not meant to be used as a standalone program
% 

function[]=analyse(num,THERM)

%clf(1);clf(5);clf(6)
for index=1:num
    %
    % file read 
    %    
    eval(['load tulokset' int2str(index)]);
    %
    % read corresponding set of parameters
    % 
    eval(['load parametrit' int2str(index)]);
    %
    % evaluate averages
    % 
    Magnetization1(index)=mean(ABS_M(MCSTEPS-THERM:MCSTEPS));
    Magnetization(index)=mean(M(MCSTEPS-THERM:MCSTEPS));
    Energy(index)= mean(E(MCSTEPS-THERM:MCSTEPS));
    Energy2(index)=mean(E2(MCSTEPS-THERM:MCSTEPS));
    if(abs(J)<0.0001)
        % set dummy mu, or new Matlab assumes I mean function mu()  
        mu = 2; 	            
        X(index)=B(ib)/mu; 
        TperJ(index)=T(it); 
    else
        TperJ(index)=abs(T(it)/J); 
    end
    Ball(index)=B(ib);
end
% energy statistical error
Energyerror = sqrt(Energy2-Energy.^2)  ;

if(abs(J)<0.0001)
    % ideal paramagnet
    %
    % exact M = mu*tanh(mu*B/kB*T) =now= tanh(x)
    %
    Nx = 300 ;
    dx = (X(num)-X(1))/(Nx-1) ;
    for ix = 1:Nx
        xtab(ix) = X(1) + (ix-1)*dx;
    end
    % 
    % plot magnetization vs. x
    %
    figure(2); clf(2)
    plot(xtab,tanh(xtab),'r ');hold on; % exact result
    plot(X,Magnetization,'*');hold on;
    %plot(X,Magnetization1,'o');hold on;;
    xlabel('x = \mu B/(K_B T)');ylabel('M(x)')
    line2=[,'L=',L,'  J=',num2str(J)];
    title({['Ideal paramagnet: Magnetization compared with exact ' ...
           'result tanh(x) (red line)'];line2})
    %
    % Plot energy vs. x
    % -----------------------------

    
    figure(3);clf(3)
    errorbar(X,Energy,Energyerror,'*');hold on;
    xlabel('x = \mu B/(K_B T)');ylabel('Energy')
    title({['Ideal Paramagnet Energy, compared with the exact '...
            'result -mu*x*tanh(x) (red line)'];line2})
    plot(xtab,-mu*xtab.*tanh(xtab),'r ');hold on; % exact result
else
    %
    % ferromagnet or paramagnet
    % 
    if(J>0)
        %
        % Onsager's exact magnetization
        % T<Tc : m = (1-(sinh(2*beta*J))^-4 )^(1/8)
        % T>Tc : m = 0
        %
        Tcrit=2*J/log(1+sqrt(2));
        for  itemp = 1:5000  % temperature
            temp(itemp) = 0.01 + (itemp-1)*0.001;
            if(temp(itemp)<Tcrit) 
                apu1=sinh(2*J/temp(itemp));
                K=1/(apu1^2);
                m(itemp)=(1-K^2).^(1/8);
            else
                m(itemp) = 0;
            end
        end
    end
    if(size(B,2)<5) 
        line2=['L=',num2str(L),'   J=',num2str(J),'   B=',num2str(B)];
    else
        line2=['L=',num2str(L),'   J=',num2str(J),'   B=', ...
               num2str(B(1)),'->',num2str(B(size(B,2)))];
    end
    %
    % Plot magnetization vs. T/J
    % --------------------------- 
    figure(2); clf(2)
    if(J>0) 
        plot(temp,m,'r ');hold on;
    end
    plot(TperJ,Magnetization,'*');hold on;
    plot(TperJ,Magnetization1,'o');hold on;
    if(J>0) 
        %figure(2);plot([Tcrit Tcrit],[0 1],'g');hold on;
        title({'2D Ising model Magnetization';['compared with exact ' ...
               'B=0 result in the thermodynamical limit (red line, ' ...
                            'spin-up side)'];line2})
        ylabel('M,|M|')
    else
        title({'2D Ising model: Magnetization';line2})
        ylabel('M,|M|')
    end
    xlabel('|T/J|')
    %
    % Plot energy and specific heat  
    % -----------------------------
    
    figure(3);clf(3)
    if(J>0) 
        plot([Tcrit Tcrit],[-2 0],'g');hold on;
    end
    %plot(TperJ,Energy,'*');hold on;
    errorbar(TperJ,Energy,Energyerror,'*');hold on;
    xlabel('|T/J|');ylabel('Energy')
    if(J>0) 
        title({'2D Ising model Energy';['The green vertical line shows ' ...
               'the location of the exact T_c in the thermodynamical limit'];line2}) 
    else
        title({'2D Ising model: Energy';line2})
    end
    % C from fluctuations
    tmp = Energy2-Energy.^2 ;
    tmp2 = TperJ.^2;    
    C = tmp./tmp2*L*L;
    figure(4);clf(4);plot(TperJ,C,'o');hold on;
    line3=['circles: from energy fluctuation'];
    % C from derivative of E
    if(size(TperJ,2)>1)
        dT = TperJ(2)-TperJ(1);
        dT = J*dT ; 
        for i = 1: size(TperJ,2)-1
            Cder(i)     = (Energy(i+1)-Energy(i))/dT;
            TperJder(i) = TperJ(i)+dT/2;
        end
        plot(TperJder,Cder,'x');hold on;
        line3=['circles: from energy fluctuation     crosses: from ' ...
               'numerical derivative dE/dT'];
    end
    xlabel('|T/J|');ylabel('C')
    title({'2D Ising model Specific heat';line2;line3}) 
    %
    % M vs. B
    %
    if(size(B,2)>1 & size(T,2)==1)
        figure(5);clf(5);plot(B,Magnetization,'o');hold on;
        plot([B(1) B(size(B,2))],[0 0],'g')
        plot([0 0],[-1 1],'g')
        xlabel('B')
        ylabel('M')
        title(['2D Ising model T=',num2str(T),' J=',num2str(J)])
    end
end

% 
% data file output
% 

% $$$ i = 0;
% $$$ for it=1:size(T,2)
% $$$     for ib=1:size(B,2)
% $$$         i = i + 1;
% $$$         outdat(i,1) = T(it)/J ;
% $$$         outdat(i,2) = B(ib) ;
% $$$         outdat(i,3) = Energy(i) ;
% $$$         outdat(i,4) = Magnetization(i) ;
% $$$         outdat(i,5) = Magnetization1(i) ;
% $$$     end
% $$$ end

outdat(:,1) = TperJ.' ; % has to take transpose
outdat(:,2) = Ball.';
outdat(:,3) = Energy.';
outdat(:,4) = Magnetization.';
outdat(:,5) = Magnetization1.';

save results.dat -ascii outdat 
%
% delete temporary data files
%
cleanup(1) ;
disp('All done')

