Re: elliptic integrals

[Vince Hradil <hradilv(at)yahoo.com>]

I wrote this using Numerical Recipes in C.  NB: I just wrote this last week,
so I haven't tested it all out yet.  Check NumRec for the reference to the
algorithm.

FUNCTION cel, qqc, pp, aa, bb
  ;; from numerical recipes in C p. 201

  ca = 0.0003D
  pio2 = !DPI/2.0

  IF ( qqc EQ 0) THEN BEGIN
    print, 'Bad qqc in cel'
    return, 0
  ENDIF

  qc = abs(qqc)
  a = aa
  b = bb
  p = pp
  e = qc
  em = 1.0D

  IF ( p GT 0.0 ) THEN BEGIN
    p = sqrt(p)
    b = b/p
  ENDIF ELSE BEGIN
    f = qc*qc
    q = 1.0-f
    g = 1.0-p
    f = f-p
    q = q*(b-a*p)
    p = sqrt(f/g)
    a = (a-b)/g
    b = -q/(g*g*p)+a*p
  ENDELSE

  WHILE ( 1 ) DO BEGIN
    f = a
    a = a+(b/p)
    g = e/p
    b = b+(f*g)
    b = b+b
    p = g+p
    g = em
    em = em+qc
    IF ( abs(g-qc) LE g*ca ) THEN BEGIN
      rval = pio2*(b+a*em)/(em*(em+p))
      return, rval
    ENDIF
    qc = sqrt(e)
    qc = qc+qc
    e = qc*em
  ENDWHILE

END

FUNCTION ellint1, k

  kc = sqrt(1.0-double(k))

  e1 = cel(kc, 1.0D, 1.0D, 1.0D)
  return, e1

END

FUNCTION ellint2, k

  kc = sqrt(1.0-double(k))

  e2 = cel(kc, 1.0D, 1.0D, kc*kc)
  return, e2

END



Tom Intrator wrote:

> do  you have any thing that would compute elliptic integrals?
>
> I need the complete ones K,E but for arguments m>1 (terminology as per
> Abramowitz and Stegun) This requires computation of some incomplete ones
> as well
>
> I am using IDL for this
>
> thanks
>
> Tom Intrator