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

rf = [pi/4 pi/4];	% Complex sequence of RF rotations.
			% This will make more sense in later exercises.

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

xpos = [-20:.1:20];	% mm.
sig = 0*xpos;		% Allocate space.
grad = .1;		% G/cm.

% Starting magnetization.

for x=1:length(xpos)
    M = [0;0;1];
    dT = 2.3;	% ms between pulses.
    M = throt(abs(rf(1)),angle(rf(1))) * M;	% RF Rotation.

    for k = 2:length(rf)
	[A,B] = freeprecess(dT,T1,T2,0);	% No off-resonance.
	M = A*M+B;				% Propagate to next pulse.
    
	gradrot = 4258*2*pi*(dT/1000)*grad*(xpos(x)/10);  
						% Rotation from gradient.
	M = zrot(gradrot)*M;
    	M = throt(abs(rf(k)),angle(rf(k))) * M;	% RF Rotation.
    end;

    % --- These lines added from f1b ---
	
    [A,B] = freeprecess(dT,T1,T2,0);	% No off-resonance.
    M = A*M+B;				% Propagate .

    gradrot = 4258*2*pi*(dT/1000)*(-0.5*grad)*(xpos(x)/10);  
    						% Rotation from gradient.
    M = zrot(gradrot)*M;
    sig(x) = M(1)+i*M(2);

end;


subplot(2,1,1);
plot(xpos,abs(sig));
ylabel('Signal (fraction of M_0)');
xlabel('Position (mm)');
title('Signal Amplitude');
grid on;

subplot(2,1,2);
plot(xpos,angle(sig));
ylabel('Phase (rad)'); 
xlabel('Position (mm)');
title('Signal Phase');
grid on;







