home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ellf / ellpj.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  3.0 KB  |  153 lines
  1. /*                            ellpj.c
  2.  *
  3.  *    Jacobian Elliptic Functions
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double u, m, sn, cn, dn, phi, ellpj();
  10.  *
  11.  * ellpj( u, m, _&sn, _&cn, _&dn, _&phi );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  *
  18.  * Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
  19.  * and dn(u|m) of parameter m between 0 and 1, and real
  20.  * argument u.
  21.  *
  22.  * These functions are periodic, with quarter-period on the
  23.  * real axis equal to the complete elliptic integral
  24.  * ellpk(1.0-m).
  25.  *
  26.  * Relation to incomplete elliptic integral:
  27.  * If u = ellik(phi,m), then sn(u|m) = sin(phi),
  28.  * and cn(u|m) = cos(phi).  Phi is called the amplitude of u.
  29.  *
  30.  * Computation is by means of the arithmetic-geometric mean
  31.  * algorithm, except when m is within 1e-9 of 0 or 1.  In the
  32.  * latter case with m close to 1, the approximation applies
  33.  * only for phi < pi/2.
  34.  *
  35.  * ACCURACY:
  36.  *
  37.  * Tested at random points with u between 0 and 10, m between
  38.  * 0 and 1.
  39.  *
  40.  *            Absolute error (* = relative error):
  41.  * arithmetic   function   # trials      peak         rms
  42.  *    DEC       sn           1800       4.5e-16     8.7e-17
  43.  *    IEEE      phi         10000       9.2e-16*    1.4e-16*
  44.  *    IEEE      sn          50000       4.1e-15     4.6e-16
  45.  *    IEEE      cn          40000       3.6e-15     4.4e-16
  46.  *    IEEE      dn          10000       1.3e-12     1.8e-14
  47.  *
  48.  *  Peak error observed in consistency check using addition
  49.  * theorem for sn(u+v) was 4e-16 (absolute).  Also tested by
  50.  * the above relation to the incomplete elliptic integral.
  51.  * Accuracy deteriorates when u is large.
  52.  *
  53.  */
  54.  
  55. /*                            ellpj.c        */
  56.  
  57.  
  58. /*
  59. Cephes Math Library Release 2.0:  April, 1987
  60. Copyright 1984, 1987 by Stephen L. Moshier
  61. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  62. */
  63.  
  64.  
  65. extern double PIO2, MACHEP;
  66.  
  67. double ellpj( u, m, sn, cn, dn, ph )
  68. double u, m;
  69. double *sn, *cn, *dn, *ph;
  70. {
  71. double ai, b, phi, t, twon;
  72. double sqrt(), fabs(), sin(), cos(), asin(), tanh();
  73. double sinh(), cosh(), atan(), exp();
  74. double a[9], c[9];
  75. int i;
  76.  
  77.  
  78. /* Check for special cases */
  79.  
  80. if( m < 0.0 || m > 1.0 )
  81.     return( puts( "ellpj m out of range" ) );
  82.  
  83. if( m < 1.0e-9 )
  84.     {
  85.     t = sin(u);
  86.     b = cos(u);
  87.     ai = 0.25 * m * (u - t*b);
  88.     *sn = t - ai*b;
  89.     *cn = b + ai*t;
  90.     *ph = u - ai;
  91.     *dn = 1.0 - 0.5*m*t*t;
  92.     return;
  93.     }
  94.  
  95. if( m >= 0.9999999999 )
  96.     {
  97.     ai = 0.25 * (1.0-m);
  98.     b = cosh(u);
  99.     t = tanh(u);
  100.     phi = 1.0/b;
  101.     twon = b * sinh(u);
  102.     *sn = t + ai * (twon - u)/(b*b);
  103.     *ph = 2.0*atan(exp(u)) - PIO2 + ai*(twon - u)/b;
  104.     ai *= t * phi;
  105.     *cn = phi - ai * (twon - u);
  106.     *dn = phi + ai * (twon + u);
  107.     return;
  108.     }
  109.  
  110.  
  111. /*    A. G. M. scale        */
  112. a[0] = 1.0;
  113. b = sqrt(1.0 - m);
  114. c[0] = sqrt(m);
  115. twon = 1.0;
  116. i = 0;
  117.  
  118. while( fabs(c[i]/a[i]) > MACHEP )
  119.     {
  120.     if( i > 7 )
  121.         {
  122.         printf( "ellpj() array full" );
  123.         goto done;
  124.         }
  125.     ai = a[i];
  126.     ++i;
  127.     c[i] = ( ai - b )/2.0;
  128.     t = sqrt( ai * b );
  129.     a[i] = ( ai + b )/2.0;
  130.     b = t;
  131.     twon *= 2.0;
  132.     }
  133.  
  134. done:
  135.  
  136. /* backward recurrence */
  137. phi = twon * a[i] * u;
  138. do
  139.     {
  140.     t = c[i] * sin(phi) / a[i];
  141.     b = phi;
  142.     phi = (asin(t) + phi)/2.0;
  143.     }
  144. while( --i );
  145.  
  146. *sn = sin(phi);
  147. t = cos(phi);
  148. *cn = t;
  149. *dn = t/cos(phi-b);
  150. *ph = phi;
  151. return;
  152. }
  153.