Return to elco2.c CVS log

Up to [Research Unix] / researchv10no / cmd / map / libmap

Annotation of researchv10no/cmd/map/libmap/elco2.c, revision 1.1.1.1

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.