%% HW3 Problem 1 Solution. Eric Aliotta. 01.25.2017

% 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]);