% WA3:  T1 relaxation and Ernst angle
%
% Problem 1.  If I apply an excitation pulse to a sample, the equilibrium
% magnetization is rotated alpha degrees, where alpha is the flip angle.
% If I apply no other RF pulses, the magnetization will relax back
% to equilibrium with a time constant T1, according to the equation:
%
%     M = M0 - (M0 - M0*cos(alpha))*e^(-t/T1)
%
%     where M0 is the magnitude of the equilibrium magnetization
%
% The following script will plot the recovery of the magnetization for
% three different T1 values:

t = 0:.1:5; % make a vector of time values, in 0.1s increments, from 0 to 5
T1_GM = 1.1; % T1 of gray matter at 3T is ~1100ms
T1_WM = 0.8; % T1 of white matter is slightly shorter
T1_CSF = 3; % T1 of CSF is always the longest

% assume that all three tissue types have equilibrium magnetization, M0,
% equal to 49 (arbitrary units).
M0 = 49; 
alpha = 60; % degrees
alpha_rad = alpha*pi/180; % Matlab assumes you're giving your angles in radians
plot(t,M0 - M0*(1-cos(alpha_rad))*exp(-t/T1_GM),'k'); hold on
plot(t,M0 - M0*(1-cos(alpha_rad))*exp(-t/T1_WM),'g'); 
plot(t,M0 - M0*(1-cos(alpha_rad))*exp(-t/T1_CSF),'r'); 
plot([0 5],[0 0],'k:'); hold off % add a dashed line to show M0 = 0
legend('GM','WM','CSF')

% A) What is the magnitude of the longitudinal magnetization at t = 1.0s
% for each of the three "tissue" types? (Hint: just replace t and T1 in the
% following equation with the appropriate value.  Edit the script to assign 
% these values to 3 different variables (e.g. M_GM, M_WM and M_CSF) because
% you'll use these later.)

M = M0 - M0*(1-cos(alpha_rad))*exp(-t/T1)


% B) If you apply a 2nd 60 degree pulse at time t = 1s, what would be the 
% magnitude of the transverse magnetization for GM, WM, and CSF?  
% (Hint:  Mt = M*sin(alpha), where Mt is the transverse
% mag., M is the longitudinal mag., and alpha is the flip angle, e.g.



% C) Plot the transverse magntization decay curves according to the form:
%         Mt = Mt_0*exp(-t/T2)
% where Mt_0 for each tissue type was calculated in B), and the different 
% decay constants are:
t = 0:.001:.1; % note the much shorter time scale!
T2_GM = .050;  T2_WM = .040;  T2_CSF = 1.000;

plot(t,Mt_GM*exp(-t/T2_GM),'k'); hold on
plot(t,Mt_WM*exp(-t/T2_WM),'g'); 
plot(t,Mt_CSF*exp(-t/T2_CSF),'r'); hold off

% D) Copy the first plotting routine, and re-plot the longitudinal "decay" 
% curves with a different M0 for each tissue type, according to your answer 
% in A).




% Problem 2:  Write a function that, when passed values for TR and T1
% returns the Ernst angle (in degrees), using the given equation:
%   cos(alpha_ernst) = exp(-TR/T1)
% Hints: the Matlab function for inverse cosine is acos(); don't forget to
% convert from radians to degrees.


















