%%% Homework on Bloch Eq. with rotation matrices
%%% We can compute the affects of an RF pulse (nutation), chemical shift (precession) and relaxation
%%%		using matrix operations.

%%% ---------------------------------------------------------------------------------------------
%%% The spin system is represented by a 3-element vector givng its x,y,z components.
%%% Start with magnetization at equilibrium (pointing along Z):
M0 = [0.0 ; 0.0; 1.0] 

%%% ---------------------------------------------------------------------------------------------
%%% An RF pulse is modeled as a rotation of "alpha" degrees about the X axis
alpha = 30.;					% degrees
Fx    = [ 1      0            0       ; 
	      0  cosd(alpha)  sind(alpha) ; 
	      0 -sind(alpha)  cosd(alpha)];

M1 = Fx * M0	% The magnetization after the RF pulse	

%%% ---------------------------------------------------------------------------------------------
%%% Chemical shift is just a rotation about the Z axis by an angle = delta_omega * time
d_omega = 17;					% this means our spins are 17 Hz off-resonance
TR      = 0.05;					% this is a 50 millisecond TR (or delay) time
theta   = d_omega * TR * 360;	% total precession angle (in degrees)
Pz = [ cosd(theta)  sind(theta)  0 ; 
	  -sind(theta)  cosd(theta)  0 ;
	       0            0        1]; 

M2 = Pz * M1	% The magnetization after TR seconds of precession	

%%% ---------------------------------------------------------------------------------------------
%%% Relaxation needs a matrix for the "loss" terms and a vector for "recovery"
T2 = 0.030;						% in seconds
T1 = 0.200;						% (ditto)
E2 = exp(-TR/T2);
E1 = exp(-TR/T1);
Rl = [ E2  0    0 ;				% the T1 & T2 "loss" operator
       0   E2   0 ;
	   0   0   E1];
Rr = [ 0 ; 0 ; 1-E1] ;			% the T1 "recovery" term

M3 = Rl * M2 + Rr	% The magnetization after TR seconds of relaxation

%%% ---------------------------------------------------------------------------------------------
%%% In Matlab, we can do all 3 steps in one line, like:
M4 = Rl*(Pz*(Fx*M0)) + Rr		% This is the same as M3 above


%%% ---------------------------------------------------------------------------------------------
%%% Homework problems
%%%		Use T2 = 40 msec, T1 = 300 msec, d_omega = 35 Hz for all problems
%%%		Start with M0 = [0.0 ; 0.0; 1.0] for all problems

%%% ---------------------------------------------------------------------------------------------
%%% Problem 1) Compute the M-vector for the sequence: "45x - 100msec" 
%%%			(this means RF pulse along X with alpha=45 degrees, then TR = 100 msec)
%%%			Also, compute the magnitude of the transverse component of M (i.e. |Mxy|)

%%% Problem 2) Compute the M-vector (and |Mxy|) after 10 repetitions of the same sequence

%%% Problem 3) Computer M for 50 repetitions of the sequence (i.e. 45x-100ms-45x-100ms-45x-...), and 
%%%				a) plot the magnitude of the M-vector versus repetition number 
%%%				b) describe what would happen if we kept repeating the sequence 1000 more times

%%% Problem 4) Repeat problem 3 with the RF pulse = 20x, and overplot the new M magnitude vs. rep

%%% Problem 5) Find the flip angle that maximizes the transverse magnetization after 50 pulses
%%%				Hint: test alpha = [1:90], compute sqrt(M(1)^2 + M(2)^2)
%%%				(If you do it right, you've found the "Ernst angle")

%%% Problem 6) So far, we've only modeled a single spin at a single frequency. Now model a
%%%			bunch of spins with a random distribution of frequencies. For the same "45x - 100msec"
%%%			sequence, compute the net M vector after one TR (like Problem #1), from a group of
%%%			200 spins with mean freq = 35 HZ and standard deviation = 10 Hz.
%%%			Matlab tips:
%%%				1) Use a "for" loop to call the code below 200 times
%%%				2) Each time compute M and add them all up
%%%					omega_mean  = 35.;
%%%					omega_sigma = 10.;
%%%					omega = omega_mean + randn() * omega_sigma;
%%%
%%%			What is the magnitude of the M-vector now, compared to Problem #1?

%%% Problem 7) Repeat problem 6, but calculate the time course of the net M-vector from after the
%%%				RF pulse until the end of TR. Plot the x-component of M (i.e. Mx) versus time. 
%%%				What happens if "omega_sigma" increases or decreases?
%%%			Matlab tips:
%%%				1) Add another loop over time, use something like:
%%%					t = 0:0.001:TR;
%%%					nt = size(t,2);
%%%					nspins = 200;
%%%					for i = 1:nspins
%%%						omega = omega_mean + randn() * omega_sigma;
%%%						for j = 1:nt
%%%							theta = omega * t(j) * 360;	
%%%								.
%%%								.
%%%								.
%%%						end
%%%					end

%%% Bonus questions:
%%%			In the examples, we computed precession before we computed relaxation, even though they 
%%%			happen at the same time. 
%%%			A) What would happen if we computed relaxation before precession?
%%%			B) Why? (i.e. what property do Pz and R have that mean this in mathematical terms?)
%%%			C) What assumption did we make about RF pulses such that we ignore relaxation & 
%%%			precession during them?