5 /* elliptic integral routine, R.Bulirsch,
6 * Numerische Mathematik 7(1965) 78-90
7 * calculate integral from 0 to x+iy of
8 * (a+b*t^2)/((1+t^2)*sqrt((1+t^2)*(1+kc^2*t^2)))
9 * yields about D valid figures, where CC=10e-D
10 * for a*b>=0, except at branchpoints x=0,y=+-i,+-i/kc;
11 * there the accuracy may be reduced.
12 * fails for kc=0 or x<0
13 * return(1) for success, return(0) for fail
15 * special case a=b=1 is equivalent to
16 * standard elliptic integral of first kind
17 * from 0 to atan(x+iy) of
18 * 1/sqrt(1-k^2*(sin(t))^2) where k^2=1-kc^2
21 #define ROOTINF 10.e18
25 elco2(double x, double y, double kc, double a, double b, double *u, double *v)
27 double c,d,dn1,dn2,e,e1,e2,f,f1,f2,h,k,m,m1,m2,sy;
32 sy = y>0? 1: y==0? 0: -1;
38 cdiv2(1+d*c,d*e2,e1,e2,&f1,&f2);
40 csqr(f1,f2,&dn1,&dn2);
51 cmul(e1,e2,c,dn2,&f1,&f2);
52 cdiv(x,y,f1,f2,&d1[0],&d2[0]);
65 cdiv2(kc+m*dn1,m*dn2,c,dn2,&f1,&f2);
66 csqr(f1/m1,k*dn2*2/f2,&dn1,&dn2);
67 cmul(dn1,dn2,x,y,&f1,&f2);
74 cmul(x,y,x,y,&e1,&e2);
77 cmul(c,dn2,1+e1*m2,e2*m2,&f1,&f2);
78 cdiv(d*x,d*y,f1,f2,&d1[i],&d2[i]);
92 cdiv2(1-y,x,1+y,-x,&e1,&e2);
100 *v = (-1-log(e1*e1+e2*e2))*a + f2*h*sy + a;
105 cdiv2(double c1, double c2, double d1, double d2, double *e1, double *e2)
108 if(fabs(d2)>fabs(d1)) {
109 t = d1, d1 = d2, d2 = t;
110 t = c1, c1 = c2, c2 = t;
113 *e2 = ROOTINF*ROOTINF;
117 *e1 = (c1+t*c2)/(d1+t*d2); /* (c1*d1+c2*d2)/(d1*d1+d2*d2) */
120 /* complex square root of |x|+iy */
122 csqr(double c1, double c2, double *e1, double *e2)
130 *e1 = sqrt((sqrt(r2) + fabs(c1))/2);