123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132 |
- #include <u.h>
- #include <libc.h>
- #include "map.h"
- /* elliptic integral routine, R.Bulirsch,
- * Numerische Mathematik 7(1965) 78-90
- * calculate integral from 0 to x+iy of
- * (a+b*t^2)/((1+t^2)*sqrt((1+t^2)*(1+kc^2*t^2)))
- * yields about D valid figures, where CC=10e-D
- * for a*b>=0, except at branchpoints x=0,y=+-i,+-i/kc;
- * there the accuracy may be reduced.
- * fails for kc=0 or x<0
- * return(1) for success, return(0) for fail
- *
- * special case a=b=1 is equivalent to
- * standard elliptic integral of first kind
- * from 0 to atan(x+iy) of
- * 1/sqrt(1-k^2*(sin(t))^2) where k^2=1-kc^2
- */
- #define ROOTINF 10.e18
- #define CC 1.e-6
- int
- elco2(double x, double y, double kc, double a, double b, double *u, double *v)
- {
- double c,d,dn1,dn2,e,e1,e2,f,f1,f2,h,k,m,m1,m2,sy;
- double d1[13],d2[13];
- int i,l;
- if(kc==0||x<0)
- return(0);
- sy = y>0? 1: y==0? 0: -1;
- y = fabs(y);
- csq(x,y,&c,&e2);
- d = kc*kc;
- k = 1-d;
- e1 = 1+c;
- cdiv2(1+d*c,d*e2,e1,e2,&f1,&f2);
- f2 = -k*x*y*2/f2;
- csqr(f1,f2,&dn1,&dn2);
- if(f1<0) {
- f1 = dn1;
- dn1 = -dn2;
- dn2 = -f1;
- }
- if(k<0) {
- dn1 = fabs(dn1);
- dn2 = fabs(dn2);
- }
- c = 1+dn1;
- cmul(e1,e2,c,dn2,&f1,&f2);
- cdiv(x,y,f1,f2,&d1[0],&d2[0]);
- h = a-b;
- d = f = m = 1;
- kc = fabs(kc);
- e = a;
- a += b;
- l = 4;
- for(i=1;;i++) {
- m1 = (kc+m)/2;
- m2 = m1*m1;
- k *= f/(m2*4);
- b += e*kc;
- e = a;
- cdiv2(kc+m*dn1,m*dn2,c,dn2,&f1,&f2);
- csqr(f1/m1,k*dn2*2/f2,&dn1,&dn2);
- cmul(dn1,dn2,x,y,&f1,&f2);
- x = fabs(f1);
- y = fabs(f2);
- a += b/m1;
- l *= 2;
- c = 1 +dn1;
- d *= k/2;
- cmul(x,y,x,y,&e1,&e2);
- k *= k;
- cmul(c,dn2,1+e1*m2,e2*m2,&f1,&f2);
- cdiv(d*x,d*y,f1,f2,&d1[i],&d2[i]);
- if(k<=CC)
- break;
- kc = sqrt(m*kc);
- f = m2;
- m = m1;
- }
- f1 = f2 = 0;
- for(;i>=0;i--) {
- f1 += d1[i];
- f2 += d2[i];
- }
- x *= m1;
- y *= m1;
- cdiv2(1-y,x,1+y,-x,&e1,&e2);
- e2 = x*2/e2;
- d = a/(m1*l);
- *u = atan2(e2,e1);
- if(*u<0)
- *u += PI;
- a = d*sy/2;
- *u = d*(*u) + f1*h;
- *v = (-1-log(e1*e1+e2*e2))*a + f2*h*sy + a;
- return(1);
- }
- void
- cdiv2(double c1, double c2, double d1, double d2, double *e1, double *e2)
- {
- double t;
- if(fabs(d2)>fabs(d1)) {
- t = d1, d1 = d2, d2 = t;
- t = c1, c1 = c2, c2 = t;
- }
- if(fabs(d1)>ROOTINF)
- *e2 = ROOTINF*ROOTINF;
- else
- *e2 = d1*d1 + d2*d2;
- t = d2/d1;
- *e1 = (c1+t*c2)/(d1+t*d2); /* (c1*d1+c2*d2)/(d1*d1+d2*d2) */
- }
- /* complex square root of |x|+iy */
- void
- csqr(double c1, double c2, double *e1, double *e2)
- {
- double r2;
- r2 = c1*c1 + c2*c2;
- if(r2<=0) {
- *e1 = *e2 = 0;
- return;
- }
- *e1 = sqrt((sqrt(r2) + fabs(c1))/2);
- *e2 = c2/(*e1*2);
- }
|