elco2.c (2532B)
1 #include <u.h> 2 #include <libc.h> 3 #include "map.h" 4 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 14 * 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 19 */ 20 21 #define ROOTINF 10.e18 22 #define CC 1.e-6 23 24 int 25 elco2(double x, double y, double kc, double a, double b, double *u, double *v) 26 { 27 double c,d,dn1,dn2,e,e1,e2,f,f1,f2,h,k,m,m1,m2,sy; 28 double d1[13],d2[13]; 29 int i,l; 30 if(kc==0||x<0) 31 return(0); 32 sy = y>0? 1: y==0? 0: -1; 33 y = fabs(y); 34 csq(x,y,&c,&e2); 35 d = kc*kc; 36 k = 1-d; 37 e1 = 1+c; 38 cdiv2(1+d*c,d*e2,e1,e2,&f1,&f2); 39 f2 = -k*x*y*2/f2; 40 csqr(f1,f2,&dn1,&dn2); 41 if(f1<0) { 42 f1 = dn1; 43 dn1 = -dn2; 44 dn2 = -f1; 45 } 46 if(k<0) { 47 dn1 = fabs(dn1); 48 dn2 = fabs(dn2); 49 } 50 c = 1+dn1; 51 cmul(e1,e2,c,dn2,&f1,&f2); 52 cdiv(x,y,f1,f2,&d1[0],&d2[0]); 53 h = a-b; 54 d = f = m = 1; 55 kc = fabs(kc); 56 e = a; 57 a += b; 58 l = 4; 59 for(i=1;;i++) { 60 m1 = (kc+m)/2; 61 m2 = m1*m1; 62 k *= f/(m2*4); 63 b += e*kc; 64 e = a; 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); 68 x = fabs(f1); 69 y = fabs(f2); 70 a += b/m1; 71 l *= 2; 72 c = 1 +dn1; 73 d *= k/2; 74 cmul(x,y,x,y,&e1,&e2); 75 k *= k; 76 77 cmul(c,dn2,1+e1*m2,e2*m2,&f1,&f2); 78 cdiv(d*x,d*y,f1,f2,&d1[i],&d2[i]); 79 if(k<=CC) 80 break; 81 kc = sqrt(m*kc); 82 f = m2; 83 m = m1; 84 } 85 f1 = f2 = 0; 86 for(;i>=0;i--) { 87 f1 += d1[i]; 88 f2 += d2[i]; 89 } 90 x *= m1; 91 y *= m1; 92 cdiv2(1-y,x,1+y,-x,&e1,&e2); 93 e2 = x*2/e2; 94 d = a/(m1*l); 95 *u = atan2(e2,e1); 96 if(*u<0) 97 *u += PI; 98 a = d*sy/2; 99 *u = d*(*u) + f1*h; 100 *v = (-1-log(e1*e1+e2*e2))*a + f2*h*sy + a; 101 return(1); 102 } 103 104 void 105 cdiv2(double c1, double c2, double d1, double d2, double *e1, double *e2) 106 { 107 double t; 108 if(fabs(d2)>fabs(d1)) { 109 t = d1, d1 = d2, d2 = t; 110 t = c1, c1 = c2, c2 = t; 111 } 112 if(fabs(d1)>ROOTINF) 113 *e2 = ROOTINF*ROOTINF; 114 else 115 *e2 = d1*d1 + d2*d2; 116 t = d2/d1; 117 *e1 = (c1+t*c2)/(d1+t*d2); /* (c1*d1+c2*d2)/(d1*d1+d2*d2) */ 118 } 119 120 /* complex square root of |x|+iy */ 121 void 122 csqr(double c1, double c2, double *e1, double *e2) 123 { 124 double r2; 125 r2 = c1*c1 + c2*c2; 126 if(r2<=0) { 127 *e1 = *e2 = 0; 128 return; 129 } 130 *e1 = sqrt((sqrt(r2) + fabs(c1))/2); 131 *e2 = c2/(*e1*2); 132 }