# include # include # include # include "sncndn.h" /******************************************************************************/ double jacobi_cn ( double u, double m ) /******************************************************************************/ /* Purpose: jacobi_cn() evaluates the Jacobi elliptic function CN(U,M). Licensing: This code is distributed under the MIT license. Modified: 26 June 2025 Author: Original ALGOL version by Roland Bulirsch. This version by John Burkardt Reference: Roland Bulirsch, Numerical calculation of elliptic integrals and elliptic functions, Numerische Mathematik, Volume 7, Number 1, 1965, pages 78-90. Input: double U, M, the arguments. Output: double JACOBI_CN, the function value. */ { double cn; double dn; double sn; sncndn ( u, m, &sn, &cn, &dn ); return cn; } /******************************************************************************/ double jacobi_dn ( double u, double m ) /******************************************************************************/ /* Purpose: jacobi_dn() evaluates the Jacobi elliptic function DN(U,M). Licensing: This code is distributed under the MIT license. Modified: 26 June 2025 Author: Original ALGOL version by Roland Bulirsch. This version by John Burkardt Reference: Roland Bulirsch, Numerical calculation of elliptic integrals and elliptic functions, Numerische Mathematik, Volume 7, Number 1, 1965, pages 78-90. Input: double U, M, the arguments. Output: double JACOBI_DN, the function value. */ { double cn; double dn; double sn; sncndn ( u, m, &sn, &cn, &dn ); return dn; } /******************************************************************************/ double jacobi_sn ( double u, double m ) /******************************************************************************/ /* Purpose: jacobi_sn() evaluates the Jacobi elliptic function SN(U,M). Licensing: This code is distributed under the MIT license. Modified: 26 June 2025 Author: Original ALGOL version by Roland Bulirsch. This version by John Burkardt Reference: Roland Bulirsch, Numerical calculation of elliptic integrals and elliptic functions, Numerische Mathematik, Volume 7, Number 1, 1965, pages 78-90. Input: double U, M, the arguments. Output: double JACOBI_SN, the function value. */ { double cn; double dn; double sn; sncndn ( u, m, &sn, &cn, &dn ); return sn; } /******************************************************************************/ void sncndn ( double u, double m, double *sn, double *cn, double *dn ) /******************************************************************************/ /* Purpose: sncndn() evaluates Jacobi elliptic functions SN, CN, and DN. Discussion: Evaluation a = - *cn / *sn; corrected to a = *cn / *sn; after comparing to reference, 26 June 2025. Licensing: This code is distributed under the MIT license. Modified: 26 June 2025 Author: Original ALGOL version by Roland Bulirsch. This version by John Burkardt Reference: Roland Bulirsch, Numerical calculation of elliptic integrals and elliptic functions, Numerische Mathematik, Volume 7, Number 1, 1965, pages 78-90. Input: double U, M, the arguments. Output: double *SN, *CN, *DN, the value of the Jacobi elliptic functions sn(u,m), cn(u,m), and dn(u,m). */ { # define MAXIT 25 double a; double b; double c; double ca; double d; double m_array[MAXIT]; double n_array[MAXIT]; int i; int l; double m_comp; double u_copy; m_comp = 1.0 - m; u_copy = u; if ( m_comp == 0.0 ) { *cn = 1.0 / cosh ( u_copy ); *dn = *cn; *sn = tanh ( u_copy ); return; } if ( 1.0 < m ) { d = 1.0 - m_comp; m_comp = - m_comp / d; d = sqrt ( d ); u_copy = d * u_copy; } ca = sqrt ( DBL_EPSILON ); a = 1.0; *dn = 1.0; l = MAXIT; for ( i = 1; i <= MAXIT; i++ ) { m_array[i-1] = a; m_comp = sqrt ( m_comp ); n_array[i-1] = m_comp; c = 0.5 * ( a + m_comp ); if ( fabs ( a - m_comp ) <= ca * a ) { l = i; break; } m_comp = a * m_comp; a = c; } u_copy = c * u_copy; *sn = sin ( u_copy ); *cn = cos ( u_copy ); if ( *sn != 0.0 ) { a = *cn / *sn; c = a * c; for ( i = l; 1 <= i; i-- ) { b = m_array[i-1]; a = c * a; c = *dn * c; *dn = ( n_array[i-1] + a ) / ( b + a ); a = c / b; } a = 1.0 / sqrt ( c * c + 1.0 ); if ( *sn < 0.0 ) { *sn = - a; } else { *sn = a; } *cn = c * *sn; } if ( 1.0 < m ) { a = *dn; *dn = *cn; *cn = a; *sn = *sn / d; } return; # undef MAXIT } /******************************************************************************/ void jacobi_cn_values ( int *n_data, double *u, double *a, double *k, double *m, double *fx ) /******************************************************************************/ /* Purpose: jacobi_cn_values() returns some values of the Jacobi elliptic function CN(U,M). Discussion: In Mathematica, the function can be evaluated by: JacobiCN[ u, m ] Licensing: This code is distributed under the MIT license. Modified: 19 November 2020 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *U, the argument of the function. double *A, *K, *M, the parameters of the function. double *FX, the value of the function. */ { # define N_MAX 20 static double m_vec[N_MAX] = { 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00 }; static double fx_vec[N_MAX] = { 0.9950041652780258E+00, 0.9800665778412416E+00, 0.8775825618903727E+00, 0.5403023058681397E+00, -0.4161468365471424E+00, 0.9950124626090582E+00, 0.9801976276784098E+00, 0.8822663948904403E+00, 0.5959765676721407E+00, -0.1031836155277618E+00, 0.9950207489532265E+00, 0.9803279976447253E+00, 0.8868188839700739E+00, 0.6480542736638854E+00, 0.2658022288340797E+00, 0.3661899347368653E-01, 0.9803279976447253E+00, 0.8868188839700739E+00, 0.6480542736638854E+00, 0.2658022288340797E+00 }; static double u_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 4.0E+00, -0.2E+00, -0.5E+00, -1.0E+00, -2.0E+00 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *a = 0.0; *k = 0.0; *m = 0.0; *u = 0.0; *fx = 0.0; } else { *m = m_vec[*n_data-1]; *k = sqrt ( *m ); *a = asin ( *k ); *u = u_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ void jacobi_dn_values ( int *n_data, double *u, double *a, double *k, double *m, double *fx ) /******************************************************************************/ /* Purpose: jacobi_dn_values() returns some values of the Jacobi elliptic function DN(U,M). Discussion: In Mathematica, the function can be evaluated by: JacobiDN[ u, m ] Licensing: This code is distributed under the MIT license. Modified: 19 November 2020 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *U, the argument of the function. double *A, *K, *M, the parameters of the function. double *FX, the value of the function. */ { # define N_MAX 20 static double m_vec[N_MAX] = { 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00 }; static double fx_vec[N_MAX] = { 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.9975093485144243E+00, 0.9901483195224800E+00, 0.9429724257773857E+00, 0.8231610016315963E+00, 0.7108610477840873E+00, 0.9950207489532265E+00, 0.9803279976447253E+00, 0.8868188839700739E+00, 0.6480542736638854E+00, 0.2658022288340797E+00, 0.3661899347368653E-01, 0.9803279976447253E+00, 0.8868188839700739E+00, 0.6480542736638854E+00, 0.2658022288340797E+00 }; static double u_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 4.0E+00, -0.2E+00, -0.5E+00, -1.0E+00, -2.0E+00 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *a = 0.0; *k = 0.0; *m = 0.0; *u = 0.0; *fx = 0.0; } else { *m = m_vec[*n_data-1]; *k = sqrt ( *m ); *a = asin ( *k ); *u = u_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ void jacobi_sn_values ( int *n_data, double *u, double *a, double *k, double *m, double *fx ) /******************************************************************************/ /* Purpose: jacobi_sn_values() returns some values of the Jacobi elliptic function SN(U,M). Discussion: In Mathematica, the function can be evaluated by: JacobiSN[ u, m ] Licensing: This code is distributed under the MIT license. Modified: 19 November 2020 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *U, the argument of the function. double *A, *K, *M, the parameters of the function. double *FX, the value of the function. */ { # define N_MAX 20 static double m_vec[N_MAX] = { 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00 }; static double fx_vec[N_MAX] = { 0.9983341664682815E-01, 0.1986693307950612E+00, 0.4794255386042030E+00, 0.8414709848078965E+00, 0.9092974268256817E+00, 0.9975068547462484E-01, 0.1980217429819704E+00, 0.4707504736556573E+00, 0.8030018248956439E+00, 0.9946623253580177E+00, 0.9966799462495582E-01, 0.1973753202249040E+00, 0.4621171572600098E+00, 0.7615941559557649E+00, 0.9640275800758169E+00, 0.9993292997390670E+00, -0.1973753202249040E+00, -0.4621171572600098E+00, -0.7615941559557649E+00, -0.9640275800758169E+00 }; static double u_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 0.1E+00, 0.2E+00, 0.5E+00, 1.0E+00, 2.0E+00, 4.0E+00, -0.2E+00, -0.5E+00, -1.0E+00, -2.0E+00 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *a = 0.0; *k = 0.0; *m = 0.0; *u = 0.0; *fx = 0.0; } else { *m = m_vec[*n_data-1]; *k = sqrt ( *m ); *a = asin ( *k ); *u = u_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX }