/* % -*- mode: C; mode: fold -*-
Copyright (C) 2001 Leopoldo Cerbaro <redbliss@libero.it>
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 2 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, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Compute the Jacobi elliptic functions sn(u|m), cn(u|m) and dn(u|m) for
argument u (real or complex) and parameter m.
usage: [sn,cn,dn] = ellipj(u,m[,tol])
u and can be complex.
m is restricted to 0 <= m <= 1.
They can be scalars, matrix and scalar, scalar and matrix,
column and row, conformant matrices.
modified so u can be complex. Leopoldo Cerbaro redbliss@libero.it
Ref: Abramowitz, Milton and Stegun, Irene A
Handbook of Mathematical Functions, Dover, 1965
Chapter 16 (Sections 16.4, 16.13 and 16.15)
Based upon ellipj.m made by David Billinghurst <David.Billinghurst@riotinto.com>
and besselj.cc
Author: Leopoldo Cerbaro <redbliss@libero.it>
Created: 15 December 2001
*/
#include "oct.h"
#include "lo-ieee.h" /* for octave_NaN */
#ifdef USE_OCTAVE_NAN
#define lo_ieee_nan_value() octave_NaN
#endif
static void
gripe_ellipj_arg ( const char *arg)
{
error ("ellipj: expecting scalar or matrix as %s argument", arg);
}
const double eps = 2.220446049e-16;
const int Nmax = 16;
static void
sncndn ( double u, double m, double& sn, double& cn, double& dn, double& err) {
/* real */
double sqrt_eps, m1, t=0., si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi;
int n, Nn, ii;
if (m < 0. || m > 1.) {
warning ("ellipj: expecting 0. <= m <= 1."); /* -lc- */
sn = cn = dn = lo_ieee_nan_value ();
return;
}
sqrt_eps = sqrt(eps);
if (m < sqrt_eps) {
/* # For small m, ( Abramowitz and Stegun, Section 16.13 ) */
/*{{{*/
si_u = sin(u);
co_u = cos(u);
t = 0.25*m*(u-si_u*co_u);
sn = si_u - t * co_u;
cn = co_u + t * si_u;
dn = 1.0 - 0.5*m*si_u*si_u;
/*}}}*/
} else if ( (1.0 - m) < sqrt_eps ) {
/* For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 ) */
/*{{{*/
m1 = 1.0-m;
si_u = sinh(u);
co_u = cosh(u);
ta_u = tanh(u);
se_u = 1.0/co_u;
sn = ta_u + 0.25*m1*(si_u*co_u-u)*se_u*se_u;
cn = se_u - 0.25*m1*(si_u*co_u-u)*ta_u*se_u;
dn = se_u + 0.25*m1*(si_u*co_u+u)*ta_u*se_u;
/*}}}*/
} else {
/*{{{*/
/*
// Arithmetic-Geometric Mean (AGM) algorithm
// ( Abramowitz and Stegun, Section 16.4 )
*/
a[0] = 1.0;
b = sqrt(1.0-m);
c[0] = sqrt(m);
for (n = 1; n<Nmax; ++n) {
a[n] = (a[n-1]+b)/2;
c[n] = (a[n-1]-b)/2;
b = sqrt(a[n-1]*b);
if ( c[n]/a[n] < eps) break;
}
if ( n >= Nmax-1) {
// fprintf(stderr, "Not enough workspace\n");
err = 1.;
return;
}
Nn = n;
for ( ii = 1; n>0; ii = ii*2, --n) ; // pow(2, Nn)
phi = ii*a[Nn]*u;
for ( n = Nn; n > 0; --n) {
t = phi;
phi = (asin((c[n]/a[n])* sin(phi))+phi)/2.;
}
sn = sin(phi);
cn = cos(phi);
dn = cn/cos(t-phi);
/*}}}*/
}
return;
}
static void
sncndn ( Complex& u, double m,
Complex& sn, Complex& cn, Complex& dn, double& err) {
double m1 = 1.-m, ss1, cc1, dd1;
sncndn( imag(u), m1, ss1, cc1, dd1, err);
if ( real(u) == 0.) {
/* u is pure imag: Jacoby imag. transf. */
/*{{{*/
sn = Complex (0. , ss1/cc1);
cn = 1/cc1; // cn.imag = 0.;
dn = dd1/cc1; // dn.imag = 0.;
/*}}}*/
} else {
/* u is generic complex */
/*{{{*/
double ss, cc, dd, ddd;
sncndn( real(u), m, ss, cc, dd, err);
ddd = cc1*cc1 + m*ss*ss*ss1*ss1;
sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd);
cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd);
dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
/*}}}*/
}
return;
}
/*
## tests taken from test_ellipj.m
%!test
%! u1 = pi/3; m1 = 0;
%! res1 = [sin(pi/3), cos(pi/3), 1];
%! [sn,cn,dn]=ellipj(u1,m1);
%! assert([sn,cn,dn], res1, 10*eps);
%!test
%! u2 = log(2); m2 = 1;
%! res2 = [ 3/5, 4/5, 4/5 ];
%! [sn,cn,dn]=ellipj(u2,m2);
%! assert([sn,cn,dn], res2, 10*eps);
%!test
%! u3 = log(2)*1i; m3 = 0;
%! res3 = [3i/4,5/4,1];
%! [sn,cn,dn]=ellipj(u3,m3);
%! assert([sn,cn,dn], res3, 10*eps);
%!test
%! u4 = -1; m4 = tan(pi/8)^4;
%! res4 = [-0.8392965923,0.5436738271,0.9895776106];
%! [sn,cn,dn]=ellipj(u4, m4);
%! assert([sn,cn,dn], res4, 1e-10);
%!test
%! u5 = -0.2 + 0.4i; m5 = tan(pi/8)^4;
%! res5 = [ -0.2152524522 + 0.402598347i, ...
%! 1.059453907 + 0.08179712295i, ...
%! 1.001705496 + 0.00254669712i ];
%! [sn,cn,dn]=ellipj(u5,m5);
%! assert([sn,cn,dn], res5, 1e-9);
%!test
%! u6 = 0.2 + 0.6i; m6 = tan(pi/8)^4;
%! res6 = [ 0.2369100139 + 0.624633635i, ...
%! 1.16200643 - 0.1273503824i, ...
%! 1.004913944 - 0.004334880912i ];
%! [sn,cn,dn]=ellipj(u6,m6);
%! assert([sn,cn,dn], res6, 1e-8);
%!test
%! u7 = 0.8 + 0.8i; m7 = tan(pi/8)^4;
%! res7 = [0.9588386397 + 0.6107824358i, ...
%! 0.9245978896 - 0.6334016187i, ...
%! 0.9920785856 - 0.01737733806i ];
%! [sn,cn,dn]=ellipj(u7,m7);
%! assert([sn,cn,dn], res7, 1e-10);
%!test
%! u=[0,pi/6,pi/4,pi/2]; m=0;
%! res = [0,1/2,1/sqrt(2),1;1,cos(pi/6),1/sqrt(2),0;1,1,1,1];
%! [sn,cn,dn]=ellipj(u,m);
%! assert([sn;cn;dn],res, 100*eps);
%! [sn,cn,dn]=ellipj(u',0);
%! assert([sn,cn,dn],res', 100*eps);
## XXX FIXME XXX
## need to check [real,complex]x[scalar,rowvec,colvec,matrix]x[u,m]
*/
DEFUN_DLD (ellipj, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{sn}, @var{cn}, @var{dn}] =} \
ellipj (@var{u}, @var{m}, @var{err})\n\
Compute the Jacobi elliptic functions sn, cn, dn of complex argument and real parameter.\n\
\n\
If @var{m} is a scalar, the results are the same size as @var{u}.\n\
If @var{u} is a scalar, the results are the same size as @var{m}.\n\
If @var{u} is a column vector and @var{m} is a row vector, the\n\
results are matrices with @code{length (@var{u})} rows and\n\
@code{length (@var{m})} columns. Otherwise, @var{u} and\n\
@var{m} must conform and the results will be the same size.\n\
\n\
The value of @var{u} may be complex.\n\
The value of @var{m} must be 0 <= m <= 1. .\n\
\n\
If requested, @var{err} contains the following status information\n\
and is the same size as the result.\n\
\n\
@enumerate 0\n\
@item\n\
Normal return.\n\
@item\n\
Error---no computation, algorithm termination condition not met,\n\
return @code{NaN}.\n\
@end enumerate\n\
@end deftypefn")
{
octave_value_list retval;
int nargin = args.length ();
if (nargin == 2 ) {
octave_value u_arg = args(0);
octave_value m_arg = args(1);
if (m_arg.is_scalar_type ()) { // m is scalar
double m = args(1).double_value ();
if (! error_state) {
if (u_arg.is_scalar_type ()) { /* u scalar */
/*{{{*/
if (u_arg.is_real_type ()) { // u real
double u = args(0).double_value ();
if (! error_state) {
double sn, cn, dn;
double err=0;
octave_value result;
sncndn(u, m, sn, cn, dn, err);
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3)
retval(3) = err;
} else
gripe_ellipj_arg ( "first");
} else { // u complex
Complex u = u_arg.complex_value ();
if (! error_state) {
Complex sn, cn, dn;
double err;
octave_value result;
sncndn( u, m, sn, cn, dn, err);
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3) retval(3) = err;
} else
gripe_ellipj_arg ( "second");
}
/*}}}*/
} else { /* u is matrix ( m is scalar ) */
/*{{{*/
ComplexMatrix u = u_arg.complex_matrix_value ();
if (! error_state) {
octave_value result;
int nr = u.rows ();
int nc = u.cols ();
ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
Matrix err (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
sncndn (u(i,j), m, sn(i,j), cn(i,j), dn(i,j), err(i,j));
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3) retval(3) = err;
} else
gripe_ellipj_arg ( "first");
/*}}}*/
}
} else
gripe_ellipj_arg ( "second");
} else { // m is matrix
Matrix m = args(1).matrix_value ();
if (! error_state) {
int mr = m.rows ();
int mc = m.cols ();
if (u_arg.is_scalar_type ()) { /* u is scalar */
/*{{{*/
octave_value result;
int nr = m.rows ();
int nc = m.cols ();
Matrix err (nr, nc);
if (u_arg.is_real_type ()) {
double u = u_arg.double_value ();
Matrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
if (! error_state) {
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3) retval(3) = err;
} else
gripe_ellipj_arg ( "first");
} else {
Complex u = u_arg.complex_value ();
ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
if (! error_state) {
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3) retval(3) = err;
} else
gripe_ellipj_arg ( "first");
}
/*}}}*/
} else { // u is matrix (m is matrix)
/*{{{*/
if (u_arg.is_real_type ()) { // u real matrix
Matrix u = u_arg.matrix_value ();
if (! error_state) {
int ur = u.rows ();
int uc = u.cols ();
if (mr == 1 && uc == 1) { // u column, m row
RowVector rm = m.row ((octave_idx_type)0);
ColumnVector cu = u.column ((octave_idx_type)0);
Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
Matrix err(ur,mc);
// octave_value result;
for (int j = 0; j < mc; j++)
for (int i = 0; i < ur; i++)
sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3) retval(3) = err;
} else if (ur == mr && uc == mc) {
Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
Matrix err(ur,mc);
// octave_value result;
for (int j = 0; j < uc; j++)
for (int i = 0; i < ur; i++)
sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3) retval(3) = err;
} else
error("u m invalid");
} else
gripe_ellipj_arg ( "first ");
} else { // u complex matrix
ComplexMatrix u = u_arg.complex_matrix_value ();
if (! error_state) {
int ur = u.rows ();
int uc = u.cols ();
if (mr == 1 && uc == 1) {
RowVector rm = m.row ((octave_idx_type)0);
ComplexColumnVector cu = u.column ((octave_idx_type)0);
ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
Matrix err(ur,mc);
// octave_value result;
for (int j = 0; j < mc; j++)
for (int i = 0; i < ur; i++)
sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3) retval(3) = err;
} else if (ur == mr && uc == mc) {
ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
Matrix err(ur,mc);
// octave_value result;
for (int j = 0; j < uc; j++)
for (int i = 0; i < ur; i++)
sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
retval (0) = sn;
retval (1) = cn;
retval (2) = dn;
if (nargout > 3) retval(3) = err;
} else
error("u m invalid");
} else
gripe_ellipj_arg ( "second");
}
/*}}}*/
}
} else
gripe_ellipj_arg ( "second");
} // m matrix
} else // wrong n. of argin
print_usage ("ellipj");
return retval;
}
/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
syntax highlighted by Code2HTML, v. 0.9.1