
function output = mt_fit(x, m, W_B, W_F, delta, R1, t_m, t_s, t_r);
% This is the function to which the parameters are passed during MT fit
% All parameters are named exactly as they appear in the paper.
k = x(1); % the k-parameter is passed through x
f = x(2); % the f-parameter is passed through x

C = [cos(10*pi/180) 0;
	0 1];
T2_B = 11e-6;

R1_B = 1;
R1_F = R1 - k*(R1_B - R1)/(R1_B - R1 + k*(1-f)/f);

R = [(-R1_F - k), (k*(1-f)/f);
		(k), (-R1_B - k*(1-f)/f)];

A = R1_F*R1_B + R1_F*k*(1-f)/f + R1_B*k;
E_s = exp(R*t_s);
E_r = exp(R*t_r);

I = eye(2); %identity matrix

%Find the Mz for each offset frequency value
for i = 1:length(W_B);
    W = [W_F(i), 0;
          0,   W_B(i)];
    E_m = exp((R + W)*t_m); 
    D = A + (R1_F + k)*W_B(i);
    M_eq = [1-f f]';
    M_ss = 1/D*[(1-f)*(A + R1_F*W_B(i)) f*(A + R1_B*W_F(i))]';

    term_1 = inv(I - E_s*E_m*E_r*C);
    term_2 = (E_s*E_m*(I - E_r) + I-E_s)*M_eq;
    term_3 = E_s*(I-E_m)*M_ss;
    
    M_z(i,:) = term_1*(term_2 + term_3);
end

%The measurements are normalized with S0, the theoretical estimate should also be normalized
%but with another theoretical Mz at zero offset frequency.
E_m = exp((R)*t_m);
D = A ;
M_eq = [1-f f]';
M_ss = 1/D*[(1-f)*(A) f*(A)]';
term_1 = inv(I - E_s*E_m*E_r*C);
term_2 = (E_s*E_m*(I - E_r) + I-E_s)*M_eq;
term_3 = E_s*(I-E_m)*M_ss;
m_norm = term_1*(term_2 + term_3);

%output is the term to be minimized
output = m - M_z(:,1)./m_norm(1);
