zadanie10.m

clear;
page_output_immediately(1);
% obliczanie energii lancucha
function E = energy(chain, J, H)
	E = -sum(chain)*H;
	E += -J*sum(chain.*[chain(end) ; chain(1:end-1)]);
endfunction

% generator LCG
function value = LCGrand()
	choice = 1;

	global state;
	% parametry generatora LCG
	a = [69069 ; 16807 ; 1664525 ; 	1103515245];
	c = [1 ; 0 ; 1013904223 ; 	12345];
	m = [ power(2,32); power(2,32) - 1; power(2,32); power(2,32)];
	state = mod(a(choice)*state+c(choice),m(choice));
	value = state/(m(choice)-1);
endfunction

function [E,M] = generateChain(H,J,T,N,iterations, steps, generator)
	E = 0;
	M = 0;
	%lancuch losowych spinow
	chain = round(rand(N,1)).*2-1;
	for cnt1=1:steps
		for cnt2=1:iterations
			% losowanie nru odwracanego spinu
			i = round(generator()*(N-1)+1);
			% odwrocenie pojedynczego spinu
			trial_chain = chain;
			trial_chain(i) *= -1;
			% roznica energii pomiedzy nowym i starym lancuchem
			dE = energy(trial_chain, J, H) - energy(chain, J, H);
			if(dE <= 0)
				% jezeli nizsza energia nowego lancucha - odwracamy
				chain = trial_chain;
			else
				if(exp(-dE/T) >= generator())
					% jezeli energia jest wyzsza
					% to z p-bstwem w/w przyjmujemy nowy lancuch
					chain = trial_chain;
				endif
			endif
		endfor
		E += energy(chain,J,H)/steps;
		M += abs(sum(chain))/N/steps;
	endfor
endfunction

%{---------------------------------------------------------------%}
imgParam = '-S1000,800';
imgExt = 'png';
% ilosc krokow Monte Carlo
MCsteps = 50;
% ilosc iteracji do termalizacji ukladu
iterations = 1500;
% ilosc spinow
N = 100;
T = 0.0001:0.1:12;
%%
% ziarno generatora LCG
global state = 1;
count = length(T);
E = zeros(count,1);
M = zeros(count,1);

disp(':> Uruchomiono generowanie wykresow');
t0 = cputime();
%%
H = 0;
J = 1;
for i=1:count
	[E(i), M(i)] = generateChain(H,J,T(i),N,iterations,MCsteps,@rand);
endfor
disp(':> Zakonczono generowanie ukladu 1.');
disp([':> Czas generacji: ' num2str((t1 = cputime()) - t0)]);
figure(1);
plot(T,E)
xlabel('Temperatura');
ylabel('Energia');
title(['Energia w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator rand()']);
grid;
print(['EJ' num2str(J) 'H' num2str(H) 'rand.' imgExt], imgParam);
close(1);


figure(2);
plot(T,M)
xlabel('Temperatura');
ylabel('Srednia magnetyzacja');
title(['Srednia magnetyzacja w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator rand()']);
grid;
print(['MJ' num2str(J) 'H' num2str(H) 'rand.' imgExt], imgParam);
close(2);

H = 0;
J = -1;
for i=1:count
	[E(i), M(i)] = generateChain(H,J,T(i),N,iterations,MCsteps,@rand);
endfor
disp(':> Zakonczono generowanie ukladu 2.');
disp([':> Czas generacji: ' num2str((t2 = cputime()) - t1)]);
figure(3);
plot(T,E)
xlabel('Temperatura');
ylabel('Energia');
title(['Energia w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator rand()']);
grid;
print(['EJ' num2str(J) 'H' num2str(H) 'rand.' imgExt], imgParam);
close(3);


figure(4);
plot(T,M)
xlabel('Temperatura');
ylabel('Srednia magnetyzacja');
title(['Srednia magnetyzacja w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator rand()']);
grid;
print(['MJ' num2str(J) 'H' num2str(H) 'rand.' imgExt], imgParam);
close(4);


H = 0;
J = 1;
for i=1:count
	[E(i), M(i)] = generateChain(H,J,T(i),N,iterations,MCsteps,@LCGrand);
endfor
disp(':> Zakonczono generowanie ukladu 3.');
disp([':> Czas generacji: ' num2str((t3 = cputime()) - t2)]);
figure(5);
plot(T,E)
xlabel('Temperatura');
ylabel('Energia');
title(['Energia w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator LCGrand()']);
grid;
print(['EJ' num2str(J) 'H' num2str(H) 'LCGrand.' imgExt], imgParam);
close(5);



figure(6);
plot(T,M)
xlabel('Temperatura');
ylabel('Srednia magnetyzacja');
title(['Srednia magnetyzacja w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator LCGrand()']);
grid;
print(['MJ' num2str(J) 'H' num2str(H) 'LCGrand.' imgExt], imgParam);
close(6);

H = 0;
J = -1;
for i=1:count
	[E(i), M(i)] = generateChain(H,J,T(i),N,iterations,MCsteps,@LCGrand);
endfor
disp(':> Zakonczono generowanie ukladu 4.');
disp([':> Czas generacji: ' num2str((t4 = cputime()) - t3)]);
figure(7);
plot(T,E)
xlabel('Temperatura');
ylabel('Energia');
title(['Energia w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator LCGrand()']);
grid;
print(['EJ' num2str(J) 'H' num2str(H) 'LCGrand.' imgExt], imgParam);
close(7);

figure(8);
plot(T,M)
xlabel('Temperatura');
ylabel('Srednia magnetyzacja');
title(['Srednia magnetyzacja w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator LCGrand()']);
grid;
print(['MJ' num2str(J) 'H' num2str(H) 'LCGrand.' imgExt], imgParam);
close(8);

disp([':> Calkowity czas generacji: ' num2str(tt = cputime() - t0)]);
disp([':> Pojedynczy krok MC: ' num2str(tt/4/MCsteps/count)]);