Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / linear-algebra-2.2.0 / thfm.m

diff --git a/octave_packages/linear-algebra-2.2.0/thfm.m b/octave_packages/linear-algebra-2.2.0/thfm.m
new file mode 100644 (file)
index 0000000..8702d0f
--- /dev/null
+++ b/octave_packages/linear-algebra-2.2.0/thfm.m
@@ -0,0 +1,98 @@
+## Copyright (C) 2001 Rolf Fabian <fabian@tu-cottbus.de>
+## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
+## Copyright (C) 2011 Philip Nienhuis <pr.nienhuis@hccnet.nl>
+## Copyright (C) 2011 Carnë Draug <carandraug+dev@gmail.com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {@var{y} =} thfm (@var{x}, @var{mode})
+## Trigonometric/hyperbolic functions of square matrix @var{x}.
+##
+## @var{mode} must be the name of a function. Valid functions are 'sin', 'cos',
+## 'tan', 'sec', 'csc', 'cot' and all their inverses and/or hyperbolic variants,
+## and 'sqrt', 'log' and 'exp'.
+##
+## The code @code{thfm (x, 'cos')} calculates matrix cosinus @emph{even if} input
+## matrix @var{x} is @emph{not} diagonalizable.
+##
+## @emph{Important note}:
+## This algorithm does @emph{not} use an eigensystem similarity transformation. It
+## maps the @var{mode} functions to functions of @code{expm}, @code{logm} and
+## @code{sqrtm}, which are known to be robust with respect to non-diagonalizable
+## ('defective') @var{x}.
+##
+## @seealso{funm}
+## @end deftypefn
+
+function y = thfm (x,M)
+  ## minimal arg check only
+  if ( nargin != 2 || !ischar (M) || ischar (x) )
+    print_usage;
+  endif
+
+  ## look for known functions of sqrt, log, exp
+  I = eye (size (x));
+
+  switch (M)
+    case {'cos'}
+      if (isreal(x))  y = real( expm( i*x ) );
+      else            y = ( expm( i*x ) + expm( -i*x ) ) / 2;
+      endif
+    case {'sin'}
+      if (isreal(x))  y = imag( expm( i*x ) );
+      else            y = ( expm( i*x ) - expm( -i*x ) ) / (2*i);
+      endif
+    case {'tan'}
+      if (isreal(x))  t = expm( i*x );    y = imag(t)/real(t);
+      else            t = expm( -2*i*x ); y = -i*(I-t)/(I+t);
+      endif
+    case {'cot'}
+      if (isreal(x))  t = expm( i*x );    y = real(t)/imag(t);
+      else            t = expm( -2*i*x ); y = i*(I+t)/(I-t);
+      endif
+    case {'sec'}
+      if (isreal(x))  y = inv( real(expm(i*x)) );
+      else            y = inv( expm(i*x)+expm(-i*x) )*2 ;
+      endif
+    case {'csc'}
+      if (isreal(x))  y = inv( imag(expm(i*x)) );
+      else            y = inv( expm(i*x)-expm(-i*x) )*2i;
+      endif
+    case {'log'}      y = logm(x);
+    case {'exp'}      y = expm(x);
+    case {'cosh'}     y = ( expm(x)+expm(-x) )/2;
+    case {'sinh'}     y = ( expm(x)-expm(-x) )/2;
+    case {'tanh'}     t = expm( -2*x ); y = (I - t)/(I + t);
+    case {'coth'}     t = expm( -2*x ); y = (I + t)/(I - t);
+    case {'sech'}     y = 2 * inv( expm(x) + expm(-x) );
+    case {'csch'}     y = 2 * inv( expm(x) - expm(-x) );
+    case {'asin'}     y = -i * logm( i*x + sqrtm(I - x*x) );
+    case {'acos'}     y =  i * logm( x - i*sqrtm(I - x*x) );
+    case {'atan'}     y = -i/2 * logm( (I + i*x)/(I - i*x) );
+    case {'acot'}     y =  i/2 * logm( (I + i*x)/(i*x - I) );
+    case {'asec'}     y = i * logm( ( I - sqrtm(I - x*x) ) / x );
+    case {'acsc'}     y = -i * logm( i*( I + sqrtm(I - x*x) ) / x );
+    case {'sqrt'}     y = sqrtm(x);
+    case {'asinh'}    y = logm( x + sqrtm (x*x + I) );
+    case {'acosh'}    y = logm( x + sqrtm (x*x - I) );
+    case {'atanh'}    y = logm( (I + x)/(I - x) ) / 2;
+    case {'acoth'}    y = logm( (I + x)/(x - I) ) / 2;
+    case {'asech'}    y = logm( (I + sqrtm (I - x*x)) / x );
+    case {'acsch'}    y = logm( (I + sqrtm (I + x*x)) / x );
+    otherwise
+      error ("thfm doesn't support function %s - try to use funm instead.", M);
+  endswitch
+
+endfunction