function [newA, newrhs, nblock]=rbmatreform(Amat, ymeas, ymodel, form)
%
% [newA, newrhs, nblock]=rbmatreform(Amat, ymeas, ymodel, output)
%
% Reformat the matrix equation A*x=(ymeas-ymodel) to choose between the log-amplitude/phase
% form or the complex form
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     Amat: the LHS of the matrix equation
%     ymeas: the vector that stores the measured data for the model to fit 
%     ymodel: the model predicted measurements at all source detector pairs, 
%           with a length matching that of ymeas
%     form: a string to indicate which form of output system, can be one of
%         'complex': no transformation, Amat and RHS keep original forms,
%             if inputs are real(complex), output are real(or complex)
%         'real': newA and newrhs are real matrices, x is assumed to be
%             real
%         'reim': newA and newrhs are real matrices, x is assumed to be
%             complex and is expanded to [real(x), imag(x)]' and newA is
%             expanded accordingly
%         'logphase': newA and newrhs are real matrices, x is assumed to be
%             complex and is expanded to [real(x) imag(x)]'; newrhs is
%             expanded as [log(ymeas-ymodel), angle(ymeans)-angle(ymodel)]'
%             and newA is transformed accordingly, please see Dr. Fang's
%             PhD Thesis, Section 3.2, Eq. 3.25
%
% output:
%     newA: the reformed LHS matrix 
%     newrhs: the reformed RHS matrix
%     nblock: nblock=length(newrhs)/size(Amat,1)
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details 
%
% -- this function is part of Redbird-m toolbox
%

if(nargin<4)
    form='complex';
end
nblock=1;

if isstruct(Amat)
    Amat = cell2mat(squeeze(struct2cell(Amat)));
    ymeas = struct2cell(ymeas);
    ymodel = struct2cell(ymodel);
    if (size(ymeas,1) < size(ymeas,2))
        ymeas = ymeas.';
        ymodel = ymodel.';
    end
    ymeas = cell2mat(ymeas);
    ymodel = cell2mat(ymodel);
end

rhs=ymeas-ymodel;

if(strcmp(form,'real') || strcmp(form,'reim'))
    newA=real(Amat);
    newrhs=real(rhs);

    if(~isreal(rhs) && ~isreal(Amat))
        if(strcmp(form,'reim'))
           newA=[real(Amat) -imag(Amat); 
                 imag(Amat) real(Amat)];
        else
           newA=[newA; imag(Amat)];
        end
	    newrhs=[newrhs; imag(rhs)];
        nblock=1;
    end
    return;
end

if(strcmp(form,'logphase'))
    temp=repmat(conj(ymodel)./abs(ymodel.*ymodel),1,size(Amat,2)).*Amat;
    if(isreal(ymodel))
        newA=real(temp);
        newrhs=log(abs(ymeas))-log(abs(ymodel));
    else
        newA=[real(temp) ; imag(temp)];
	    newrhs=[log(abs(ymeas)) - log(abs(ymodel)); 
                angle(ymeas) - angle(ymodel)];
        nblock=1;
    end
    return;
end


%{
if (~isstruct(Amat))
    Amat = struct('J',Amat);
    ymeas = struct('detphi',ymeas);
    ymodel = struct('detphi',ymodel);
end

rfcw = length(Amat);
Aout = [];
rhsout = [];

for ii = rfcw:-1:1
    y0 = ymeas(ii).detphi;
    y1 = ymodel(ii).detphi;
    A = Amat(ii).J;
    
    rhs=y0-y1;

    if(strcmp(form,'complex'))
        newA=A;
        newrhs=rhs;
    end

    if(strcmp(form,'real') || strcmp(form,'reim'))
        newA=real(A);
        newrhs=real(rhs);

        if(~isreal(rhs) && ~isreal(A))
            if(strcmp(form,'reim'))
               newA=[real(A) -imag(A); 
                     imag(A) real(A)];
            else
               newA=[newA; imag(A)];
            end
            newrhs=[newrhs; imag(rhs)];
            nblock=1;
        end
    end

    if(strcmp(form,'logphase'))
        temp=repmat(conj(y1)./abs(y1.*y1),1,size(A,2)).*A;
        if(isreal(y1))
            newA=real(temp);
            newrhs=log(abs(y0))-log(abs(y1));
        else
            newA=[real(temp) ; imag(temp)];
            newrhs=[log(abs(y0))-log(abs(y1)); angle(y0)-angle(y1)];
            nblock=1;
        end
    end
    
    Aout = [Aout; newA];
    rhsout = [rhsout; newrhs];
end

newA = Aout;
newrhs = rhsout;
%}
