sys/src/cmd/map/libmap/elco2.c

   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 }