Return to elco2.c CVS log

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

Annotation of researchv8dc/cmd/map/libmap/elco2.c, revision 1.1

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

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