\ elip12     Complete Elliptic Integrals        ACM Algorithms #55 and # 56

\ Forth Scientific Library Algorithm #28

\ Evaluates the Complete Elliptic Integral of the first kind,
\     K[k] = int_0^{\pi/2} 1/Sqrt{ 1 - k^2 Sin^2(v)} dv
\ and of the second kind,
\     E[k] = int_0^{\pi/2}   Sqrt{ 1 - k^2 Sin^2(v)} dv

\ Note:
\       Uses the MODULUS k  (the parameter m = k^2).
\       These algorithms are not suitable for k = 1, and the accuracy
\         breaks down very near k = 1 ( 0.999999 )
\       These evaluations are by polynomial expansions, the accuracy is
\         controlled by the polynomial coefficients to about 7 places.
 
\ This is an ANS Forth program requiring:
\      1. The Floating-Point word set
\      3. Uses the words 'FLOAT' and Array to create floating point arrays.
\      4. The word '}' to dereference a one-dimensional array.
\      5. Uses the word '}Horner' (FSL #3) for fast polynomial evaluation.
\      6. The compilation of the test code is controlled by
\         the conditional compilation words in the Programming-Tools wordset
\      7. This program requires a FLOATING-STACK of at least 16 floats.

\ Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
\ 1980; Association for Computing Machinery Inc., New York,
\ ISBN 0-89791-017-6

\ (c) Copyright 1994 Everett F. Carter.  Permission is granted by the
\ author to use this software for any application provided this
\ copyright notice is preserved.

\ ( ELIP12     V1.1                   3 December 1994   EFC )

\ 16 constant floating-stack
\ include lib/ansfloat.4th
\ include lib/fsl-util.4th

[UNDEFINED] K[k]      [IF]
[UNDEFINED] }horner   [IF] include lib/horner.4th  [THEN]
[UNDEFINED] fln       [IF] include lib/flnflog.4th [THEN]


4 FLOAT [*] 1 [+] ARRAY K-Coefs1{
4 FLOAT [*] 1 [+] ARRAY K-Coefs2{
3 FLOAT [*] 1 [+] ARRAY E-Coefs1{
4 FLOAT [*] 1 [+] ARRAY E-Coefs2{
FLOAT K-Coefs1{ FSL-ARRAY
FLOAT K-Coefs2{ FSL-ARRAY
FLOAT E-Coefs1{ FSL-ARRAY
FLOAT E-Coefs2{ FSL-ARRAY
:THIS K-Coefs1{ DOES> (FSL-ARRAY) ;
:THIS K-Coefs2{ DOES> (FSL-ARRAY) ;
:THIS E-Coefs1{ DOES> (FSL-ARRAY) ;
:THIS E-Coefs2{ DOES> (FSL-ARRAY) ;


: init-coefs
  S" 0.5"         S>FLOAT K-Coefs1{ 0 } F!
  S" 0.12475074"  S>FLOAT K-Coefs1{ 1 } F!
  S" 0.060118519" S>FLOAT K-Coefs1{ 2 } F! 
  S" 0.010944912" S>FLOAT K-Coefs1{ 3 } F!

  S" 1.3862944"   S>FLOAT K-Coefs2{ 0 } F!
  S" 0.097932891" S>FLOAT K-Coefs2{ 1 } F!
  S" 0.054544409" S>FLOAT K-Coefs2{ 2 } F!
  S" 0.032024666" S>FLOAT K-Coefs2{ 3 } F!

  S" 0.24969795"  S>FLOAT E-Coefs1{ 0 } F!
  S" 0.08150224"  S>FLOAT E-Coefs1{ 1 } F!
  S" 0.01382999"  S>FLOAT E-Coefs1{ 2 } F!

     1            S>F     E-Coefs2{ 0 } F!
  S" 0.44479204"  S>FLOAT E-Coefs2{ 1 } F!
  S" 0.085099193" S>FLOAT E-Coefs2{ 2 } F!
  S" 0.040905094" S>FLOAT E-Coefs2{ 3 } F!
;


: K[k] ( -- ) ( F: k -- K[k] )                  \ ACM Algorithm #55

       FDUP F* 1 S>F FSWAP F-

       FDUP K-Coefs2{ 3 }Horner
       FSWAP FDUP K-Coefs1{ 3 }Horner
       FSWAP FLN F*
       F-       
;


: E[k] ( -- ) ( F: k -- K[k] )                  \ ACM Algorithm #56

       FDUP F* 1 S>F FSWAP F-

       FDUP E-Coefs2{ 3 }Horner
       FSWAP FDUP E-Coefs1{ 2 }Horner
       FOVER F* FSWAP FLN F*
       F-       
;

fclear
init-coefs

false [IF]   \ test code ==========================================
\ convert a modulus angle in degrees to the  modulus
\ : modulus   PI F* 180 S>F F/ FCOS FDUP F* 1 S>F FSWAP F- FSQRT ;

\ test driver,  calculates the complete elliptic integral of the first
\ and second kind compare with Abramowitz & Stegun,
\ Handbook of Mathematical Functions, Table 17.1


: EK_test ( -- )
        fclear
        10 set-precision

        CR
        ."  m     k         E(k) exact  K(k) exact      E(k)   K(k) " CR

       ." 0.0  0.0         1.57079633  1.57079633     "
       0 S>F FDUP E[k] F.  K[k] F. CR

      ." 0.44 0.66332495  1.38025877  1.80632756     "
      S" 0.66332495" S>FLOAT FDUP E[k] F.  K[k] F. CR

      ." 0.75 0.86602539  1.21105603  2.15651565     "
      S" 0.86602539"S>FLOAT  FDUP E[k] F.  K[k] F. CR

      ." 0.96 0.97979589  1.05050223  3.01611249     "
      S" 0.97979589" S>FLOAT FDUP E[k] F.  K[k] F. CR
;

EK_test
[THEN]

[DEFINED] 4TH# [IF]
hide K-Coefs1{
hide K-Coefs2{
hide E-Coefs1{
hide E-Coefs2{
hide init-coefs
[THEN]
[THEN]