/*
 * FIG : Facility for Interactive Generation of figures
 * Copyright (c) 1985-1988 by Supoj Sutanthavibul
 * Parts Copyright (c) 1989-2002 by Brian V. Smith
 * Parts Copyright (c) 1991 by Paul King
 * Parts Copyright (c) 2004 by Chris Moller
 *
 * Any party obtaining a copy of these files is granted, free of charge, a
 * full and unrestricted irrevocable, world-wide, paid up, royalty-free,
 * nonexclusive right and license to deal in this software and documentation
 * files (the "Software"), including without limitation the rights to use,
 * copy, modify, merge, publish distribute, sublicense and/or sell copies of
 * the Software, and to permit persons who receive copies from any such
 * party to do so, with the only requirement being that the above copyright
 * and this permission notice remain intact.
 *
 */

#include "fig.h"
#include "figx.h"
#include "resources.h"
#include "object.h"
#include "mode.h"
#include "w_snap.h"
#include "w_intersect.h"
#include "u_quartic.h"
#include <alloca.h>
#include <math.h>
#if !defined(HAVE_NO_COMPLEX) 	       
#include <complex.h>
#endif
#undef I
  
#define ISET_P1 (1 << 0)
#define ISET_P2 (1 << 1)

intersect_state_e intersect_state = INTERSECT_INITIAL;

static INLINE Boolean
is_point_on_segment(double x1, double y1, double xc, double yc, double x2, double y2)
{
  return ((((x1 <= xc) && (xc <= x2)) || ((x1 >= xc) && (xc >= x2))) && 
	  (((y1 <= yc) && (yc <= y2)) || ((y1 >= yc) && (yc >= y2)))) ? True : False;
}

void
insert_isect(isect_cb_s * isect_cb, double ix, double iy, int seg_idx)
{
  if (isect_cb->max_isects <= isect_cb->nr_isects) {
    isect_cb->max_isects += ISECTS_INCR;
    isect_cb->isects = realloc(isect_cb->isects,
			       isect_cb->max_isects * sizeof(isects_s));
  }
  isect_cb->isects[isect_cb->nr_isects].x = (int)lrint(ix);
  isect_cb->isects[isect_cb->nr_isects].y = (int)lrint(iy);
  isect_cb->isects[isect_cb->nr_isects++].seg_idx = seg_idx;
}


/* returns coeffs of the form Ax + By + C = 0 */

static void
fget_line_from_points(double * c,
		      double s1x, double s1y,
		      double s2x, double s2y)
{
  c[0] = s1y - s2y;
  c[1] = s2x - s1x;
  c[2] = (s1x * s2y) - (s1y * s2x);
}

/* returns True if the bounding box defined by corners p1a - p1b	*/
/* overlaps the box defined by p2a - p2b				*/

static Boolean
boxes_overlap(struct f_point * p1a, struct f_point * p1b,
	      struct f_point * p2a, struct f_point * p2b)
{
  Boolean rc;
  int p1_min_x = (p1a->x < p1b->x) ? p1a->x : p1b->x;
  int p1_max_x = (p1a->x > p1b->x) ? p1a->x : p1b->x;
  int p1_min_y = (p1a->y < p1b->y) ? p1a->y : p1b->y;
  int p1_max_y = (p1a->y > p1b->y) ? p1a->y : p1b->y;
  int p2_min_x = (p2a->x < p2b->x) ? p2a->x : p2b->x;
  int p2_max_x = (p2a->x > p2b->x) ? p2a->x : p2b->x;
  int p2_min_y = (p2a->y < p2b->y) ? p2a->y : p2b->y;
  int p2_max_y = (p2a->y > p2b->y) ? p2a->y : p2b->y;
  rc = ((p1_max_x < p2_min_x) || (p1_min_x > p2_max_x) ||
	(p1_max_y < p2_min_y) || (p1_min_y > p2_max_y)) ? False : True;
  return rc;
}

static void
do_circle_ellipse_intersect(r, X, Y, e, x, y, arc, isect_cb)
     double r, X, Y;
     F_ellipse * e;
     int x, y;
     F_arc * arc;
     isect_cb_s * isect_cb;
{
  double K, L, M, N;
  double hix, hiy;
  
  double * poly;
  double * solr;
  double * soli;
  double mind = HUGE_VAL;
  int nsol, i;

  double a = (double)(e->radiuses.x);
  double b = (double)(e->radiuses.y);
  
  /* translate so center of ellipse is at origin */
  double EX = (double)(x - e->center.x);
  double EY = (double)(y - e->center.y);
  
  /* rotate around ellipse angle so ellipse semi-axes are ortho to space */
  snap_rotate_vector (&EX, &EY, EX, EY, (double)(e->angle));
  snap_rotate_vector (&X,  &Y,  X,  Y,  (double)(e->angle));

  poly = alloca(5 * sizeof(double));
  solr = alloca(4 * sizeof(double));
  soli = alloca(4 * sizeof(double));

  K = 1.0 - pow(b/a, 2.0);
  L = 2.0 * X;
  M = (pow(X, 2.0) + pow(Y, 2.0) + pow(b, 2.0)) - pow(r, 2.0);
  N = (2.0 * Y * b)/a;

  poly[4] = pow(K, 2.0);
  poly[3] = -2.0 * K * L;
  poly[2] = (2.0 * K * M) + pow(L, 2.0) + pow(N, 2.0);
  poly[1] = -2.0 * L * M;
  poly[0] = pow(M, 2.0) - pow(N * a, 2.0);

  nsol = quartic(poly, solr, soli);

  for (i = 0; i < nsol; i++) {
    if (1.0 > fabs(soli[i])) {	/* ignore complex roots */
      double dd;
      double ix, cy, ey;
      int ix_idx, cy_idx, ey_idx;
      
      for (ix_idx = 0; ix_idx < 2; ix_idx++) {			/* for +/- ix vals... */
	double rxc, rxe, dt, de;
	
	ix = (ix_idx & 1) ? -solr[i] : solr[i];
	rxc = pow(r, 2.0) - pow(ix - X, 2.0);
	rxe = pow(a, 2.0) - pow(ix, 2.0);
	if ((0.0 <= rxc) && (0.0 <= rxe)) {
	  for (cy_idx = 0; cy_idx < 2; cy_idx++) {		/* and +/- circle slns */
	    cy = Y + ((cy_idx & 1) ? -sqrt(rxc) : sqrt(rxc));
	    for (ey_idx = 0; ey_idx < 2; ey_idx++) {		/* and +/- ellipse slns */
	      Boolean point_ok;
	    
	      ey = (b/a) * ((ey_idx & 1) ? -sqrt(rxe) : sqrt(rxe));

	      if (arc) {
		double ttx, tty;
		int ittx, itty;
		snap_rotate_vector (&ttx,  &tty,  ix,  cy,  -(double)(e->angle));
		ittx = ((int)lrint(ttx)) + e->center.x;
		itty = ((int)lrint(tty)) + e->center.y;
		point_ok = is_point_on_arc(arc, ittx, itty);
	      }
	      else point_ok = True;
	      
	      if (True == point_ok) {
		dt = fabs(cy - ey);
		if (dt < 1.0) {
		  if (isect_cb) {
		    double isx, isy;
		    snap_rotate_vector (&isx,  &isy,  ix,  cy,  -(double)(e->angle));
		    isx += (double)(e->center.x);
		    isy += (double)(e->center.y);
		    insert_isect(isect_cb, isx, isy, -1);
		  }
		  dd = hypot(cy - EY, ix - EX);
		  if (dd < mind) {
		    mind = dd;
		    hix = ix;
		    hiy = cy;
		    snap_found = True;
		  }
		}
	      }
	    }
	  }
	}
      }
    }
  }

  if (True == snap_found) {
    snap_rotate_vector (&hix,  &hiy,  hix,  hiy,  -(double)(e->angle));
    snap_gx = ((int)lrint(hix)) + e->center.x;
    snap_gy = ((int)lrint(hiy)) + e->center.y;
  }
}
     
static void
circle_ellipse_intersect(c, e, x, y, isect_cb)
     F_ellipse * c;
     F_ellipse * e;
     int x, y;
     isect_cb_s * isect_cb;
{
  double r = (double)(c->radiuses.x);
  double X  = (double)(c->center.x - e->center.x);
  double Y  = (double)(c->center.y - e->center.y);

  do_circle_ellipse_intersect(r, X, Y, e, x, y, isect_cb);
}

static void
do_circle_circle(PX, PY, X2, Y2, R1, R2, OX, OY, arc, arc2, isect_cb)
     double PX;				/* event point	*/
     double PY;
     double X2;				/* translated arc center */
     double Y2;
     double R1;				/* circle radius */
     double R2;
     double OX;				/* circle -> arc displacement */
     double OY;
     F_arc * arc;
     F_arc * arc2;
     isect_cb_s * isect_cb;
{
  double dc = hypot(Y2, X2);

  if ((0.0 == dc) && (R1 == R2)) {
    put_msg("Selected circles/arcs are congruent.");
    beep();
    snap_msg_set = True;
  }
  else if (dc < (R1 + R2)) {			/* real intersection */
    double K = (pow(R1, 2.0) - pow(R2, 2.0)) + pow(dc, 2.0);
    double L = 2.0 * Y2;
    double M = 2.0 * X2;
    double A = pow(L, 2.0) + pow(M, 2.0);
    double B = -2.0 * K * L;
    double C = pow(K, 2.0) - (pow(M, 2.0) * pow(R1, 2.0));
    double rx = pow(B, 2.0) - (4.0 * A * C);
	  
    if (0.0 <= rx) {			/* should albays be true...*/
      double ix, iy;
      double hix, hiy;
      double mind = HUGE_VAL;
      int iy_idx, ix_idx, ixa_idx;

      rx = sqrt(rx)/(2.0 * A);

      for (iy_idx = 0; iy_idx < 2; iy_idx++) {
	double tix, tixa;
	iy = (-B/(2.0 * A)) + ((iy_idx & 1) ? -rx : rx);
	tix = sqrt(pow(R1, 2.0) - pow(iy, 2.0));	/* sln of circle or arc1 */
	tixa = sqrt(pow(R2, 2.0) - pow(iy - Y2, 2.0));
	for (ixa_idx = 0; ixa_idx < 2; ixa_idx++) {
	  double ixa = X2 + ((ixa_idx & 1) ? -tixa : tixa);

	  for (ix_idx = 0; ix_idx < 2; ix_idx++) {
	    double dd;
	    Boolean point_ok;
	  
	    ix = (ix_idx & 1) ? -tix : tix;
	    if (fabs(ix - ixa) < 1.0) {
	      if (arc) {
		int ittx, itty;
	  
		ittx = (int)lrint(ix + OX);
		itty = (int)lrint(iy + OY);
		point_ok = is_point_on_arc(arc, ittx, itty);

		if (arc2 && (True == point_ok))
		  point_ok = is_point_on_arc(arc2, ittx, itty);
	    
	      }
	      else point_ok = True;
	  
	      if (True == point_ok) {
		if (isect_cb) insert_isect(isect_cb, ix + OX, iy + OY, -1);
		dd = hypot(iy - PY, ix - PX);
		if (dd < mind) {
		  mind = dd;
		  hix = ix;
		  hiy = iy;
		}
	      }
	    }
	  }
	}
      }
      
      snap_gx = (int)lrint(hix + OX);
      snap_gy = (int)lrint(hiy + OY);
      snap_found = True;
    }
  }
  else if (dc == (R1 + R2)) {		/* single point tangency */
    double theta = atan2(Y2, X2);
    double ix = R1 * cos(theta);
    double iy = R1 * sin(theta);
    
    if (isect_cb) insert_isect(isect_cb, ix, iy, -1);
    snap_gx = (int)lrint(ix);
    snap_gy = (int)lrint(iy);
    snap_found = True;
  }
  else {
    /* no intersection */
  }
}

static void
circle_circle_intersect(e1, e2, x, y, isect_cb)
     F_ellipse * e1;
     F_ellipse * e2;
     int x, y;
     isect_cb_s * isect_cb;
{
  double PX = (double)(x - e1->center.x);
  double PY = (double)(y - e1->center.y);
  double X2 = (double)(e2->center.x - e1->center.x);
  double Y2 = (double)(e2->center.y - e1->center.y);
  double R1 = (double)(e1->radiuses.x);
  double R2 = (double)(e2->radiuses.x);
  double OX = (double)(e1->center.x);
  double OY = (double)(e1->center.y);
  
  do_circle_circle(PX, PY, X2, Y2, R1, R2, OX, OY, NULL, NULL, isect_cb);
}



static void
non_ortho_ellipse_ellipse_intersect(e1, e2, x, y, isect_cb)
     F_ellipse * e1;
     F_ellipse * e2;
     int x, y;
     isect_cb_s * isect_cb;
{						/* ellipse-ellipse -- non-orthogonal */
  
  double A, B, C, D;
  double SIN, COS;
  double KV, KW;
  double E, F, G, H, J, K, L, M, N, P, Q, R, S, T, U, V, W;
  double  Z0,  Z1,  Z2,  Z3,  Z4, Z5, Z6, Z7, Z8, Z9;
  double Z10, Z11, Z12, Z13, Z14;
  double * poly;
  double * solr;
  double * soli;
  double ix, iye1, iye2;
  double PX = (double)(x - e1->center.x);
  double PY = (double)(y - e1->center.y);
  double a  = (double)(e1->radiuses.x);
  double b  = (double)(e1->radiuses.y);
  double c  = (double)(e2->radiuses.x);
  double d  = (double)(e2->radiuses.y);
  double X  = (double)(e2->center.x - e1->center.x);
  double Y  = (double)(e2->center.y - e1->center.y);
  double mind = HUGE_VAL;
  double dist;
  double hix, hiy;
  int nsol, i, ix_idx, iye1_idx, iye2_idx;
  
  snap_rotate_vector (&PX, &PY, PX, PY, (double)(e1->angle));
  snap_rotate_vector (&X,  &Y,  X,  Y,  (double)(e1->angle));
  
  poly = alloca(5 * sizeof(double));
  solr = alloca(4 * sizeof(double));
  soli = alloca(4 * sizeof(double));
  
  /*
   * x^2 / a^2   +   y^2 / b^2  =  1			Eq  1,  origin ctrd, ortho, ellipse 1
   *
   * A = a^2
   * B = b^2
   *
   */
  
  A = pow(a, 2.0);
  B = pow(b, 2.0);

  /*
   * x^2 / A  +  y^2 / B  = 1				Eq  2, from Eq 1
   *
   *   B x^2  +  A y^2    =  AB				Eq  3, from Eq 2
   *
   *             A y^2    =  AB  -  B x^2		Eq  4, from Eq 3
   *
   *               y^2    =   B  -  (B/A) x^2		Eq  5, from Eq 4
   *
   *               y      = sqrt(B  -  (B/A) x^2)	Eq  6, from Eq 5
   *
   */
  
  /*
   * v^2 / c^2   +   w^2 / d^2  =  1			Eq  7,  [X Y] origin, non-ortho, ellipse 2
   *
   * C = c^2
   * D = d^2
   *
   */
  
  C = pow(c, 2.0);
  D = pow(d, 2.0);
  
  /*	
   * v^2 / C  +  w^2 / D  =  1				Eq  8, from Eq 7
   *
   *  D v^2   +  C w^2    =  CD				Eq  9, from Eq 8
   *
   * [v w] = ([x y] - [X Y]) rot theta			[v w] is a rotation of [x y] around
   *							the center of ellipse 2
   *
   * v = (x - X) cos T + (y - Y) sin T			Eq 10
   * w = (x - X) sin T - (y - Y) cos T			Eq 11
   *
   * COS = cos T
   * SIN = sin T
   *
   */

  COS = cos((double)(e1->angle - e2->angle));
  SIN = sin((double)(e1->angle - e2->angle));

  /*
   * v = x COS - X COS  +  y SIN - Y SIN		Eq 12, from Eq 10
   * w = x SIN - X SIN  -  y COS + Y COS		Eq 13, from Eq 11
   *
   * v = (x COS + y SIN)  -   (X COS + Y SIN)		Eq 14, from Eq 12
   * w = (x SIN - y COS)  -   (X SIN - Y COS)		Eq 15, from Eq 13
   *
   * KV = X COS + Y SIN
   * KW = X SIN - Y COS
   *
   */
  
  KV = (X * COS) + (Y * SIN);
  KW = (X * SIN) - (Y * COS);

  
  /*
   * v = (x COS + y SIN) -  KV				Eq 16, from Eq 14
   * w = (x SIN - y COS) -  KW				Eq 17, from Eq 15
   *
   *
   *   D (x COS + y SIN -  KV)^2   +  C  (x SIN - y COS -  KW)^2  =  CD		Eq 18, from Eq  9
   *
   *   D ( COS^2 x^2 + COS SIN xy - COS KV x					Eq 19, from Eq 18
   *                 + COS SIN xy             + SIN^2 y^2 - SIN KV y
   *                              - COS KV x              - SIN KV y    + KV^2)
   * + C ( SIN^2 x^2 - COS SIN xy - SIN KW x
   *                 - COS SIN xy             + COS^2 y^2 + COS KW y
   *                              - SIN KW x              + COS KW y    + KW^2)  = CD
   *
   *   + (   D COS^2    +    C SIN^2   ) x^2					Eq 21, from Eq 20
   *   - ( 2 D COS KV   +  2 C SIN KW  )  x
   *   + ( 2 D COS SIN  -  2 C COS SIN ) xy
   *   - ( 2 D SIN KV   -  2 C COS KW  )  y
   *   + (   D SIN^2    +    C COS^2   ) y^2
   *   + (   D KV^2     +    C KW^2  - CD )  = 0
   *
   */
  
  E =  (      D * pow(COS, 2.0)    +        C * pow(SIN, 2.0) );
  F = -(2.0 * D * COS * KV         +  2.0 * C * SIN * KW      );
  G =  (2.0 * D * COS * SIN        -  2.0 * C * COS * SIN     );
  H = -(2.0 * D * SIN * KV         -  2.0 * C * COS * KW      );
  K =  (      D * pow(SIN, 2.0)    +        C * pow(COS, 2.0) );
  L =  (    ( D * pow(KV, 2.0)     +        C * pow(KW, 2.0) ) - C * D);

   
  /*
   *
   *   E x^2  +  F x  +  G xy  +  H y  +  K y^2  +  L  = 0			Eq 22, from Eq 21
   *
   *
   *     E x^2
   *  +  F x
   *  +  K (B - (B/A) x^2)
   *  +  L
   *             =
   *  -  G x sqrt(B - (B/A) x^2)				Eq 23, from Eq 22
   *  -  H sqrt(B - (B/A) x^2)
   *
   *
   *
   *     E x^2
   *  +  F x
   *  +  K B - K(B/A) x^2
   *  +  L
   *             =
   *  -  (G x + H) sqrt(B - (B/A) x^2)				Eq 23, from Eq 22
   *
   *
   *
   *     (E - KB/A) x^2
   *  +  F x
   *  +  (KB +  L)
   *             =
   *  -  (G x + H) sqrt(B - (B/A) x^2)				Eq 23, from Eq 22
   *
   *
   */

  M = E - (K * B / A);
  N = (K * B) + L;

  /*
   *    M x^2  +  F x  +  N  =  -(G x + H) sqrt(B - (B/A) x^2)			Eq 27, from Eq 26
   *
   *   (M x^2  +  F x  +  N)^2  =  (G x + H)^2 (B - (B/A) x^2)			Eq 28, from Eq 27
   *
   *         M^2 x^4    +     FM x^3         +    MN x^2			Eq 29, from Eq 28
   *                    +     FM x^3         +   F^2 x^2  +   FN x          
   *                                         +    MN x^2  +   FN x  + N^2
   *		   = (G^2 x^2  + 2GH x  + H^2)  (B - (B/A) x^2)
   *
   *         M^2 x^4  +      2FM x^3  +  (2MN + F^2) x^2  +  2FN x  + N^2	Eq 30, from Eq 29
   *                                        -  B G^2 x^2  - 2BGH x  - B H^2
   * + (B/A) G^2 x^4  + 2(B/A)GH x^3  +    (B/A) H^2 x^2  =  0
   *
   */

  poly[4] = pow(M, 2.0)  +  (B/A) * pow(G, 2.0);
  poly[3] = 2.0 * F * M  +  2.0 * (B/A) * G * H;
  poly[2] = (2.0 * M * N  +  pow(F, 2.0)  +  (B/A) * pow(H, 2.0)) - B * pow(G, 2.0);
  poly[1] = 2.0 * F * N  -  2.0 * B * G * H;
  poly[0] = pow(N, 2.0)  -  B * pow(H, 2.0);

  nsol = quartic(poly, solr, soli);
  
  for (i = 0; i < nsol; i++) {
    if (1.0 > fabs(soli[i])) {
      for (ix_idx = 0; ix_idx < 2; ix_idx++) {
	double rx;
	ix = (ix_idx & 1) ? -solr[i] : solr[i];
	rx = A - pow(ix, 2.0);						/* for origin ellipse 1 */
	if (0.0 <= rx) {
	  for (iye1_idx = 0; iye1_idx < 2; iye1_idx++) {
	    double v, w, res;
	    iye1 = (b/a) * ((iye1_idx & 1) ? -sqrt(rx) : sqrt(rx));		/* e1 */
	    
	    v = (ix - X) * COS + (iye1 - Y) * SIN;
	    w = (ix - X) * SIN - (iye1 - Y) * COS;
	    res = (pow(v, 2.0) / C)  +  (pow(w, 2.0) / D);

	    if (fabs(1.0 - res) < 0.0001) {
	      if (isect_cb) {
		double isx, isy;
		snap_rotate_vector (&isx,  &isy,  ix,  iye1,  -(double)(e1->angle));
		isx += (double)(e1->center.x);
		isy += (double)(e1->center.y);
		insert_isect(isect_cb, isx, isy, -1);
	      }
	      dist = hypot(iye1 - PY, ix - PX);
	      if (dist < mind) {
		mind = dist;
		hix = ix;
		hiy = iye1;
		snap_found = True;
	      }
	    }
	  }
	}
      }
    }
  }
  if (True == snap_found) {
    snap_rotate_vector (&hix, &hiy, hix, hiy, -(double)(e1->angle));
    snap_gx = (int)lrint(hix) + e1->center.x;
    snap_gy = (int)lrint(hiy) + e1->center.y;
  }
}


static void
ortho_ellipse_ellipse_intersect(e1, e2, x, y, isect_cb)
     F_ellipse * e1;
     F_ellipse * e2;
     int x, y;
     isect_cb_s * isect_cb;
{						/* ellipse-ellipse -- orthogonal */
  /*
   * x^2 / a^2   +   y^2 / b2  =  1					Eq  1,  origin ctrd, eclipse 1
   *
   * (x - X)^2 / c^2  +  (y - Y)^2 / d^2  =  1				Eq  2.  [X Y] ctrd eclipse 2
   *
   * A = a^2
   * B = b^2
   * C = c^2
   * D = d^2
   *
   * x^2 / A  +  y^2 / B  = 1						Eq  3, from Eq 1
   *	
   * (x - X)^2 / C  +  (y - Y)^2 / D  =  1				Eq  4, from Eq 2
   *
   * D(x - X)^2 / C  +  (y - Y)^2  =  D					Eq  4a, from Eq 4
   *
   *                    (y - Y)^2  =  D - D(x - X)^2 / C			Eq  4b, from Eq 4a
   *
   *                    (y - Y)^2  =  D(1 - (x - X)^2 / C)		Eq  4c, from Eq 4b
   *
   *		      (y - Y)^2  =  (D/C)(C - (x - X)^2)		Eq  4d, from Eq 4c
   *
   *		      (y - Y)    =  (d/c)sqrt(C - (x - X)^2)		Eq  4e, from Eq 4d
   *
   *		      y  =  Y  + (d/c)sqrt(C - (x - X)^2)		Eq  4f, from Eq 4e
   * 
   * B x^2  +  A y^2  =  AB						Eq  5, from Eq 3
   * 
   *           A y^2  =  AB  -  B x^2					Eq  6, from Eq 5
   *
   *             y^2  =  B   -  (B/A) x^2	 =  B(A - x^2)/A		Eq  7, from Eq 6
   *
   *	       y    =  (b/a)sqrt(A - x^2)				Eq  8, from Eq 7b
   *
   * D(x - X)^2  +  C(y - Y)^2  = CD					Eq  9, from Eq 4
   *
   * D(x^2 - 2Xx + X^2)  +  C(y^2 - 2Yy + Y^2)  =  CD			Eq 10, from Eq 9
   *
   * Dx^2 - 2DXx + DX^2  +  Cy^2 - 2CYy + CY^2  =  CD			Eq 11, from Eq 10
   *
   *	plug Eq 7a and Eq 8 into Eq 11
   *
   * Dx^2 - 2DXx + DX^2  +  C(B - (B/A) x^2) - 2CY((b/a)sqrt(A - x^2)) + CY^2  =  CD	Eq 12, from Eq 11
   *
   * Dx^2 - 2DXx + DX^2  +  CB - C(B/A) x^2 - 2CY((b/a)sqrt(A - x^2)) + CY^2  =  CD	Eq 13, from Eq 12
   *
   * (D - CB/A)x^2 - 2DXx + (DX^2 + CY^2 +  CB - CD)  =   2CY((b/a)sqrt(A - x^2))	Eq 14, from Eq 13
   *
   * E = D - CB/A
   * F = 2DX
   * G = DX^2 + CY^2 +  CB - CD
   * H = 2CY((b/a)
   *
   * Ex^2 - Fx + G  =   H sqrt(A - x^2)						Eq 15, from Eq 14
   *
   * (Ex^2 - Fx + G)^2  =   H^2 (A - x^2)						Eq 16, from Eq 15
   *
   * E^2 x^4  -  EF x^3  +   EG x^2
   *          -  EF x^3  +  F^2 x^2  -  FG x
   *                     +   EG x^2  -  FG x  +  G^2  = H^2 A  -  H^2 x^2		Eq 17, from Eq 16
   *
   * E^2 x^4 - 2EF x^3 + (2EG + F^2) x^2 - 2FG x + G^2 = H^2 A - H^2 x^2		Eq 18, from Eq 17
   *
   * E^2 x^4 - 2EF x^3 + (2EG + F^2 + H^2) x^2 - 2FG x + (G^2 - H^2 A) = 0	Eq 19, from Eq 18
   *
   * J = E^2
   * K = 2EF
   * L = 2EG + F^2 + H^2
   * M = 2FG
   * N = G^2 - H^2 A
   *
   * J x^4 - K x^3 + L x^2 - M x + N = 0					Eq 20, from Eq 19
   *
   */
  
  double A, B, C, D;
  double E, F, G, H, J, K, L, M, N;
  double * poly;
  double * solr;
  double * soli;
  double ix, iye1, iye2;
  double PX = (double)(x - e1->center.x);
  double PY = (double)(y - e1->center.y);
  double a = (double)(e1->radiuses.x);
  double b = (double)(e1->radiuses.y);
  double c = (double)(e2->radiuses.x);
  double d = (double)(e2->radiuses.y);
  double X = (double)(e2->center.x - e1->center.x);
  double Y = (double)(e2->center.y - e1->center.y);
  double mind = HUGE_VAL;
  double dist;
  double hix, hiy;
  int nsol, i, ix_idx, iye1_idx, iye2_idx;
  
  snap_rotate_vector (&PX, &PY, PX, PY, (double)(e1->angle));
  snap_rotate_vector (&X,  &Y,  X,  Y,  (double)(e1->angle));
  
  poly = alloca(5 * sizeof(double));
  solr = alloca(4 * sizeof(double));
  soli = alloca(4 * sizeof(double));
  
  A = pow(a, 2.0);
  B = pow(b, 2.0);
  C = pow(c, 2.0);
  D = pow(d, 2.0);
  
  E = D - (C * B)/A;
  F = 2.0 * D * X;
  G = ((D * pow(X, 2.0)) + (C * pow(Y, 2.0)) +  (C * B)) - (C * D);
  H = 2.0 * C * Y * (b/a);
  
  J = pow(E, 2.0);
  K = 2.0 * E * F;
  L = (2.0 * E * G) + pow(F, 2.0) + pow(H, 2.0);
  M = 2.0 * F * G;
  N = pow(G, 2.0) - (pow(H, 2.0) * A);

  poly[4] =  J;
  poly[3] = -K;
  poly[2] =  L;
  poly[1] = -M;
  poly[0] =  N;

  nsol = quartic(poly, solr, soli);
  
  for (i = 0; i < nsol; i++) {
    if (1.0 > fabs(soli[i])) {
      for (ix_idx = 0; ix_idx < 2; ix_idx++) {
	double rx, ex;
	ix = (ix_idx & 1) ? -solr[i] : solr[i];
	rx = A - pow(ix, 2.0);						/* for origin ellipse 1 */
	if (0.0 <= rx) {
	  for (iye1_idx = 0; iye1_idx < 2; iye1_idx++) {
	    iye1 = (b/a) * ((iye1_idx & 1) ? -sqrt(rx) : sqrt(rx));		/* e1 */
	    
	    /* y  =  Y  + (d/c)sqrt(C - (x - X)^2) */
	    
	    ex = C - pow(ix - X, 2.0);					/* for [X Y] ellipse 2 */
	    if (0.0 <= ex) {
	      for (iye2_idx = 0; iye2_idx < 2; iye2_idx++) {
		iye2 = Y + ((d/c) * ((iye2_idx & 1) ? -sqrt(ex) : sqrt(ex)));
		if (fabs(iye1 - iye2) < 1.0) {
		  if (isect_cb) {
		    double isx, isy;
		    snap_rotate_vector (&isx,  &isy,  ix,  iye1,  -(double)(e1->angle));
		    isx += (double)(e1->center.x);
		    isy += (double)(e1->center.y);
		    insert_isect(isect_cb, isx, isy, -1);
		  }
		  dist = hypot(iye1 - PY, ix - PX);
		  if (dist < mind) {
		    mind = dist;
		    hix = ix;
		    hiy = iye1;
		    snap_found = True;
		  }
		}
	      }
	    }
	  }
	}
      }
    }
  }
  if (True == snap_found) {
    snap_rotate_vector (&hix, &hiy, hix, hiy, -(double)(e1->angle));
    snap_gx = (int)lrint(hix) + e1->center.x;
    snap_gy = (int)lrint(hiy) + e1->center.y;
  }
}

void
intersect_ellipse_ellipse_handler(F_ellipse *  e1, F_ellipse *  e2, int x, int y, isect_cb_s * isect_cb )
{
  switch(e1->type) {
  case T_ELLIPSE_BY_RAD:
  case T_ELLIPSE_BY_DIA:
    switch(e2->type) {
    case T_ELLIPSE_BY_RAD:
    case T_ELLIPSE_BY_DIA:
      if (e1->angle == e2->angle) ortho_ellipse_ellipse_intersect(e1, e2, x, y, isect_cb);
      else non_ortho_ellipse_ellipse_intersect(e1, e2, x, y, isect_cb);
      break;
    case T_CIRCLE_BY_RAD:
    case T_CIRCLE_BY_DIA:
      circle_ellipse_intersect(e2, e1, x, y, isect_cb);
      break;
    }
    break;
  case T_CIRCLE_BY_RAD:
  case T_CIRCLE_BY_DIA:
    switch(e2->type) {
    case T_ELLIPSE_BY_RAD:
    case T_ELLIPSE_BY_DIA:
      circle_ellipse_intersect(e1, e2, x, y, isect_cb);
      break;
    case T_CIRCLE_BY_RAD:
    case T_CIRCLE_BY_DIA:
      circle_circle_intersect(e1, e2, x, y, isect_cb);
      break;
    }
    break;
  }
}

static void
do_intersect_ellipse_polyline(ecx, ecy, ea, eb, theta, l, x, y, arc, isect_cb)
     double ecx;
     double ecy;
     double ea;
     double eb;
     double theta;
     F_line * l;
     int x, y;
     F_arc * arc;
     isect_cb_s * isect_cb;
{
  /* if arc is non-null, ecx, ecy, ea, and eb arc from arc.*/
  /* if isect_cb is non-null, return the intersects */
  struct f_point * p;
  struct f_point * p_start;
  double dist;
  double mind;
  double c[3];
  double ix1, iy1;
  double ix2, iy2;
  int seg_idx;

  /* fixme -- handle ellipse v. circle */
  
  mind = HUGE_VAL;
  p_start = NULL;
  for (seg_idx = -1, p = l->points; p != NULL; seg_idx++, p = p->next) {
    if (p_start) {
      double rx;
      /* translate to ellipse origin */
      double psx = (double)(p_start->x) - ecx;
      double psy = (double)(p_start->y) - ecy;
      double pex = (double)(p->x) - ecx;
      double pey = (double)(p->y) - ecy;
      
      /* rotate around ellipse angle so ellipse semi-axes are ortho to space */
      snap_rotate_vector (&psx, &psy, psx, psy, theta);
      snap_rotate_vector (&pex, &pey, pex, pey, theta);
      
      if  (1) {					/* fixme find a quick screening test */
	int iset = 0;
	double M = pow(eb, 2.0);		
	double N = pow(ea, 2.0);
	fget_line_from_points(c, psx, psy, pex, pey);

	if ((0.0 == c[0]) && (0.0 != c[1]) && (0.0 != M)) {	/* a == 0; horzontal */
	  iy1 = iy2 = -c[2]/c[1];
	  rx = N - ((N/M) * pow(iy1, 2.0));
	  if (0.0 <= rx) {			/* tangent or normal two-point intersection */
	    ix1 =  sqrt(rx);
	    ix2 = -sqrt(rx);
	    iset = ISET_P1 | ISET_P2;
	  }
	  else if (-1.0 < rx) {			/* near miss sometimes caused by roundoff */
	    ix1 = ix2 = 0.0;
	    iset = ISET_P1;
	  }
	}
	else if ((0.0 != c[0]) && (0.0 == c[1]) && (0.0 != N)) {	/* b == 0; vertical */
	  ix1 = ix2 = -c[2]/c[0];
	  rx = M - ((M/N) * pow(ix1, 2.0));
	  if (0.0 <= rx) {			/* tangent or normal two-point intersection */
	    iy1 =  sqrt(rx);
	    iy2 = -sqrt(rx);
	    iset = ISET_P1 | ISET_P2;
	  }
	  else if (-1.0 < rx) {			/* near miss sometimes caused by roundoff */
	    iy1 = iy2 = 0.0;
	    iset = ISET_P1;
	  }
	}
	else if ((0.0 != c[0]) && (0.0 != c[1])) {
	  /* simplify things... */
	  /* normalise line coeffs... */
	  double K = -c[1]/c[0];
	  double L = -c[2]/c[0];
	  
	  double P = M * N;
	  
	  /* plug line into ellipse... */
	  double W = (M * pow(K, 2.0)) + N;
	  double X = 2.0 * K * L * M;
	  double Y = (M * pow(L, 2.0)) - P;
	  if (0.0 != W) {
	    rx = pow(X, 2.0) - (4.0 * W * Y);
	    if (0.0 <= rx) {
	      iy1 = ((-X) + sqrt(rx))/(2.0 * W);
	      ix1 = (K * iy1) + L;
	      iy2 = ((-X) - sqrt(rx))/(2.0 * W);
	      ix2 = (K * iy2) + L;
	      iset = ISET_P1 | ISET_P2;
	    }
	    else if (fabs(rx) < 10.0) {	/* fixme -- (too?) small arbitrary tolerance */
	      /* due to roundoff errors, the intersection of tangents to ellipses and
	       * circles may be slightly complex.  this code catches such near misses.
	       */
#if defined(HAVE_NO_COMPLEX)
              double dd;
	      double x1mag, y1mag;
	      double x2mag, y2mag;

	      iy1 = ((-X) )/(2.0 * W);
	      ix1 = (K * iy1) + L;
	      iy2 = ((-X) )/(2.0 * W);
	      ix2 = (K * iy2) + L;

	      y1mag = hypot(iy1, dd=(   (sqrt(-rx)))/(2.0 * W) );
	      x1mag = hypot(ix1, (K * dd) + L );

	      y2mag = hypot(iy2, dd=( - (sqrt(-rx)))/(2.0 * W) );
	      x2mag = hypot(ix2, (K * dd) + L );
	      
#else       
	      double complex ix1c;
	      double complex iy1c;
	      double complex ix2c;
	      double complex iy2c;
	      double x1mag, y1mag;
	      double x2mag, y2mag;

	      iy1c = ((-X) + (sqrt(-rx) *  _Complex_I))/(2.0 * W);
	      ix1c = (K * iy1c) + L;
	      
	      iy2c = ((-X) - (sqrt(-rx) *  _Complex_I))/(2.0 * W);
	      ix2c = (K * iy2c) + L;

	      ix1 = creal(ix1c);
	      iy1 = creal(iy1c);
	      ix2 = creal(ix2c);
	      iy2 = creal(iy2c);

	      x1mag = hypot(ix1, cimag(ix1c));
	      y1mag = hypot(iy1, cimag(iy1c));

	      x2mag = hypot(ix2, cimag(ix2c));
	      y2mag = hypot(iy2, cimag(iy2c));
#endif	    
	      
	      if ((1.0 > fabs(ix1) - x1mag) &&
		  (1.0 > fabs(iy1) - y1mag) &&
		  (1.0 > fabs(ix2) - x2mag) &&
		  (1.0 > fabs(iy2) - y2mag)) iset = ISET_P1 | ISET_P2;
	    }
	  }
	}

	if (0 != iset) {
	  Boolean point_ok;
	
	  if ((iset & ISET_P1) && (!is_point_on_segment(psx, psy, ix1, iy1, pex, pey)))
	    iset &= ~ISET_P1;
	  if ((iset & ISET_P2) && (!is_point_on_segment(psx, psy, ix2, iy2, pex, pey)))
	    iset &= ~ISET_P2;
	  
	  /* rotate back to screen space */
	  snap_rotate_vector (&ix1, &iy1, ix1, iy1, -theta);
	  snap_rotate_vector (&ix2, &iy2, ix2, iy2, -theta);
	      
	  /* translate back to screen space */
	  ix1 +=  ecx;
	  iy1 +=  ecy;
	  ix2 +=  ecx;
	  iy2 +=  ecy;

	  if (iset & ISET_P1) {
	    point_ok = (NULL == arc)
	      ? True
	      : is_point_on_arc(arc, (int)lrint(ix1), (int)lrint(iy1));
	    if (True == point_ok) {
	      if (isect_cb) insert_isect(isect_cb, ix1, iy1, seg_idx);
	      dist = hypot(iy1 - (double)y, ix1 - (double)x);
	      if (dist < mind) {
		mind = dist;
		snap_gx = (int)lrint(ix1);
		snap_gy = (int)lrint(iy1);
		snap_found = True;
	      }
	    }
	  }
	  
	  if (iset & ISET_P2) {
	    point_ok = (NULL == arc)
	      ? True
	      :	 is_point_on_arc(arc, (int)lrint(ix2), (int)lrint(iy2));
	    if (True == point_ok) {
	      if (isect_cb) insert_isect(isect_cb, ix2, iy2, seg_idx);
	      dist = hypot(iy2 - (double)y, ix2 - (double)x);
	      if (dist < mind) {
		mind = dist;
		snap_gx = (int)lrint(ix2);
		snap_gy = (int)lrint(iy2);
		snap_found = True;
	      }
	    }
	  }
	  
	  if (True == snap_found) {
	    /* check if we're near a polyline segment endpoint.  if so, snap to it. */
	    double de1 = hypot((double)(p_start->y - snap_gy),
			       (double)(p_start->x - snap_gx));
	    double de2 = hypot((double)(p->y - snap_gy),
				 (double)(p->x - snap_gx));
	    if (5.0 > de1) {			/* arbitarary tolerance... */
	      snap_gx = p_start->x;
	      snap_gy = p_start->y;
	    }
	    else if (5.0 > de2) {		/* arbitarary tolerance... */
	      snap_gx = p->x;
	      snap_gy = p->y;
	    }
	    snap_found = (NULL == arc) ? True : is_point_on_arc(arc, snap_gx, snap_gy);
	  }
	}
      }
    }
    p_start = p;
  }
  if (!isect_cb) {
    if (False == snap_found) {
      put_msg("Selected ellipse and polyline do not intersect.");
      beep();
      snap_msg_set = True;
    }
  }
}

void
intersect_ellipse_polyline_handler(F_ellipse * e, F_line *  l, int x, int y, isect_cb_s * isect_cb)
{
  double ecx = (double)(e->center.x);
  double ecy = (double)(e->center.y);
  double ea  = (double)(e->radiuses.x);
  double eb  = (double)(e->radiuses.y);
  double theta = (double)(e->angle);

  do_intersect_ellipse_polyline(ecx, ecy, ea, eb, theta, l, x, y, NULL, isect_cb);
}

static void
intersect_ellipse_spline_handler(obj1, obj2, x, y)
     void * obj1;
     void * obj2;
     int x, y;
{
  put_msg("Ellipse-spline intersections not yet implemented");
  beep();
  snap_msg_set = True;
}

void
delete_text_bounding_box(F_line * l)
{
  if (l) {
    if (l->points) free(l->points);
    free(l);
  }
}

F_line *
build_text_bounding_box(F_text * t) 
{
  /*
   * do everything wrt to a (rotated) square bounding box around the text.
   * the lower corners are at the descent, the upper at the ascent.
   *
   */
  struct f_point * points;
  int i;
  F_line * f_line_p = malloc(sizeof(F_line));
  f_line_p->type = T_BOX;
  f_line_p->points = points = malloc(5 * sizeof(struct f_point));

  points[0].x = 0;		points[0].y =   t->descent;
  points[1].x = t->length;	points[1].y =   t->descent;
  points[2].x = t->length;	points[2].y = -(t->ascent);
  points[3].x = 0;		points[3].y = -(t->ascent);
  points[4] = points[0];

  for (i = 0; i < 5; i++) {
    double tx, ty;
    switch(t->type) {
    case T_LEFT_JUSTIFIED:
      /* do nothing */
      break;
    case T_CENTER_JUSTIFIED:
      points[i].x -= (t->length)/2;
      break;
    case T_RIGHT_JUSTIFIED:
      points[i].x -= t->length;
      break;
    }
    tx = (double)points[i].x;
    ty = (double)points[i].y;
    snap_rotate_vector(&tx, &ty, tx, ty, -(double)(t->angle));
    points[i].x = t->base_x + (int)lrint(tx);
    points[i].y = t->base_y + (int)lrint(ty);
    points[i].next = (i < 4) ? &points[i+1] : NULL;
  }
  return f_line_p;
}

static void
intersect_ellipse_text_handler(e, t, x, y)
     F_ellipse * e;
     F_text * t;
     int x, y;
{
  F_line * f_line_p = build_text_bounding_box(t);
  intersect_ellipse_polyline_handler(e, f_line_p, x, y, NULL);
  delete_text_bounding_box(f_line_p);
}

void
intersect_ellipse_arc_handler(F_ellipse * e, F_arc * a, int x, int y, isect_cb_s * isect_cb)
{
  switch(e->type) {
  case T_ELLIPSE_BY_RAD:
  case T_ELLIPSE_BY_DIA:
    {
      double r = hypot((double)(a->center.y) - (double)(a->point[1].y),
		       (double)(a->center.x) - (double)(a->point[1].x));
      double X  = (double)(a->center.x) - (double)(e->center.x);
      double Y  = (double)(a->center.y) - (double)(e->center.y);
 
      do_circle_ellipse_intersect(r, X, Y, e, x, y, a, isect_cb);
    }
    break;
  case T_CIRCLE_BY_RAD:
  case T_CIRCLE_BY_DIA:
    {
      double PX = (double)(x - e->center.x);
      double PY = (double)(y - e->center.y);
      double X2 = (double)(a->center.x) - (double)(e->center.x);
      double Y2 = (double)(a->center.y) - (double)(e->center.y);
      double R1 = (double)(e->radiuses.x);
      double R2 =  hypot((double)(a->center.y) - (double)(a->point[1].y),
			 (double)(a->center.x) - (double)(a->point[1].x));
      double OX = (double)(e->center.x);
      double OY = (double)(e->center.y);
  
      do_circle_circle(PX, PY, X2, Y2, R1, R2, OX, OY, a, isect_cb);
    }
  }
}

void
intersect_polyline_polyline_handler(F_line * l1, F_line * l2, int x, int y, isect_cb_s * isect_cb)
{
  struct f_point * p1;
  struct f_point * p2;
  struct f_point * p1_start;
  struct f_point * p2_start;
  double c1[3];
  double c2[3];
#define lA c1[0]
#define lB c1[1]
#define lC c1[2]
#define lD c2[0]
#define lE c2[1]
#define lF c2[2]
  double ci[2][2];
#define ciA ci[0][0]  
#define ciB ci[0][1]  
#define ciC ci[1][0]  
#define ciD ci[1][1]
  double det;
  double ix, iy;
  double dist;
  double mind;
  int seg_idx;

  mind = HUGE_VAL;
  p1_start = NULL;
  for (p1 = l1->points; p1 != NULL; p1 = p1->next) {
    if (p1_start) {
      get_line_from_points(c1, p1_start, p1);
      p2_start = NULL;
      for (seg_idx = -1, p2 = l2->points; p2 != NULL; seg_idx++, p2 = p2->next) {
	if (p2_start) {
	  if (True == boxes_overlap(p1_start, p1, p2_start, p2)) {
	    get_line_from_points(c2, p2_start, p2);
	    det = (lB * lD) - (lA * lE);
	    ciA = lE/det;
	    ciB = -lD/det;
	    ciC = -lB/det;
	    ciD = lA/det;
	    ix = (lC * ciA) + (lF * ciC);
	    iy = (lC * ciB) + (lF * ciD);
	    snap_found = True;
	    if (isect_cb) {
	      if ((True == is_point_on_segment((double)(p1_start->x), (double)(p1_start->y),
					       ix, iy,
					       (double)(p1->x), (double)(p1->y))) && 
		  (True == is_point_on_segment((double)(p2_start->x), (double)(p2_start->y),
					       ix, iy,
					       (double)(p2->x), (double)(p2->y))))
		insert_isect(isect_cb, ix, iy, seg_idx);
	    }
	    else {
	      dist = hypot(iy - (double)y, ix - (double)x);
	      if (dist < mind) {
		mind = dist;
		snap_gx = (int)lrint(ix);
		snap_gy = (int)lrint(iy);
	      }
	    }
	  }
	}
	p2_start = p2;
      }
    }
    p1_start = p1;
  }
  if (!isect_cb) {
    if (False == snap_found) {
      put_msg("Selected polylines do not intersect.");
      beep();
      snap_msg_set = True;
    }
  }
}

static void
intersect_polyline_spline_handler(obj1, obj2, x, y)
     void * obj1;
     void * obj2;
     int x, y;
{
  put_msg("Polyline-spline intersections not yet implemented");
  beep();
  snap_msg_set = True;
}

static void
intersect_polyline_text_handler(l, t, x, y)
     F_line * l;
     F_text * t;
     int x, y;
{
  F_line * f_line_p = build_text_bounding_box(t);
  intersect_polyline_polyline_handler(l, f_line_p, x, y, NULL);
  delete_text_bounding_box(f_line_p);
}

void
intersect_polyline_arc_handler(F_line * l,  F_arc *a,  int x, int y, isect_cb_s * isect_cb)
{
  double ecx = (double)(a->center.x);
  double ecy = (double)(a->center.y);
  double ea  =  hypot((double)(a->center.y) - (double)(a->point[1].y),
		      (double)(a->center.x) - (double)(a->point[1].x));
  double theta = 0.0;

  do_intersect_ellipse_polyline(ecx, ecy, ea, ea, theta, l, x, y, a, isect_cb);
}

static void
intersect_spline_spline_handler(obj1, obj2, x, y)
     void * obj1;
     void * obj2;
     int x, y;
{
  put_msg("Spline-spline intersections not yet implemented");
  beep();
  snap_msg_set = True;
}

static void
intersect_spline_text_handler(s, t, x, y)
     F_spline * s;
     F_text * t;
     int x, y;
{
  put_msg("Spline-text intersections not yet implemented");
  beep();
  snap_msg_set = True;
}

static void
intersect_spline_arc_handler(obj1, obj2, x, y)
     void * obj1;
     void * obj2;
     int x, y;
{
  put_msg("Spline-arc intersections not yet implemented");
  beep();
  snap_msg_set = True;
}

static void
intersect_text_text_handler(t1, t2, x, y)
     F_text * t1;
     F_text * t2;
     int x, y;
{
  F_line * f1_line_p = build_text_bounding_box(t1);
  F_line * f2_line_p = build_text_bounding_box(t2);
  intersect_polyline_polyline_handler(f1_line_p, f2_line_p, x, y, NULL);
  delete_text_bounding_box(f1_line_p);
  delete_text_bounding_box(f2_line_p);
}

static void
intersect_text_arc_handler(t, a, x, y)
     F_text * t;
     F_arc * a;
     int x, y;
{
  F_line * f_line_p = build_text_bounding_box(t);
  intersect_polyline_arc_handler(f_line_p, a, x, y, NULL);
  delete_text_bounding_box(f_line_p);
}

void
intersect_arc_arc_handler(F_arc * a1, F_arc * a2, int x, int y, isect_cb_s * isect_cb)
{
  /* translate to center of a1 */
  double PX = (double)(x) - (double)(a1->center.x);	/* translated event point */
  double PY = (double)(y) - (double)(a1->center.y);
  double X2 = (double)(a2->center.x - a1->center.x);	/* translated center of a2 */
  double Y2 = (double)(a2->center.y - a1->center.y);
  double R1 = hypot((double)(a1->center.y) - (double)(a1->point[1].y),
		    (double)(a1->center.x) - (double)(a1->point[1].x));
  double R2 = hypot((double)(a2->center.y) - (double)(a2->point[1].y),
		    (double)(a2->center.x) - (double)(a2->point[1].x));
  double OX = (double)(a1->center.x);
  double OY = (double)(a1->center.y);
  
  do_circle_circle(PX, PY, X2, Y2, R1, R2, OX, OY, a1, a2, isect_cb);
}

static void
intersect_ellipse_handler(obj1, obj2, type2, x, y)
     void * obj1;
     void * obj2;
     int type2;
     int x, y;
{
  switch (type2) {
  case O_ELLIPSE:
    intersect_ellipse_ellipse_handler(obj1, obj2, x, y, NULL);
    break;
  case O_POLYLINE:
    intersect_ellipse_polyline_handler(obj1, obj2, x, y, NULL);
    break;
  case O_SPLINE:
    intersect_ellipse_spline_handler(obj1, obj2, x, y);
    break;
  case O_TEXT:
    intersect_ellipse_text_handler(obj1, obj2, x, y);
    break;
  case O_ARC:
    intersect_ellipse_arc_handler(obj1, obj2, x, y, NULL);
    break;
  }
}

static void
intersect_polyline_handler(obj1, obj2, type2, x, y)
     void * obj1;
     void * obj2;
     int type2;
     int x, y;
{
  switch (type2) {
  case O_ELLIPSE:
    intersect_ellipse_polyline_handler(obj2, obj1, x, y, NULL);
    break;
  case O_POLYLINE:
    intersect_polyline_polyline_handler(obj1, obj2, x, y, NULL);
    break;
  case O_SPLINE:
    intersect_polyline_spline_handler(obj1, obj2, x, y);
    break;
  case O_TEXT:
    intersect_polyline_text_handler(obj1, obj2, x, y);
    break;
  case O_ARC:
    intersect_polyline_arc_handler(obj1, obj2, x, y, NULL);
    break;
  }
}

static void
intersect_spline_handler(obj1, obj2, type2, x, y)
     void * obj1;
     void * obj2;
     int type2;
     int x, y;
{
  switch (type2) {
  case O_ELLIPSE:
    intersect_ellipse_spline_handler(obj2, obj1, x, y);
    break;
  case O_POLYLINE:
    intersect_polyline_spline_handler(obj2, obj2, x, y);
    break;
  case O_SPLINE:
    intersect_spline_spline_handler(obj1, obj2, x, y);
    break;
  case O_TEXT:
    intersect_spline_text_handler(obj1, obj2, x, y);
    break;
  case O_ARC:
    intersect_spline_arc_handler(obj1, obj2, x, y);
    break;
  }
}

static void
intersect_text_handler(obj1, obj2, type2, x, y)
     void * obj1;
     void * obj2;
     int type2;
     int x, y;
{
  switch (type2) {
  case O_ELLIPSE:
    intersect_ellipse_text_handler(obj2, obj1, x, y);
    break;
  case O_POLYLINE:
    intersect_polyline_text_handler(obj2, obj2, x, y);
    break;
  case O_SPLINE:
    intersect_spline_text_handler(obj2, obj1, x, y);
    break;
  case O_TEXT:
    intersect_text_text_handler(obj1, obj2, x, y);
    break;
  case O_ARC:
    intersect_text_arc_handler(obj1, obj2, x, y);
    break;
  }
}

static void
intersect_arc_handler(obj1, obj2, type2, x, y)
     void * obj1;
     void * obj2;
     int type2;
     int x, y;
{
  switch (type2) {
  case O_ELLIPSE:
    intersect_ellipse_arc_handler(obj2, obj1, x, y, NULL);
    break;
  case O_POLYLINE:
    intersect_polyline_arc_handler(obj2, obj1, x, y, NULL);
    break;
  case O_SPLINE:
    intersect_spline_arc_handler(obj2, obj1, x, y);
    break;
  case O_TEXT:
    intersect_text_arc_handler(obj2, obj1, x, y);
    break;
  case O_ARC:
    intersect_arc_arc_handler(obj1, obj2, x, y, NULL);
    break;
  }
}

void
snap_intersect_handler(obj1, type1, obj2, type2, x, y)
     void * obj1;
     int type1;
     void * obj2;
     int type2;
     int x, y;
{
  switch (type1) {
  case O_ELLIPSE:
    intersect_ellipse_handler(obj1, obj2, type2, x, y);
    break;
  case O_POLYLINE:
    intersect_polyline_handler(obj1, obj2, type2, x, y);
    break;
  case O_SPLINE:
    intersect_spline_handler(obj1, obj2, type2, x, y);
    break;
  case O_TEXT:
    intersect_text_handler(obj1, obj2, type2, x, y);
    break;
  case O_ARC:
    intersect_arc_handler(obj1, obj2, type2, x, y);
    break;
  }
  if ((False == snap_found) && (False == snap_msg_set)) {
    put_msg("No intersection found.");
    beep();
    snap_msg_set = True;
  }
}