N = 500; % number of spins
M = 1000; % number of time-point to simulate
SE = 200; % sample interval between spin-echoes 
D = 0; % diffusion coefficient

dt = 4e-6; % time interval for bloch simulator
mx = ones(N,1); % Initial magnetization is mxy = 1x;
my = zeros(N,1);
mz = zeros(N,1);
dp = (rand(N,2)-0.5)*1.8; % generate random spin distribution

G = [4,0]; % set gradient to 4 G/cm

S = zeros(M,1);

for n=1:M
    
    [mx,my,mz] = bloch(0,G,dt,600e-3,40e-3,0,dp,2,mx,my,mz); % simulate one step
    dp = dp + randn(N,2)*D; % move spins randomly
    S(n) = mean(mx(:) + i*my(:)); % compute the FID
    
    if mod(n,SE)==SE/2 % simulate a SE by flipping phase
        mxy = mx(:) + i*my(:);
        mxy = abs(mxy).*exp(i*(pi-angle(mxy)));
        mx = real(mxy(:));
        my = imag(mxy(:));
    end
    
    figure(100), subplot(211),quiver(dp(:,1),dp(:,2),mx(:),my(:)); axis([-1,1,-1,1]);
    subplot(212), plot(abs(S)); drawnow;
end
