function yamachuchi_matrix = yamaguchi(C)
%对于给定的C矩阵做yamaguchi分解，返回结果是个各个散射的能量是一个row X col X 4的矩阵
%P_s（单词散射）,P_d（二面角）,P_v（体散射）,P_c（旋状体）
row=size(C.C11,1);
col=size(C.C11,2);
i=sqrt(-1);
yamachuchi_matrix = zeros(row,col,4);
for p=1:col
    for q=1:row
        TP = C.C11(q,p)+C.C22(q,p)+C.C33(q,p); 
        P_c = 2*abs(imag((C.C12_real(q,p)+i*C.C12_imag(q,p)-conj(C.C23_real(q,p)+i*C.C23_imag(q,p)))/sqrt(2)));
        a = 10*log10(C.C33(q,p)/C.C11(q,p));
        
        %if a<-2
           % P_v = (15*C.C22(q,p))/4 - (15*P_c)/8;
           % tS = (C.C11(q,p)+C.C33(q,p)+2*C.C13_real(q,p)-P_v)/2;
           %tD = (C.C11(q,p)+C.C33(q,p)-2*C.C13_real(q,p))/2-(7*C.C22(q,p))/8-P_c/16;
           %tC = (C.C11(q,p)-C.C33(q,p)-(C.C13_real(q,p)+i*C.C13_imag(q,p))+(C.C13_real(q,p)-i*C.C13_imag(q,p)))/2-P_v/6;
        %elseif a<2
            P_v = 4*C.C22(q,p) - 2*P_c;
            tS = (C.C11(q,p)+C.C33(q,p)+2*C.C13_real(q,p))/2-2*C.C22(q,p)+P_c;
            tD = (C.C11(q,p)+C.C33(q,p)-2*C.C13_real(q,p))/2-C.C22(q,p);
            tC = (C.C11(q,p)-C.C33(q,p)-(C.C13_real(q,p)+i*C.C13_imag(q,p))+(C.C13_real(q,p)-i*C.C13_imag(q,p)))/2;
        %elseif a>2
            %P_v = (15*C.C22(q,p))/4 - (15*P_c)/8;
            %tS = (C.C11(q,p)+C.C33(q,p)+2*C.C13_real(q,p))/2-2*C.C22(q,p)+P_c;
            %tD = (C.C11(q,p)+C.C33(q,p)-2*C.C13_real(q,p))/2-C.C22(q,p);
            %tC = (C.C11(q,p)-C.C33(q,p)-(C.C13_real(q,p)+i*C.C13_imag(q,p))+(C.C13_real(q,p)-i*C.C13_imag(q,p)))/2;
        %end
        
        if (P_v+P_c)<TP
            C_0 = C.C13_real(q,p)+i*C.C13_imag(q,p)-C.C22(q,p)/2-P_c/2;
            if real(C_0)<0
                P_s = tS-((abs(tC))^2)/tD;
                P_d = tD+((abs(tC))^2)/tD;
            else
                P_s = tS+((abs(tC))^2)/tS;
                P_d = tD-((abs(tC))^2)/tS;
            end
            
            if (P_s>0)&&(P_d<0)
                P_d = 0;
                P_s = TP-P_v-P_c;
            elseif (P_s<0)&&(P_d>0)
                P_s = 0;
                P_d = TP-P_v-P_c;
            end
            
        else
            P_s = 0;
            P_d = 0;
            P_v = TP-P_c;
        end
        
        
        yamachuchi_matrix(q,p,1) = abs(P_s);
        yamachuchi_matrix(q,p,2) = abs(P_d);
        yamachuchi_matrix(q,p,3) = abs(P_v);
        yamachuchi_matrix(q,p,4) = abs(P_c);
        
    end
end