| Re: elliptic integrals [message #18631] |
Wed, 26 January 2000 00:00 |
Vince Hradil
Messages: 574
Registered: December 1999
|
Senior Member |
|
|
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
|
|
|
|
| Re: elliptic integrals [message #18632 is a reply to message #18631] |
Wed, 26 January 2000 00:00  |
Vince Hradil
Messages: 574
Registered: December 1999
|
Senior Member |
|
|
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
|
|
|
|
| Re: elliptic integrals [message #18633 is a reply to message #18631] |
Wed, 26 January 2000 00:00 |
Vince Hradil
Messages: 574
Registered: December 1999
|
Senior Member |
|
|
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
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
