//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [covm] = CL_cor2cov(corm, sd)
// Correlation to covariance matrix
//
// Calling Sequence
// covm = CL_cor2cov(corm, sd)
//
// Description
// <itemizedlist><listitem>
// <p> Computes the covariance matrix from the correlation matrix and the standard deviation vector. </p>
// <p> The covariance matrix is defined as follows: </p>
// <p> - <b>covm</b>(i,i) = <b>sd</b>(i)^2 </p>
// <p> - <b>covm</b>(i,j) = <b>corm</b>(i,j) * <b>sd</b>(i) * <b>sd</b>(j) </p>
// <p></p>
// <p><b>Notes</b>:</p>
// <p>- The correlation matrix should be symetrical (not checked). </p>
// <p>- An error is raised if one diagonal term is not 1. </p>
// <p>- An error is raised if one standard deviation value is not positive. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// corm: Correlation matrix (NxNxK)
// sd: Standard deviation vector (NxK)
// covm: Covariance matrix (NxNxK)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_cov2cor
//
// Examples
// corm = [1, 0, 0.5; 0, 1, -0.8; 0.5, -0.8, 1];  
// sd = [1; 2; 3]; 
// covm = CL_cor2cov(corm, sd);
//
// // Consistency check: 
// [corm2, sd2] = CL_cov2cor(covm)


// Declarations:

// correlation to covariance matrix (2D case)    
function [covm] = cor2cov(corm, sd)

  I = find(sd < 0);
  if (I <> []); CL__error("Invalid term in standard deviation vector"); end
  
  I = find(abs(diag(corm) - 1) > 1 + 2*%eps);
  if (I <> []); CL__error("Invalid term in correlation matrix"); end
  
  w = diag(sd);
  covm = w * corm * w;    
endfunction
  
  
// Code:

// special case: 
// output = [] if input = []
if (corm == [] & sd == [])
  covm = []; 
  return; 
end

// add 1 => also works for matrices
s1 = [size(corm), 1];
s2 = [size(sd), 1];
  
if (s1(1) <>  s1(2) | s1(3) <> s2(2) | s2(3) <> 1)
  CL__error("Invalid arguments sizes");
end
  
// Matrix :
if (s1(3) == 1)
  
  // Matrix    
  covm = cor2cov(corm, sd);
   
else
  // Hypermatrix :
  covm = zeros(corm);
    
  for (i = 1 : s1(3))   
    covmi = cor2cov(corm(:,:,i), sd(:,i));
    covm(:,:,i) = covmi;
  end
end

endfunction