## -*- texinfo -*-
## @deftypefn {Function File} {} placet_get_sigma_matrix (@var{B})
## Return the sigma matrix of beam `B', where `B' is a 6-columns matrix in guinea-pig format, or a 17-columns matrix in sliced beam format
## @end deftypefn

function S=placet_get_sigma_matrix(B)
  if nargin==1 && ismatrix(B) && columns(B)==6
    Cx=covm(B(:,[2,5]));
    Cy=covm(B(:,[3,6]));
    S=zeros(4,4);
    S(1:2,1:2)=Cx;
    S(3:4,3:4)=Cy;
  elseif nargin==1 && ismatrix(B) && columns(B)==17
    sigma_w = sum(B(:,2));
    mu_x = B(:,2)'*B(:,4) / sigma_w;
    mu_xp = B(:,2)'*B(:,5) / sigma_w;
    mu_y = B(:,2)'*B(:,6) / sigma_w;
    mu_yp = B(:,2)'*B(:,7) / sigma_w;
    sigma_yy = B(:,2)'*(B(:,11)+((B(:,6)-mu_y).*(B(:,6)-mu_y)));
    sigma_yyp = B(:,2)'*(B(:,12)+((B(:,6)-mu_y).*(B(:,7)-mu_yp)));
    sigma_ypyp = B(:,2)'*(B(:,13)+((B(:,7)-mu_yp).*(B(:,7)-mu_yp)));
    sigma_xx = B(:,2)'*(B(:,8)+((B(:,4)-mu_x).*(B(:,4)-mu_x)));
    sigma_xxp = B(:,2)'*(B(:,9)+((B(:,4)-mu_x).*(B(:,5)-mu_xp)));
    sigma_xpxp = B(:,2)'*(B(:,10)+((B(:,5)-mu_xp).*(B(:,5)-mu_xp)));
    S(1:2,1:2) = [sigma_xx  sigma_xxp; sigma_xxp sigma_xpxp ]/ sigma_w;
    S(3:4,3:4) = [sigma_yy  sigma_yyp; sigma_yyp sigma_ypyp ]/ sigma_w;
  else
    help placet_get_sigma_matrix
  endif
endfunction
 
