//  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 [xmean] = CL_mean(x, meth, rc)
// Mean value based on quadrature
//
// Calling Sequence
// xmean = CL_mean(x, meth [, rc])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the mean value by quadrature. </p>
// <p>Assuming that the implicit (evenly spaced) abscissas are t1 .. tn, the mean value 
// is defined by: </p>
// <p>xmean = Integral from t1 to tn of x(t), divided by tn-t1 (integral approximated by quadrature). </p>
// <p></p> 
// <p>The mean value can be computed on the rows (rc = "r") or on the columns (rc = "c"), 
// exactly as the function "mean" does. If rc = "r", the result is a row vector. 
// If rc = "c", the result is a column vector. </p>
// <p></p>
// <p>Three methods are available: </p>
// <p>- trap: integral evaluated by trapezoidal method</p>
// <p>- simp: integral evaluated by Simpson's rule</p>
// <p>- boole: integral evaluated by Boole's rule</p>
// <p></p>
// <p>The required number of values depends on the method used: </p>
// <p>- trap: any number</p>
// <p>- simp: odd number (e.g. 1, 3, 5, 7...) </p>
// <p>- boole: multiple of 4 plus 1 (e.g. 1, 5, 9, 13...) </p>
// <p></p>
// <p>Notes: </p>
// <p>- rc can be omitted for a row vector or a column vector. The result is then a real number. </p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// x: Matrix of real values (PxN)
// meth: (string) Method used: "trap", "simp", "boole" (1x1)
// rc: (string, optional) Direction: "r": mean computed on rows or "c": mean  computed on columns. For default, see above.  
// xmean: Mean value (Px1) or (1xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// t = linspace(0, %pi, 101); 
// CL_mean(sin(t), "trap")  - 2/%pi
// CL_mean(sin(t), "simp")  - 2/%pi
// CL_mean(sin(t), "boole") - 2/%pi

// Declarations:

// internal functions

// Trapezoidal method
// NB: valid if number of points n = k+1 (k >= 1)
// 2 points: mean = (x1+x2)/2
// n points: k*mean = (x1+x2)/2 + (x2+x3)/2 + ... 

function [xmean] = meanTrap(x, rc)
  n = size(x, rc) - 1; // number of intervals
  if (rc == "r")
    xmean = (2*sum(x,"r") - x(1,:) - x($,:)) / (2*n); 
  else
    xmean = (2*sum(x,"c") - x(:,1) - x(:,$)) / (2*n); 
  end
endfunction


// Simpson's rule
// NB: valid if number of points n = 2*k+1 (k >= 1)
// 3 points: mean = (x1+4*x2+x3)/6
// n points: k*mean = (x1+4*x2+x3)/6 + (x3+4*x4+x5)/6 + ...

function [xmean] = meanSimpson(x, rc)
  n = size(x, rc) - 1; // number of intervals
  if (rc == "r")
    xmean = (2*sum(x(1:2:$,:),"r") + 4*sum(x(2:2:$,:),"r") - x(1,:) - x($,:)) / (3*n); 
  else
    xmean = (2*sum(x(:,1:2:$),"c") + 4*sum(x(:,2:2:$),"c") - x(:,1) - x(:,$)) / (3*n); 
  end
endfunction


// Boole's rule - formulas:
// NB: valid if number of points n = 4*k+1 (k>=1)
// 5 points: mean = (7*x1+32*x2+12*x3+32*x4+7*x5)/90
// n points: k*mean = (7*x1+32*x2+12*x3+32*x4+7*x5)/90 + (7*x5+32*x6+12*x7+32*x8+7*x9)/90 + ...

function [xmean] = meanBoole(x, rc)
  n = size(x, rc) - 1; // number of intervals
  if (rc == "r")
    xmean = (28*sum(x(1:4:$,:),"r") + 64*sum(x(2:2:$,:),"r") + ..
             24*sum(x(3:4:$,:),"r") - 14*x(1,:) - 14*x($,:)) / (45*n); 
  else
    xmean = (28*sum(x(:,1:4:$),"c") + 64*sum(x(:,2:2:$),"c") + ..
             24*sum(x(:,3:4:$),"c") - 14*x(:,1) - 14*x(:,$)) / (45*n); 
  end
endfunction


// Code:
if (~exists("rc", "local")); rc = []; end

if (rc <> "r" & rc <> "c" & rc <> [])
  CL__error("Invalid value for argument ''rc''"); 
end

if (meth <> "trap" & meth <> "simp" & meth <> "boole")
  CL__error("Invalid method"); 
end
 
if (x == [])
  xmean = []; 
  return; 
end

// if rc == [] 
// => replace by "r" if column vector or "c" if row vector

if (rc == [])
  n = size(x); 
  if (n(1) == 1); 
    // column vector
    rc = "c"; 
  elseif (n(2) == 1)
    // row vector
    rc = "r"; 
  else
    CL__error("Argument ''rc'' missing");
  end     
end

// size in the "rc" direction
nrc = size(x, rc); 

// check number of rows/columns is odd
if ((meth == "simp" & modulo(nrc,2) <> 1) | (meth == "boole" & modulo(nrc,4) <> 1))
  CL__error("Invalid size for argument ''x''");
end

if (nrc == 1)
  xmean = x; 

elseif (meth == "trap")
  xmean = meanTrap(x, rc); 
  
elseif (meth == "simp")
  xmean = meanSimpson(x, rc); 
  
else
  xmean = meanBoole(x, rc); 
end


endfunction