Warning: Attempt to read property "date" on null in /usr/local/www/websvn.planix.org/blame.php on line 247

Warning: Attempt to read property "msg" on null in /usr/local/www/websvn.planix.org/blame.php on line 247
WebSVN – planix.SVN – Blame – /os/branches/feature_posix/sys/src/cmd/map/libmap/elco2.c – Rev 2

Subversion Repositories planix.SVN

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 
}