clear all; close all;
global H gamma lambda s nu
H = [0,0,1]; %external field
gamma = 1e5; %gyromagnetic ratio
alpha = .1; %damping factor
lambda = -alpha*gamma/1 ;

N=5; %number of moments
M=zeros(N,N,N,3); %array of NxNxN moment vectors

dt=1e-5;

%heat equation constants
k=.001;
dx=.001;
eta=k*dt/dx^2;


%initial random moment vectors
for i=1:N
    for j=1:N
        for k=1:N
            m = [rand()-.5,rand()-.5,rand()-.5];
            for l=1:3 
                M(i,j,k,l) = m(l)/norm(m);
            end
        end
    end
end

%initial temp distribution
T=zeros(3*N,3*N,3*N);
for i=N+1:2*N
    for j=N+1:2*N
        for k=N+1:2*N
            T(i,j,k)=1;
        end
    end
end

x=1:N;
[X,Y,Z]=meshgrid(x,x,x);

C=colormap(jet(128));

t=0;
maxit=300;
TM=zeros(1,maxit); %store magnetization to plot
tt=dt:dt:maxit*dt; 

for it=1:maxit
    totmag=0;
    for i=1:N
        for j=1:N
            for k=1:N
                %evolve each moment vector according to
                %Landau-Lifshitz, using RK4
                s=[rand()-.5,rand()-.5,rand()-.5]; %random torque vector
                s=s/norm(s);
                nu = T(i+N,j+N,k+N); %random torque magnitude given by temp
                m  = [M(i,j,k,1)  M(i,j,k,2) M(i,j,k,3)];
                f1 = LL( m );
                f2 = LL( m + f1 .* (dt/2));
                f3 = LL( m + f2 .* (dt/2));
                f4 = LL( m + f3 .*  dt   );
                m  = m + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
                m  = m/norm(m);
                M(i,j,k,:) = m;
                totmag=totmag+m; 
                
                %plot moment vector, with color given by temperature
                quiver3(i,j,k,m(1)/2,m(2)/2,m(3)/2,'LineWidth',2,'maxheadsize',1,...
                        'Color',C(round(T(i+N,j+N,k+N)*115),:));
                axis([0 N+1 0 N+1 0 N+1]);
                hold on

            end
        end
    end
    
    drawnow
    hold off
    
    %update temp
    for i=2:3*N-1
        for j=2:3*N-1
            for k=2:3*N-1
                T(i,j,k)=eta*(T(i-1,j,k)+T(i+1,j,k)+T(i,j-1,k)+T(i,j+1,k)+T(i,j,k-1)+T(i,j,k+1)...
                         -6*T(i,j,k)) +T(i,j,k) ;
            end
        end
    end
    
    t=t+dt;
    TM(it)=norm(totmag)/N^3;
    
end

figure
plot(tt,TM)
xlabel('t')
ylabel('Magnetization (normalized)')
            