% ICE3:  T1 relaxation in repeated experiments
%
% At steady state, the magnetization after the nth pulse is
% equal to the magnetization after the (n+1)th pulse:
%
%     M = M0 - (M0 - Mss*cos(alpha))*e^(-t/T1)
%
%     where M0 is the magnitude of the equilibrium magnetization
%       and Mss is the steady-state magnetization
%
% The following script will plot the magnetization for
% three different T1 values when an experiment is repeated:
M0 = 49; % give them all the same proton density
alpha = 24; % degrees, the Ernst angle
TR = .1; % seconds
nRep = 50;
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
 % t = .01:.01:TR;
%M_GM= M0 - M0_GM*(1-cos(alpha_rad))*exp(-t/T1_GM);
%plot([0 t],[M0_GM M_GM])

alpha_rad = alpha*pi/180; % Matlab assumes you're giving your angles in radians
M0_GM = M0;
M0_WM = M0;
M0_CSF = M0; % start all at same intensity
for n = 1:nRep
    t = .01:.01:TR; % make a vector of time values, in 0.1s increments, from 0 to 5
    M_GM((n-1)*length(t)+(1:length(t))) = ...
        M0 - (M0 - M0_GM*cos(alpha_rad))*exp(-t/T1_GM);
    M0_GM = M_GM(n*length(t));
    M_WM((n-1)*length(t)+(1:length(t))) = ...
        M0 - (M0-M0_WM*cos(alpha_rad))*exp(-t/T1_WM);
    M0_WM = M_WM(n*length(t));
    M_CSF((n-1)*length(t)+(1:length(t))) = ...
        M0 - (M0 - M0_CSF*cos(alpha_rad))*exp(-t/T1_CSF);
    M0_CSF = M_CSF(n*length(t));
end

plot(.01:.01:TR*nRep,M_GM,'k'); hold on
plot(.01:.01:TR*nRep,M_WM,'g'); 
plot(.01:.01:TR*nRep,M_CSF,'r'); hold off 
legend('GM','WM','CSF')













