~jaspervdg/+junk/aem-diffusion-curves

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
/* circle_circle_intersection() *
 * Determine the points where 2 circles in a common plane intersect.
 *
 * int circle_circle_intersection(
 *                                // center and radius of 1st circle
 *                                double x0, double y0, double r0,
 *                                // center and radius of 2nd circle
 *                                double x1, double y1, double r1,
 *                                // 1st intersection point
 *                                double *xi, double *yi,              
 *                                // 2nd intersection point
 *                                double *xi_prime, double *yi_prime)
 *
 * This is a public domain work. 3/26/2005 Tim Voght
 * Ported to lib2geom, 2006 Nathan Hurst
 *
 * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
 *
 * This library is free software; you can redistribute it and/or
 * modify it either under the terms of the GNU Lesser General Public
 * License version 2.1 as published by the Free Software Foundation
 * (the "LGPL") or, at your option, under the terms of the Mozilla
 * Public License Version 1.1 (the "MPL"). If you do not alter this
 * notice, a recipient may use your version of this file under either
 * the MPL or the LGPL.
 *
 * You should have received a copy of the LGPL along with this library
 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 * You should have received a copy of the MPL along with this library
 * in the file COPYING-MPL-1.1
 *
 * The contents of this file are subject to the Mozilla Public License
 * Version 1.1 (the "License"); you may not use this file except in
 * compliance with the License. You may obtain a copy of the License at
 * http://www.mozilla.org/MPL/
 *
 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
 * the specific language governing rights and limitations.
 * 
 */
#include <stdio.h>
#include <math.h>
#include <2geom/point.h>

namespace Geom{

int circle_circle_intersection(Point X0, double r0,
                               Point X1, double r1,
                               Point & p0, Point & p1)
{
  /* dx and dy are the vertical and horizontal distances between
   * the circle centers.
   */
  Point D = X1 - X0;

  /* Determine the straight-line distance between the centers. */
  double d = L2(D);

  /* Check for solvability. */
  if (d > (r0 + r1))
  {
    /* no solution. circles do not intersect. */
    return 0;
  }
  if (d <= fabs(r0 - r1))
  {
    /* no solution. one circle is contained in the other */
    return 1;
  }
  
  /* 'point 2' is the point where the line through the circle
   * intersection points crosses the line between the circle
   * centers.  
   */

  /* Determine the distance from point 0 to point 2. */
  double a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;

  /* Determine the coordinates of point 2. */
  Point p2 = X0 + D * (a/d);

  /* Determine the distance from point 2 to either of the
   * intersection points.
   */
  double h = sqrt((r0*r0) - (a*a));

  /* Now determine the offsets of the intersection points from
   * point 2.
   */
  Point r = (h/d)*rot90(D);

  /* Determine the absolute intersection points. */
  p0 = p2 + r;
  p1 = p2 - r;

  return 2;
}

};


#ifdef TEST

void run_test(double x0, double y0, double r0,
              double x1, double y1, double r1)
{
  double x3, y3, x3_prime, y3_prime;

  printf("x0=%F, y0=%F, r0=%F, x1=%F, y1=%F, r1=%F :\n",
          x0, y0, r0, x1, y1, r1);
  Geom::Point p0, p1;
  Geom::circle_circle_intersection(Geom::Point(x0, y0), r0, 
				   Geom::Point(x1, y1), r1,
				   p0, p1);
  printf("  x3=%F, y3=%F, x3_prime=%F, y3_prime=%F\n",
            p0[0], p0[1], p1[0], p1[1]);
}

int main(void)
{
  /* Add more! */    
  run_test(-1.0, -1.0, 1.5, 1.0, 1.0, 2.0);
  run_test(1.0, -1.0, 1.5, -1.0, 1.0, 2.0);
  run_test(-1.0, 1.0, 1.5, 1.0, -1.0, 2.0);
  run_test(1.0, 1.0, 1.5, -1.0, -1.0, 2.0);
  exit(0);
}
#endif

/*
  Local Variables:
  mode:c++
  c-file-style:"stroustrup"
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
  indent-tabs-mode:nil
  fill-column:99
  End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :