[z_found,a_found, R2star, z_noswap, ind]=HR0_v2_2_T2star_notconst_sim2(t,N,M,d2,lambda,betta)

% Description: performs harmonic retrieval
% Input:    - t: data 
%           - M: # of data samples
%           - N: # of harmonics
%           - d2: sampling time interval along phase-encodings (need to find R2star)
%           - lambda & betta: two regularization parameters for tunning. Lambda can be estimated if know rough estimate of T2*,
%             betta is the weighting between data and regularization terms
% Output:   - {z_found,a_found}: harmonics, from which will estimate image and field map
%           - R2star 

% For further information, see Nguyen et al., Joint estimation and
% correction of geometric distortions for EPI functional MRI

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

function [z_found,a_found, R2star, z_noswap, ind]=HR0_v2_2_T2star_notconst_sim2(t,N,M,d2,lambda,betta)
% like HR0_v2_2 but doesn't enforce magnitude of recovered {z_n} to be 1

% -------------Compute matrix T -------------------------
 c=t(1:(M+1-N),1).';
 r=t((M+1-N):M,1).';
 T=hankel(c,r);

        % ----------------- Recover zi and ai using HR --------------------
   
        [z_found]=LSProny_T2star(T,t,N,M,lambda,betta);
        
        [phi,order]=sort(angle(z_found));                      
        z_found=z_found(order); 
        
        % reverse order of z since original are in decreasing order
        z_found=z_found(end:-1:1);
  
        [var,ind]=min(abs(angle(z_found)));
        z_noswap=z_found;
        
        % ind will be first point of first half
        z_found=[z_found(ind:end);z_found(1:ind-1)];
      
        % correcting for angle(z)<-3.1415
        index=find(angle(z_found)<=-3.14);
        z_found(index)=abs(z_found(index))*exp(j*pi);
        
        % Find a's as a least squares solution ----------------------------
        A=Matrix_Z(z_found.',N,M);
        a_found=pinv(A)*t;
        
        % Find R2star map
        R2star=-log(abs(z_found))./d2;