function [sAImax,d_max,G_max]=DemodGouy(mSens,probeIndexes,cDrvNames,isitDC)
% see LIGO-T080186 Appendix B
% find demod phases and Gouy phases that maximize signals / minimize signal contamination
% isitDC: if they are DC signals, put something here. Default is off(regard as RF signals).

DCflag = 1;
try
    isitDC; 
catch
    DCflag = 0;
end

sAI=mSens(probeIndexes(1),:);
sAQ=mSens(probeIndexes(2),:);
sBI=mSens(probeIndexes(3),:);
sBQ=mSens(probeIndexes(4),:);
Ndrv=length(cDrvNames);

alpha=sAI.^2-sAQ.^2+sBI.^2-sBQ.^2;
beta=sAI.^2+sAQ.^2-sBI.^2-sBQ.^2;
delta=2*(sAI.*sAQ+sBI.*sBQ);
gamma=delta.^2+alpha.^2;
eta=2*(sAI.*sBI+sAQ.*sBQ);
d_max=atan((-alpha+sqrt(gamma))./delta);
G_max=atan((-beta+sqrt(gamma))./eta);
sAImax=(cos(d_max).*cos(G_max)).*sAI+(sin(d_max).*cos(G_max)).*sAQ+(cos(d_max).*sin(G_max)).*sBI+(sin(d_max).*sin(G_max)).*sBQ;

initcontam=100;   % initial contamination level
if DCflag
    Ngd=1;
    tableelement=ones(6,Ndrv)*initcontam;
    tableelement(1,:)=G_max/pi*180;  % Gouy phase that maximize signals
    tableelement(2,:)=sAImax; % maximum signal amplitude
    ylabels={'max Gouy','max ampl','Gouy','ampl','contami','N big sig'};
    for Gouyphase=(-90:90)/180*pi
        sAIprime=cos(Gouyphase).*sAI+sin(Gouyphase).*sBI;
        contam=1./abs(sAIprime)*sum(abs(sAIprime));
        for kk=1:Ndrv
            if contam(kk) < tableelement(5,kk)
                tableelement(3,kk)=Gouyphase/pi*180; % Gouy phase that minimizes contamination level
                tableelement(4,kk)=sAIprime(kk);     % signal amplitude for this DOF under Gouy phase above
                tableelement(5,kk)=contam(kk);       % contamination level for this DOF
                tableelement(6,kk)=sum(abs(sAIprime)/abs(sAIprime(kk))>=1);   % number of signals that is over than the signal from this DOF
            end
        end
    end
    for kk=1:Ndrv
        if tableelement(5,kk)==initcontam  % if contamination level is over than 'initcontam', returns NaN (because signal contamination is awful)
            tableelement(3:6,kk)=ones(4,1)*NaN;
        end
    end
else
    Ngd=2;
    tableelement=ones(8,Ndrv)*initcontam;
    tableelement(1,:)=G_max/pi*180;
    tableelement(2,:)=d_max/pi*180;
    tableelement(3,:)=sAImax;
    ylabels={'max Gouy','max demod','max ampl','Gouy','demod','ampl','contami','N big sig'};
    for Gouyphase=(-90:90)/180*pi
        for demodphase=(-90:90)/180*pi
            sAIprime=(cos(demodphase).*cos(Gouyphase)).*sAI+(sin(demodphase).*cos(Gouyphase)).*sAQ+(cos(demodphase).*sin(Gouyphase)).*sBI+(sin(demodphase).*sin(Gouyphase)).*sBQ;
            contam=1./abs(sAIprime)*sum(abs(sAIprime));
            for kk=1:Ndrv
                if contam(kk) < tableelement(7,kk)
                    tableelement(4,kk)=Gouyphase/pi*180;
                    tableelement(5,kk)=demodphase/pi*180;
                    tableelement(6,kk)=sAIprime(kk);
                    tableelement(7,kk)=contam(kk);
                    tableelement(8,kk)=sum(abs(sAIprime)/abs(sAIprime(kk))>=1);
                end
            end
        end
    end
    for kk=1:Ndrv
        if tableelement(7,kk)==initcontam
            tableelement(4:8,kk)=ones(5,1)*NaN;
        end
    end    
end
N=length(ylabels);

% Plot table
for ii=0:N
    line([0,Ndrv],[ii,ii],'Color','k','LineWidth',0.1);
end
for ii=0:Ndrv
    line([ii,ii],[0,N],'Color','k','LineWidth',0.1);
end
set(gca, 'XLim',[0,Ndrv],'YLim',[0,N]);
set(gca,'XTick',[0.5:Ndrv-0.5], 'YTick', [0.5:N-0.5]);
set(gca,'FontSize',10)
set(gca,'XTickLabel',cDrvNames,'YTickLabel',ylabels);
for ii=1:Ngd
    for jj=1:Ndrv
        text(jj-0.5,ii-0.5,num2str(tableelement(ii,jj), '%.2f'), 'HorizontalAlignment','center','FontSize',10);
    end
end
Nmaxamp=Ngd+1;
[m,maxjj]=max(abs(tableelement(Nmaxamp,:)));
if isnan(m)~=1
    for jj=1:Ndrv
        text(jj-0.5,Nmaxamp-0.5,exphyouji(num2str(tableelement(Nmaxamp,jj), '%1.1e')), 'HorizontalAlignment','center','FontSize',10,'color',[abs(tableelement(Nmaxamp,jj))/m,0,0]);
    end
end
for ii=Ngd+2:2*Ngd+1
    for jj=1:Ndrv
        text(jj-0.5,ii-0.5,num2str(tableelement(ii,jj), '%.0f'), 'HorizontalAlignment','center','FontSize',10);
    end
end
for ii=2*(Ngd+1):2*(Ngd+1)
    for jj=1:Ndrv
        text(jj-0.5,ii-0.5,exphyouji(num2str(tableelement(ii,jj), '%1.1e')), 'HorizontalAlignment','center','FontSize',10);
    end
end
for ii=2*(Ngd+1)+1:2*(Ngd+1)+1
    for jj=1:Ndrv
        text(jj-0.5,ii-0.5,num2str(tableelement(ii,jj), '%.2f'), 'HorizontalAlignment','center','FontSize',10);
    end
end
for ii=N:N
    for jj=1:Ndrv
        text(jj-0.5,ii-0.5,num2str(tableelement(ii,jj), '%.0f'), 'HorizontalAlignment','center','FontSize',10);
    end
end

set(gca,'XAxisLocation','top'); 
set(gca,'YDir','reverse');

function s=exphyouji(a)
s=strrep(strrep(strrep(a,'e+00','E'),'e-00','E-'),'e-0','E-');