% set a range of values to simulate
TR = 10:10:10000;
TE = 5:0.1:100;
alpha = 30*pi/180; % fip angle, deg
% Tissue A
T1a = 2000; % T1, ms
T2a = 40; % T2, ms
T2as = 25; % T2*, ms
% Tissue B
T1b = 500; % T1, ms
T2b = 40; % T2, ms
T2bs = 25; % T2*, ms
A_GRE = zeros(length(TR),length(TE),2);
A_SE = zeros(length(TR),length(TE),2);
% cycle through each TE and TR
for j = 1:length(TR)
for k = 1:length(TE)
TRtmp = TR(j);
TEtmp = TE(k);
% GRE contrast for Tissue A
A_GRE(j,k,1) = ( 1-exp(-TRtmp/T1a) )*( sin(alpha)*exp(-TEtmp/T2as) )/( 1-cos(alpha)*exp(-TRtmp/T1a) );
% GRE contrast for Tissue B
A_GRE(j,k,2) = ( 1-exp(-TRtmp/T1b) )*( sin(alpha)*exp(-TEtmp/T2bs) )/( 1-cos(alpha)*exp(-TRtmp/T1b) );
% SE contrast for Tissue A
A_SE(j,k,1) = ( 1-exp(-TRtmp/T1a) )*exp(-TEtmp/T2a);
% SE contrast for Tissue B
A_SE(j,k,2) = ( 1-exp(-TRtmp/T1b) )*exp(-TEtmp/T2b);
end
end
% find the point of maximum contrast for both SE and GRE
TRsearch = repmat(TR',[1,length(TE)]);
TEsearch = repmat(TE,[length(TR),1]);
C_GRE = abs(A_GRE(:,:,1)-A_GRE(:,:,2));
C_SE = abs(A_SE(:,:,1)-A_SE(:,:,2));
Cmax_GRE = max(C_GRE(:));
C_GRE(C_GRE<Cmax_GRE) = 0;
TR_opt_GRE = TRsearch(C_GRE~=0);
TE_opt_GRE = TEsearch(C_GRE~=0);
Cmax_SE = max(C_SE(:));
C_SE(C_SE<Cmax_SE) = 0;
TR_opt_SE = TRsearch(C_SE~=0);
TE_opt_SE = TEsearch(C_SE~=0);
fprintf('Gradient Echo : Optimal TE = %2.1fms\n',TE_opt_GRE);
fprintf(' Optimal TR = %2.1fms\n\n',TR_opt_GRE);
fprintf('Spin Echo : Optimal TE = %2.1fms\n',TE_opt_SE);
fprintf(' Optimal TR = %2.1fms\n\n',TR_opt_SE);
fprintf('For optimal T1 contrast, the spin echo sequence is %2.2fx longer! \n',TR_opt_SE/TR_opt_GRE);
%% view plots of contrast for each parameter
figure;
subplot(1,2,1);
surf(TE,TR,abs(A_GRE(:,:,1)-A_GRE(:,:,2)),'LineStyle','none'); view(0,90); title('GRE Contrast');
xlabel('TE [ms]'); ylabel('TR [ms]'); zlabel('Contrast'); ylim([0 3000]);
subplot(1,2,2);
surf(TE,TR,abs(A_SE(:,:,1)-A_SE(:,:,2)),'LineStyle','none'); view(0,90); title('SE Contrast');
xlabel('TE [ms]'); ylabel('TR [ms]'); zlabel('Contrast'); ylim([0 3000]);