[scilab-Users] Re: Elliptic integral

[Peter.Cusack at csiro.au]

Hi Jens.

Thanks for looking into this.

In fact, the code was translated from working Maxima code, and that included a check to ensure that the argument remained less than 1.

Here is the full version of my function declaration, in hope that you can see an obvious flaw:

deff('[a]=al(r,a)','a=r/a') //*alpha*/
deff('[b]=be(x,a)','b=x/a') //*beta*/
deff('[g]=ga(x,r)','g=x/r') //*gamma*/
deff('[que]=Q(x,r,a)','que=(1+al(r,a))^2 + be(x,a)^2')//may need to be square rooted
deff('[kay]=k(x,r,a)','if al(r,a)/Q(x,r,a)<0.25 then kay=4*al(r,a)/Q(x,r,a); else kay=0.999; end')//*ellipticK(1) is undefined*/
deff('[b0]=B0(a,j)','b0=j*Mu0/(2*a)') //*Field at centre of single turn with radius a and current j*/
//Elliptic integrals
deff('ek=elliptic_kc(k)','ek=delip(1,k)')//First kind is a built in function
deff('y=f(t,k)','y=sqrt((1-(k*t)^2)/(1-t^2))') //Core definition for integral
deff('ec=elliptic_ec(k)','ec=intg(0,1,list(f,k))') //Second kind, not built into SciLab
//Field calculations
//Single turns
deff('[xx]=Bx(x,r,a,j)','xx=B0(a,j)/(%pi*sqrt(Q(x,r,a)))*(elliptic_ec(k(x,r,a))*(1-al(r,a)^2-be(x,a)^2)/(Q(x,r,a)-4*al(r,a))+elliptic_kc(k(x,r,a)))') //X field
deff('[rr]=Br(x,r,a,j)','if r>0 then rr=B0(a,j)*ga(x,r)/(%pi*sqrt(Q(x,r,a)))*(elliptic_ec(k(x,r,a))*(1+al(r,a)^2+be(x,a)^2)/(Q(x,r,a)-4*al(r,a))-elliptic_kc(k(x,r,a))); else rr=0; end') //*R field. ga(x,0) is undefined*/
//Total of all turns
deff('[xt]=Bxt(x,r)','xt=sum(Bx(x-X,r,A,I))')//*x field at x,r*/
deff('[rt]=Brt(x,r)','rt=sum(Br(x-X,r,A,I))')//*r field at x,r*/
//Uniformity
Bcentre=Bxt(0,0) //*Field at centre of winding*/

Cheers;

Peter.

> ------------------------------------------------------------------------
> Hi Peter,
> You probably got the error message for arguments >1. In this case the
> sqrt
> returns an imaginary number, which cannot be handled by intg. For
> arguments <=1
> your function works  properly. However, I did not check if your code
> realizes
> the integral you intend to calculate correctly.
> Cheers, Jens
>