Annotation of researchv10no/cmd/ideal/inter.c, revision 1.1.1.1

1.1       root        1: #include "ideal.h"
                      2: #include "y.tab.h"
                      3: 
                      4: boolean llinter (x1,y1, x2,y2, z1,w1, z2,w2, alpha, beta, collinear)
                      5:        /* if this function returns TRUE,
                      6:        /* then alpha[(x1,y1),(x2,y2)] = beta[(z1,w1),(z2,w2)] */
                      7: float x1,y1, x2,y2, z1,w1, z2,w2;
                      8: float *alpha;
                      9: float *beta;
                     10: boolean *collinear;
                     11: {
                     12:        float A1, B1, C1, A2, B2, C2, D, x, y;
                     13:        dprintf "(%f,%f) -- (%f,%f)\n", x1,y1, x2,y2);
                     14:        dprintf "(%f,%f) -- (%f,%f)\n", z1,w1, z2,w2);
                     15:        A1 = y1 - y2;
                     16:        B1 = x2 - x1;
                     17:        C1 = -B1*y1 - A1*x1;
                     18:        A2 = w1 - w2;
                     19:        B2 = z2 - z1;
                     20:        C2 = -B2*w1 - A2*z1;
                     21:        D = A1*B2 - A2*B1;
                     22:        if (fabs(D) < EPSILON) {
                     23:                *collinear = arecollinear(x1,y1,x2,y2,z1,w1);
                     24:                dprintf "%s\n", (*collinear)?"coincident":"disjoint");
                     25:                return (FALSE);
                     26:        }
                     27:        *collinear = FALSE;
                     28:        x = (B1*C2 - B2*C1)/D;
                     29:        y = (A2*C1 - A1*C2)/D;
                     30:        if (fabs(x2 - x1) > EPSILON) {
                     31:                *alpha = (x - x1)/(x2 - x1);
                     32:        } else if (fabs(y2 - y1) > EPSILON) {
                     33:                *alpha = (y - y1)/(y2 - y1);
                     34:        } else fprintf (stderr, "ideal: llinter: can't happen\n");
                     35:        if (fabs(z2 - z1) > EPSILON) {
                     36:                *beta = (x - z1)/(z2 - z1);
                     37:        } else if (fabs(w2 - w1) > EPSILON) {
                     38:                *beta = (y - w1)/(w2 - w1);
                     39:        } else fprintf (stderr, "ideal: llinter: can't happen\n");
                     40:        dprintf "intersection alpha = %f; beta = %f\n", *alpha, *beta);
                     41:        return (TRUE);
                     42: } /* llinter */
                     43: 
                     44: boolean lcinter (x1,y1, x2,y2, x0,y0, r, alpha1,theta1, alpha2,theta2)
                     45:        /* if this function returns TRUE,
                     46:        /* then alpha1[(x1,y1),(x2,y2)] = (x0,y0) + r*cis(theta1)
                     47:        /* and alpha2[(x1,y1),(x2,y2)] = (x0,y0) + r*cis(theta2) */
                     48: float x1,y1, x2,y2, x0,y0, r;
                     49: float *alpha1;
                     50: float *theta1;
                     51: float *alpha2;
                     52: float *theta2;
                     53: {
                     54:        float dx1, dx2, dy1, dy2;
                     55:        float A, B, C, D;
                     56: 
                     57:        dprintf "intersection parameters:\n");
                     58:        dprintf "%f, %f -- %f, %f\n", x1, y1, x2, y2);
                     59:        dprintf "%f, %f (%f)\n", x0, y0, r);
                     60:        r = fabs(r);
                     61:        dx1 = x1 - x0;
                     62:        dx2 = x2 - x1;
                     63:        dy1 = y1 - y0;
                     64:        dy2 = y2 - y1;
                     65:        A = dx2*dx2 + dy2*dy2;
                     66:        dprintf "A=%f\n", A);
                     67:        if (A < EPSILON) {
                     68:                if (fabs (hypot (dx1, dy1) - r) < EPSILON) {
                     69:                        *alpha1 = *alpha2 = 0.0;
                     70:                        *theta1 = atan2 (dy1, dx1);
                     71:                        *theta2 = *theta1 = rprin (*theta1);
                     72:                        dprintf "alpha1 = alpha2 = %f theta1 = theta2 = %f\n",
                     73:                                *alpha1, *theta1);
                     74:                        return (TRUE);
                     75:                }
                     76:                else
                     77:                        return (FALSE);
                     78:        }
                     79:        B = 2*(dx1*dx2 + dy1*dy2);
                     80:        C = dx1*dx1 + dy1*dy1 - r*r;
                     81:        D = B*B - 4*A*C;
                     82:        dprintf "B=%f C=%f D=%f\n", B, C, D);
                     83:        if (D < -EPSILON)
                     84:                return (FALSE);
                     85:        if (fabs(D) < EPSILON)
                     86:                D = 0.0;
                     87:        D = sqrt(D);
                     88:        *alpha1 = (-B + D)/(2.0*A);
                     89:        *theta1 = rprin (atan2 (dy1 + *alpha1*dy2, dx1 + *alpha1*dx2));
                     90:        *alpha2 = (-B - D)/(2.0*A);
                     91:        *theta2 = rprin (atan2 (dy1 + *alpha2*dy2, dx1 + *alpha2*dx2));
                     92:        dprintf "intersection alpha1 = %f, theta1 = %f\n", *alpha1, *theta1);
                     93:        dprintf "intersection alpha2 = %f, theta2 = %f\n", *alpha2, *theta2);
                     94:        return (TRUE);
                     95: }
                     96: 
                     97: boolean ccinter (x0,y0,r0, x1,y1,r1, theta1,phi1, theta2,phi2)
                     98:        /* if this function returns TRUE,
                     99:        /* then (x0,y0) + r0*cis(theta1) = (x1,y1) + r1*cis(phi1)
                    100:        /* and (x0,y0) + r0*cis(theta2) = (x1,y1) + r1*cis(phi2) */
                    101: float x0,y0,r0;
                    102: float x1,y1,r1;
                    103: float *theta1;
                    104: float *phi1;
                    105: float *theta2;
                    106: float *phi2;
                    107: {
                    108:        float xcoeff, ycoeff, const;
                    109:        float u1, v1, u2, v2;
                    110:        boolean lncrc;
                    111: 
                    112:        dprintf "intersection parameters\n");
                    113:        dprintf "%f %f (%f)\n", x0, y0, r0);
                    114:        dprintf "%f %f (%f)\n", x1, y1, r1);
                    115: 
                    116:        r0 = fabs(r0);
                    117:        r1 = fabs(r1);
                    118:        xcoeff = 2*(x1 - x0);
                    119:        ycoeff = 2*(y1 - y0);
                    120:        const = r0*r0 - x0*x0 - y0*y0 - r1*r1 + x1*x1 + y1*y1;
                    121:        if (fabs(xcoeff) < EPSILON && fabs(ycoeff) < EPSILON)
                    122:                return (FALSE);
                    123:        if (fabs(xcoeff) < EPSILON) {
                    124:                u1 = 0.0;
                    125:                u2 = 1.0;
                    126:                v1 = v2 = const/ycoeff;
                    127:        } else if (fabs(ycoeff) < EPSILON) {
                    128:                v1 = 0.0;
                    129:                v2 = 1.0;
                    130:                u1 = u2 = const/xcoeff;
                    131:        } else if (fabs(const) < EPSILON) {
                    132:                u1 = 0.0;
                    133:                v1 = 0.0;
                    134:                u2 = 1.0;
                    135:                v2 = (const - 1.0/xcoeff)/ycoeff;
                    136:        } else {
                    137:                u1 = 0.0;
                    138:                v1 = const/ycoeff;
                    139:                u2 = const/xcoeff;
                    140:                v2 = 0.0;
                    141:        }
                    142:        lncrc = lcinter (u1,v1, u2,v2, x1,y1,r1, theta1,phi1, theta2,phi2);
                    143:        if (lncrc) {
                    144:                *phi1 = rprin (*phi1);
                    145:                *phi2 = rprin (*phi2);
                    146:                *theta1 = atan2 (y1 + r1*sin(*phi1) - y0, x1 + r1*cos(*phi1) - x0);
                    147:                *theta2 = atan2 (y1 + r1*sin(*phi2) - y0, x1 + r1*cos(*phi2) - x0);
                    148:                *theta1 = rprin (*theta1);
                    149:                *theta2 = rprin (*theta2);
                    150:                dprintf "intersection theta1 = %f phi1 = %f\n", *theta1, *phi1);
                    151:                dprintf "intersection theta2 = %f phi2 = %f\n", *theta2, *phi2);
                    152:                return (TRUE);
                    153:        } else
                    154:                return (FALSE);
                    155: }

unix.superglobalmegacorp.com

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