function M=mueller(S1,S2,S3,S4)
% Muiller matrix elements depends of 4 angles:
% incident and scattered zenith and azimuth angle
% M=[S11,S12,S13,S14
%    S21,S22,S23,S24
%    S31,S32,S33,S34
%    S41,S42,S43,S44]

M(1,1)=0.5*((abs(S1)).^2+(abs(S2)).^2+(abs(S3)).^2+(abs(S4)).^2);
M(1,2)=0.5*((abs(S2)).^2-(abs(S1)).^2+(abs(S4)).^2-(abs(S3)).^2);
M(1,3)=real(S2.*conj(S3)-S1.*conj(S4));
M(1,4)=imag(S2.*conj(S3)-S1.*conj(S4));
M(2,1)=0.5*((abs(S2)).^2-(abs(S1)).^2+(abs(S3)).^2-(abs(S4)).^2);
M(2,2)=0.5*((abs(S1)).^2+(abs(S2)).^2-(abs(S3)).^2-(abs(S4)).^2);
M(2,3)=real(S2.*conj(S3)-S1.*conj(S4));
M(2,4)=imag(S2.*conj(S3)+S1.*conj(S4));
M(3,1)=real(S2.*conj(S4)+S1.*conj(S3));
M(3,2)=real(S2.*conj(S4)-S1.*conj(S3));
M(3,3)=real(S1.*conj(S2)+S3.*conj(S4));
M(3,4)=imag(S2.*conj(S1)+S4.*conj(S3));
M(4,1)=imag(S4.*conj(S2)+S1.*conj(S3));
M(4,2)=imag(S4.*conj(S2)-S1.*conj(S3));
M(4,3)=imag(S1.*conj(S2)-S3.*conj(S4));
M(4,4)=real(S1.*conj(S2)-S3.*conj(S4));
Pol=sqrt(M(2,1).^2+M(3,1).^2)./M(1,1);
% elements of this matrix are in case of:
% 1. spherical particles S11=S22 and S33=S44
% 2. Rayleigh scattering S12=0

