[x,obj_hat]=simdata_T2star(obj,we,N1,N2,M1,M2,d1,d2,delta1,delta2,T2star);

% Description: simulates MRI data with T2* term included according to the forward model 
% Output: simulated data 'x'

% Hien M. Nguyen, University of Illinois at Urbana-Champaign
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,obj_hat]=simdata_T2star(obj,we,N1,N2,M1,M2,d1,d2,delta1,delta2,T2star);

n1=[0:N1-1];
n2=[0:N2-1].';

i = sqrt(-1);
j = sqrt(-1);

% k and r traj assume starts from origin, so we need modulation for symmetric k-space
L1 = (N1)/2;
L2 = (N2)/2;

for k=0:N2-1
    for l=0:N1-1
        obj_hat(k+1,l+1)=obj(k+1,l+1).*exp(j*2*pi*L1*l*delta1/N1).*exp(j*2*pi*L2*k*delta2/N2);
    end;
end;

for m1=0:M1-1   
    A(:,:,m1+1)=exp(-j*m1.*(we(n2+1,n1+1)*d1+2*pi*delta1*repmat(n1,N2,1)/N1)).*exp(-m1*d1./T2star(n2+1,n1+1));  % additional T2star term
end;

for m2=0:M2-1
    B(:,:,m2+1)=exp(-j*m2.*(we(n2+1,n1+1)*d2+2*pi*repmat(n2,1,N1)*delta2/N2)).*exp(-m2*d2./T2star(n2+1,n1+1));  % additional T2star term
end;

for m1=0:M1-1
    for m2=0:M2-1
        x(m2+1,m1+1)=sum(sum(obj_hat(:,:).*A(:,:,m1+1).*B(:,:,m2+1),1),2);;
    end;
end;