% Bloch Equation Simulation, Excercise F-1a
% -----------------------------------------
% 


t = 0:.04:6;           % in ms.
dT = t(2)-t(1);         % in ms.
b1 = .05*sinc(t-3);     % in G.
rf = 2*pi*4258*b1*(dT/1000); % Rotation in radians.

T1 = 600;       % ms.
T2 = 100;       % ms.

freq = [-1000:20:1000];	% Hz.
sig = 0*freq;		% Allocate space.

% Starting magnetization.

for f=1:length(freq)
    M = [0;0;1];
    M = throt(abs(rf(1)),angle(rf(1))) * M;	% RF Rotation.

    for k = 2:length(rf)
	[A,B] = freeprecess(dT,T1,T2,freq(f));
	M = A*M+B;				% Propagate to next pulse.
    	M = throt(abs(rf(k)),angle(rf(k))) * M;	% RF Rotation.
    end;
    sig(f) = M(1)+i*M(2);

end;

subplot(2,1,1);
plot(freq,abs(sig));
ylabel('Signal (fraction of M_0)');
xlabel('Frequency (Hz)');
title('Signal Amplitude');
grid on;

subplot(2,1,2);
plot(freq,angle(sig));
ylabel('Phase (rad)'); 
xlabel('Frequency (Hz)');
title('Signal Phase');
grid on;