1 28994509 2004-04-21 devnull #include <u.h>
  2 28994509 2004-04-21 devnull #include <libc.h>
  3 28994509 2004-04-21 devnull #include "map.h"
  4 28994509 2004-04-21 devnull 
  5 28994509 2004-04-21 devnull /* elliptic integral routine, R.Bulirsch,
  6 28994509 2004-04-21 devnull  *	Numerische Mathematik 7(1965) 78-90
  7 28994509 2004-04-21 devnull  *	calculate integral from 0 to x+iy of
  8 28994509 2004-04-21 devnull  *	(a+b*t^2)/((1+t^2)*sqrt((1+t^2)*(1+kc^2*t^2)))
  9 28994509 2004-04-21 devnull  *	yields about D valid figures, where CC=10e-D
 10 28994509 2004-04-21 devnull  *	for a*b>=0, except at branchpoints x=0,y=+-i,+-i/kc;
 11 28994509 2004-04-21 devnull  *	there the accuracy may be reduced.
 12 28994509 2004-04-21 devnull  *	fails for kc=0 or x<0
 13 28994509 2004-04-21 devnull  *	return(1) for success, return(0) for fail
 14 28994509 2004-04-21 devnull  *
 15 28994509 2004-04-21 devnull  *	special case a=b=1 is equivalent to
 16 28994509 2004-04-21 devnull  *	standard elliptic integral of first kind
 17 28994509 2004-04-21 devnull  *	from 0 to atan(x+iy) of
 18 28994509 2004-04-21 devnull  *	1/sqrt(1-k^2*(sin(t))^2) where k^2=1-kc^2
 19 28994509 2004-04-21 devnull */
 20 28994509 2004-04-21 devnull 
 21 28994509 2004-04-21 devnull #define ROOTINF 10.e18
 22 28994509 2004-04-21 devnull #define CC 1.e-6
 23 28994509 2004-04-21 devnull 
 24 28994509 2004-04-21 devnull int
 25 28994509 2004-04-21 devnull elco2(double x, double y, double kc, double a, double b, double *u, double *v)
 26 28994509 2004-04-21 devnull {
 27 28994509 2004-04-21 devnull 	double c,d,dn1,dn2,e,e1,e2,f,f1,f2,h,k,m,m1,m2,sy;
 28 28994509 2004-04-21 devnull 	double d1[13],d2[13];
 29 28994509 2004-04-21 devnull 	int i,l;
 30 28994509 2004-04-21 devnull 	if(kc==0||x<0)
 31 28994509 2004-04-21 devnull 		return(0);
 32 28994509 2004-04-21 devnull 	sy = y>0? 1: y==0? 0: -1;
 33 28994509 2004-04-21 devnull 	y = fabs(y);
 34 28994509 2004-04-21 devnull 	csq(x,y,&c,&e2);
 35 28994509 2004-04-21 devnull 	d = kc*kc;
 36 28994509 2004-04-21 devnull 	k = 1-d;
 37 28994509 2004-04-21 devnull 	e1 = 1+c;
 38 28994509 2004-04-21 devnull 	cdiv2(1+d*c,d*e2,e1,e2,&f1,&f2);
 39 28994509 2004-04-21 devnull 	f2 = -k*x*y*2/f2;
 40 28994509 2004-04-21 devnull 	csqr(f1,f2,&dn1,&dn2);
 41 28994509 2004-04-21 devnull 	if(f1<0) {
 42 28994509 2004-04-21 devnull 		f1 = dn1;
 43 28994509 2004-04-21 devnull 		dn1 = -dn2;
 44 28994509 2004-04-21 devnull 		dn2 = -f1;
 45 28994509 2004-04-21 devnull 	}
 46 28994509 2004-04-21 devnull 	if(k<0) {
 47 28994509 2004-04-21 devnull 		dn1 = fabs(dn1);
 48 28994509 2004-04-21 devnull 		dn2 = fabs(dn2);
 49 28994509 2004-04-21 devnull 	}
 50 28994509 2004-04-21 devnull 	c = 1+dn1;
 51 28994509 2004-04-21 devnull 	cmul(e1,e2,c,dn2,&f1,&f2);
 52 28994509 2004-04-21 devnull 	cdiv(x,y,f1,f2,&d1[0],&d2[0]);
 53 28994509 2004-04-21 devnull 	h = a-b;
 54 28994509 2004-04-21 devnull 	d = f = m = 1;
 55 28994509 2004-04-21 devnull 	kc = fabs(kc);
 56 28994509 2004-04-21 devnull 	e = a;
 57 28994509 2004-04-21 devnull 	a += b;
 58 28994509 2004-04-21 devnull 	l = 4;
 59 28994509 2004-04-21 devnull 	for(i=1;;i++) {
 60 28994509 2004-04-21 devnull 		m1 = (kc+m)/2;
 61 28994509 2004-04-21 devnull 		m2 = m1*m1;
 62 28994509 2004-04-21 devnull 		k *= f/(m2*4);
 63 28994509 2004-04-21 devnull 		b += e*kc;
 64 28994509 2004-04-21 devnull 		e = a;
 65 28994509 2004-04-21 devnull 		cdiv2(kc+m*dn1,m*dn2,c,dn2,&f1,&f2);
 66 28994509 2004-04-21 devnull 		csqr(f1/m1,k*dn2*2/f2,&dn1,&dn2);
 67 28994509 2004-04-21 devnull 		cmul(dn1,dn2,x,y,&f1,&f2);
 68 28994509 2004-04-21 devnull 		x = fabs(f1);
 69 28994509 2004-04-21 devnull 		y = fabs(f2);
 70 28994509 2004-04-21 devnull 		a += b/m1;
 71 28994509 2004-04-21 devnull 		l *= 2;
 72 28994509 2004-04-21 devnull 		c = 1 +dn1;
 73 28994509 2004-04-21 devnull 		d *= k/2;
 74 28994509 2004-04-21 devnull 		cmul(x,y,x,y,&e1,&e2);
 75 28994509 2004-04-21 devnull 		k *= k;
 76 28994509 2004-04-21 devnull 
 77 28994509 2004-04-21 devnull 		cmul(c,dn2,1+e1*m2,e2*m2,&f1,&f2);
 78 28994509 2004-04-21 devnull 		cdiv(d*x,d*y,f1,f2,&d1[i],&d2[i]);
 79 28994509 2004-04-21 devnull 		if(k<=CC) 
 80 28994509 2004-04-21 devnull 			break;
 81 28994509 2004-04-21 devnull 		kc = sqrt(m*kc);
 82 28994509 2004-04-21 devnull 		f = m2;
 83 28994509 2004-04-21 devnull 		m = m1;
 84 28994509 2004-04-21 devnull 	}
 85 28994509 2004-04-21 devnull 	f1 = f2 = 0;
 86 28994509 2004-04-21 devnull 	for(;i>=0;i--) {
 87 28994509 2004-04-21 devnull 		f1 += d1[i];
 88 28994509 2004-04-21 devnull 		f2 += d2[i];
 89 28994509 2004-04-21 devnull 	}
 90 28994509 2004-04-21 devnull 	x *= m1;
 91 28994509 2004-04-21 devnull 	y *= m1;
 92 28994509 2004-04-21 devnull 	cdiv2(1-y,x,1+y,-x,&e1,&e2);
 93 28994509 2004-04-21 devnull 	e2 = x*2/e2;
 94 28994509 2004-04-21 devnull 	d = a/(m1*l);
 95 28994509 2004-04-21 devnull 	*u = atan2(e2,e1);
 96 28994509 2004-04-21 devnull 	if(*u<0)
 97 28994509 2004-04-21 devnull 		*u += PI;
 98 28994509 2004-04-21 devnull 	a = d*sy/2;
 99 28994509 2004-04-21 devnull 	*u = d*(*u) + f1*h;
100 28994509 2004-04-21 devnull 	*v = (-1-log(e1*e1+e2*e2))*a + f2*h*sy + a;
101 28994509 2004-04-21 devnull 	return(1);
102 28994509 2004-04-21 devnull }
103 28994509 2004-04-21 devnull 
104 28994509 2004-04-21 devnull void
105 28994509 2004-04-21 devnull cdiv2(double c1, double c2, double d1, double d2, double *e1, double *e2)
106 28994509 2004-04-21 devnull {
107 28994509 2004-04-21 devnull 	double t;
108 28994509 2004-04-21 devnull 	if(fabs(d2)>fabs(d1)) {
109 28994509 2004-04-21 devnull 		t = d1, d1 = d2, d2 = t;
110 28994509 2004-04-21 devnull 		t = c1, c1 = c2, c2 = t;
111 28994509 2004-04-21 devnull 	}
112 28994509 2004-04-21 devnull 	if(fabs(d1)>ROOTINF)
113 28994509 2004-04-21 devnull 		*e2 = ROOTINF*ROOTINF;
114 28994509 2004-04-21 devnull 	else
115 28994509 2004-04-21 devnull 		*e2 = d1*d1 + d2*d2;
116 28994509 2004-04-21 devnull 	t = d2/d1;
117 28994509 2004-04-21 devnull 	*e1 = (c1+t*c2)/(d1+t*d2); /* (c1*d1+c2*d2)/(d1*d1+d2*d2) */
118 28994509 2004-04-21 devnull }
119 28994509 2004-04-21 devnull 
120 28994509 2004-04-21 devnull /* complex square root of |x|+iy */
121 28994509 2004-04-21 devnull void
122 28994509 2004-04-21 devnull csqr(double c1, double c2, double *e1, double *e2)
123 28994509 2004-04-21 devnull {
124 28994509 2004-04-21 devnull 	double r2;
125 28994509 2004-04-21 devnull 	r2 = c1*c1 + c2*c2;
126 28994509 2004-04-21 devnull 	if(r2<=0) {
127 28994509 2004-04-21 devnull 		*e1 = *e2 = 0;
128 28994509 2004-04-21 devnull 		return;
129 28994509 2004-04-21 devnull 	}
130 28994509 2004-04-21 devnull 	*e1 = sqrt((sqrt(r2) + fabs(c1))/2);
131 28994509 2004-04-21 devnull 	*e2 = c2/(*e1*2);
132 28994509 2004-04-21 devnull }