Get the equation of a circle through the points $(1,1), (2,4), (5,3) $.
I can solve this by simply drawing it, but is there a way of solving it (easily) without having to draw?
Get the equation of a circle through the points $(1,1), (2,4), (5,3) $.
I can solve this by simply drawing it, but is there a way of solving it (easily) without having to draw?
Follow these steps:
Based on what I've learnt from this page by Stephen R. Schmitt (thanks to @Dan_Uznanski for pointing the concept out), there's a faster way to do it using matrices:
$$\text{let }〈x_1, y_1〉, 〈x_2, y_2〉, 〈x_3, y_3〉\text{ be your 3 points, and let }〈x_0, y_0〉\text{ represent the center of the circle.}\\\ \\ \text{let }A = \left[\begin{array}{cccc} x^2+y^2 & x & y & 1\\ x_1^2+y_1^2 & x_1 & y_1 & 1\\ x_2^2+y_2^2 & x_2 & y_2 & 1\\ x_3^2+y_3^2 & x_3 & y_3 & 1\\ \end{array}\right] = \left[\begin{array}{cccc} x^2+y^2 & x & y & 1\\ 1^2+1^2 & 1 & 1 & 1\\ 2^2+4^2 & 2 & 4 & 1\\ 5^2+3^2 & 5 & 3 & 1\\ \end{array}\right] = \left[\begin{array}{cccc} x^2+y^2 & x & y & 1\\ 2 & 1 & 1 & 1\\ 20 & 2 & 4 & 1\\ 34 & 5 & 3 & 1\\ \end{array}\right]\\\ \\\ \\\ \text{Note: }M_{11} = 0 ⟹ 〈x_0, y_0〉\text{ is undefined}∧\text{the points lie on a line rather than a circle.}\\\ \\\ x_0 = \frac{1}{2} ⋅ \frac{M_{12}}{M_{11}} = \frac{1}{2} ⋅ \left|\begin{array}{cccc} 2 & 1 & 1\\ 20 & 4 & 1\\ 34 & 3 & 1\\ \end{array}\right| \Bigg{/} \left|\begin{array}{cccc} 1 & 1 & 1\\ 2 & 4 & 1\\ 5 & 3 & 1\\ \end{array}\right| = \frac{1}{2} ⋅ \frac{-60}{-10} = 3 $$
$$ y_0 = -\frac{1}{2} ⋅ \frac{M_{13}}{M_{11}} = -\frac{1}{2} ⋅ \left|\begin{array}{cccc} 2 & 1 & 1\\ 20 & 2 & 1\\ 34 & 5 & 1\\ \end{array}\right| \Bigg{/} \left|\begin{array}{cccc} 1 & 1 & 1\\ 2 & 4 & 1\\ 5 & 3 & 1\\ \end{array}\right| = -\frac{1}{2} ⋅ \frac{40}{-10} = 2 $$
The radius, $r$, can then be calculated with Pythogoras' theorem, or using matrices again:
$$r^2 = x_0^2 + y_0^2 + \frac{M_{14}}{M_{11}}$$
If you have the time, it's really worth learning about matrices — they make a lot of things much faster, and are just generally awesome! KhanAcademy has a pretty easy introduction to matrices and linear algebra.
Coming back to this question because I was looking for ways to do this quickly on a computer without matrix math or systems of equations. A derivation follows, with the result being only 3 lines of simple code (plus error checking).
This problem becomes quite easy if we map each pair of $(x, y)$ points to three points in the complex plane of the form $z = x+iy$. We now apply the linear transformation $z \rightarrow \frac{z - z_1}{z_2 - z_1}$, which transforms the unique points $[z_1, z_2, z_3]$ to the points $[0, 1, w = \frac{z_3 - z_1}{z_2 - z_1}]$.
We now square the relation $|z - c| = r$ for each transformed point to find the following set of three equations:
$$|c|^2 = r^2$$ $$1 - c - \bar c +|c|^2 = r^2$$ $$|w|^2 - \bar w c - w\bar c +|c|^2 = r^2$$
Subtracting the first from the second and third equations gives us the series of equations: $$ \begin{bmatrix} 1 & 1 \\ \bar w & w \\ \end{bmatrix} \begin{bmatrix} c \\ \bar c \\ \end{bmatrix} = \begin{bmatrix} 1 \\ |w|^2 \\ \end{bmatrix} $$ Invert the left hand side and multiply to get: $$ \begin{bmatrix} c \\ \bar c \\ \end{bmatrix} = \frac{1}{w - \bar w} \begin{bmatrix} w - |w|^2 \\ -\bar w + |w|^2 \\ \end{bmatrix} $$
Note that this is singular if and only if $w = \bar w$, which is to say if $w$ is real, which means that the points $[0, 1, w]$ are collinear on the real axis. And since we applied a linear transformation, this only happens if the original points $[z_1, z_2, z_3]$ were collinear in the complex plane to begin with!
The center of the circle can now be read off from the first row:
$$c' = \frac{w - |w|^2}{w - \bar w}$$
But remember that we need to undo the linear transformation, so here is the answer in our original coordinates:
$$c = (z_2 - z_1)\frac{w - |w|^2}{w - \bar w} + z_1$$
The radius $r$ can be found from the equation before.
Or to make this nice and tidy, here is how you solve the problem in a few lines of python (no imports needed):
def circle_from_3_points(z1:complex, z2:complex, z3:complex) -> tuple[complex, float]:
if (z1 == z2) or (z2 == z3) or (z3 == z1):
raise ValueError(f"Duplicate points: {z1}, {z2}, {z3}")
w = (z3 - z1)/(z2 - z1)
# You should change 0 to a small tolerance for floating point comparisons
if abs(w.imag) <= 0:
raise ValueError(f"Points are collinear: {z1}, {z2}, {z3}")
c = (z2 - z1)*(w - abs(w)**2)/(2j*w.imag) + z1 # Simplified denominator
r = abs(z1 - c)
return c, r
>> print(circle_from_3_points(1+1j, 2+4j, 5+3j)) # Can also generate 1+1j with `complex(1,1)`
((3+2j), 2.2360679775)
>> print(circle_from_3_points(1+1j, 1+1j, 3+3j))
ValueError: Duplicate points: (1+1j), (1+1j), (3+3j)
>> print(circle_from_3_points(1+1j, 2+2j, 3+3j))
ValueError: Points are collinear: (1+1j), (2+2j), (3+3j)
Thanks to David Wright of Oklahoma State University for this method.
I'm surprised this one hasn't been mentioned; you can find the equation by using the determinant of a matrix:
$$\left|\begin{array}{cccc} x^2+y^2&x&y&1\\ 1^2+1^2&1&1&1\\ 2^2+4^2&2&4&1\\ 5^2+3^2&5&3&1\\ \end{array}\right|=0$$ This gives the equation of the circle through those three points. This sort of thing can be used in a lot of situations: matrix-determinant solutions are available for any shape I can think of where you're given points that land on the shape.
Implementing this on a computer involves having a thing that calculates determinants; to do it numerically you'll need to apply the cofactors method to avoid multiplying the variables in.
You can also find first $R$ from the sin Law:
$$R= \frac{BC}{2 \sin (A)}= \frac{BC \cdot AB \cdot AC }{2 \| AB \times AC \|} \tag{$*$}$$
Next, write the equations of circles of radius $R$ with centre $A$ and $B$ and solve.
Note The formula $(*)$ is the well known geometric formula for the area of a triangle:
$$\mbox{Area}= \frac{abc}{4R} \,.$$
As a programmer this was the best solution for me as there is no division by zero:
For given $(x_1,y_1) , (x_2,y_2)$ and $(x_3,y_3)$ first form $A,B,C$ and $D$:
$$A=x_1(y_2-y_3) - y_1(x_2-x_3)+x_2y_3 -x_3y_2$$ $$B=(x_1^2 + y_1^2)(y_3-y_2)+(x_2^2+y_2^2)(y_1-y_3)+(x_3^2+y_3^2)(y_2-y_1)$$ $$C=(x_1^2+y_1^2)(x_2-x_3)+(x_2^2+y_2^2)(x_3-x_1)+(x_3^2+y_3^2)(x_1-x_2)$$ $$D=({x_1^2}+{y_1^2})({x_3}{y_2} - {x_2}{y_3})+(x_2^2 + y_2^2)(x_1y_3-x_3y_1)+(x_3^2+y_3^2)(x_2y_1 - x_1y_2)$$
If A=0 then the points are colinear elasewhere using A,B,C you can find center and radius of circle:
$$x_c=-B/2A$$ $$y_c=-C/2A$$
$$R= \sqrt {\frac {B^2+C^2-4AD}{4A^2}}$$
Finally the formula of the circle is:
$$(x-x_c)^2+(y-y_c)^2=R^2$$
Big hint:
Let $A\equiv (1,1)$,$B\equiv (2,4)$ and $C\equiv (5,3)$.
We know that the perpendicular bisectors of the three sides of a triangle are concurrent.Join $A$ and $B$ and also $B$ and $C$.
The perpendicular bisector of $AB$ must pass through the point $(\frac{1+2}{2},\frac{1+4}{2})$
Now find the equations of the straight lines AB and BC and after that the equation of the perpendicular bisectors of $AB$ and $BC$.Solve for the equations of the perpendicular bisectors of $AB$ and $BC$ to get the centre of your circle.
“Plücker’s mu” provides another way to do this. Choose two of the points, $\mathbf p_1$ and $\mathbf p_2$, and let the equations of two circles that pass through those points be $f(x,y)=0$ and $g(x,y)=0$. Every linear combination of these two equations is also that of a circle (perhaps degenerate) that passes through these points. Then, an equation of the circle that also passes through the point $\mathbf p_0=(x_0,y_0)$ is $f(x_0,y_0)g(x,y)-g(x_0,y_0)f(x,y)=0$. This method applies generally to point sets that are given as null sets of functions.
Sometimes the coordinates are such that one can build the two required functions by inspection, but one can always use the “degenerate circle”—the line between two of the points—in this construction. For example, taking $\mathbf p_1=(1,1)$ and $\mathbf p_2=(5,3)$, we can use the circle centered at their midpoint: $$(x-3)^2+(y-2)^2=5$$ and the double line through these points: $$(x-2y+1)^2=0.$$ The required equation is then $$[(2-3)^2+(4-2)^2-5](x-2y+1)^2-(3-2\cdot4+1)^2[(x-3)^2+(y-2)^2-5]=0,$$ which simplifies to $(x-3)^2+(y-2)^2=5$. (This could also have been computed in matrix form, but that would’ve required expanding the equation of the line. In this case it was easier to work directly with the equations.)
Given 3 points $P_1(x_1,y_1) ,P_2(x_2,y_2)$ and $P_3(x_3,y_3)$ on circle circumference whose center is (h,k). Circle equation for each pair of points can be expressed as:
$$(x_1-h)^2+(y_1-k)^2 = (x_2-h)^2+(y_2-k)^2 ,$$ and $$(x_1-h)^2+(y_1-k)^2 = (x_3-h)^2+(y_3-k)^2 ,$$
By rearranging those equations, they can easily be put in the form of "difference of squares" in the form: $$A^2 - B^2 = C^2-D^2$$
we get : $$(x_1-h)^2 - (x_2-h)^2 = (y_2-k)^2 - (y_1-k)^2 ,$$ Where $$ A = (x_1-h) , B= (x_2-h) , C = (y_2-k) , D = (y_1-k)$$ Solving $$ A^2-B^2 = (A-B)(A+B) $$ we get: $$[ (x_1-h) - (x_2-h) ] [(x_1-h) + (x_2-h) ] = [ (y_2-k)-(y_1-k) ] * [ (y_2-k) + (y_1-k) ]$$
simplifying:
$$ (x_1-x_2)(x_1+x_2-2h) = (y_2-y_1)(y_2+y_1 -2k)$$
Introducing $$ a_1=(x_1-x_2) , a_1' =(x_1+x_2) , b_1 = (y_2-y_1) , b_1' =(y_2+y_1) $$ And $$ c_1 = \frac {(b_1b_1' - a_1a_1')} { 2a_1} $$
we get for each pair of points : $$ h = \frac {b_1}{a_1} k - c_1 $$ and another reduced equation for the other pair of points would give: $$ h = \frac {b_2} {a_2} k - c_2 $$
Solving above equations , we get : $$ k = \frac {a_1a_2(c_1-c_2)} { (b_1a_2- b_2a_1)} $$
and $$ r = \sqrt (E^2 + D^2) $$
Where $$ E=(x1-h) $$ and $$ D = (y1-k) $$
Consider the general (implicit) equation that defines a circle, with parameters $\alpha$, $\beta$, $\gamma$. Substitute the coordinate of the given points and get three linear equations in the three variables $\alpha$, $\beta$, $\gamma$. Solve the system.
I just wanted to flag Scott's answer to this question which is very concise and is probably faster than the matrix-based version. I've included below a C++17 implementation of the complex arithmetic solution to the problem
#include <complex>
#include <optional>
#include <iostream>
struct point {
double x;
double y;
};
struct circle {
double x;
double y;
double r;
};
std::optional<circle> circle_through_three_points(const point& pt1, const point& pt2, const point& pt3) {
using namespace std::complex_literals;
auto pt_to_z = [](const point& pt)->std::complex<double> { return { pt.x,pt.y }; };
auto z1 = pt_to_z(pt1), z2 = pt_to_z(pt2), z3 = pt_to_z(pt3);
auto w = (z3 - z1) / (z2 - z1);
if (w.imag() == 0)
return std::nullopt; // the three points are collinear
auto magnitude_w = std::abs(w);
auto c = (z2 - z1) * (w - magnitude_w * magnitude_w) / (2i * w.imag()) + z1;
auto r = std::abs(z1 - c);
return circle{ c.real(), c.imag(), r };
}
int main()
{
auto c1 = circle_through_three_points({ 1,1 }, { 2,4 }, { 5,3 });
auto c2 = circle_through_three_points({ 1,1 }, { 2,2 }, { 3,3 });
if (c1)
std::cout << "circle 1 is (" << c1->x << "," << c1->y << "," << c1->r << ")\n";
else
std::cout << "circle 1 is degenerate\n";
if (c2)
std::cout << "circle 2 is (" << c2->x << "," << c2->y << "," << c2->r << ")\n";
else
std::cout << "circle 2 is degenerate\n";
}
This solves your problem in python...
def three_point_circle(z1,z2,z3):
a = 1j*(z1-z2)
b = 1j*(z3-z2)
m1 = a.imag/a.real
c = (z1-z2)/2
p1 = z2+c
b1 = p1.imag-m1*p1.real
m2 = b.imag/b.real
d = (z3-z2)/2
p2 = z2+d
b2 = p2.imag-m2*p2.real
x = (b2-b1)/(m1-m2)
y = (m2*b1-m1*b2)/(m2-m1)
center = x+1j*y
radius = abs(center-z1)
return x,y,radius
as long as there is not division by zero.
If there is division by zero, you can solve this by revolving your input numbers (multiply by some random power of the imaginary unit 1j) and then revolve the x,y coordinates back (set to a complex number and multiply by 1j to the negative power you used).
It must be unlucky to add a 13th answer to an 8-year-old question that has been viewed 131,000 times (I don't know why it appeared on the front page today), but I just have to point out that this particular example can be solved by a simple mental calculation.
If $A = (1,1),$ $B = (2,4)$ and $C = (5,3),$ then - making $B$ the "origin", to see if it helps (I tried $A$ first, and would have tried $C$ next, if $B$ hadn't worked) - we have $A - B = (-1, -3)$ and $C - B = (3, - 1),$ whence $(A - B) \cdot (C - B) = -3 + 3 = 0,$ i.e., the angle at $B$ is a right angle, therefore $AC$ is a diameter, therefore the centre is $(A + C)/2 = (3, 2),$ and a radius vector is $(2, 1),$ therefore the square of the radius is $5,$ therefore the equation is $(x -3)^2 + (y - 2)^2 = 25.$
Presumably the values of $A, B, C$ were chosen to enable this simple solution.
I'll get me coat. :)
(Rephrasing Scott's answer.) Let the three points be $z_1$, $z_2$, $z_3$ as complex numbers. These lie on a circle $C_{old}$ given by $|z-c|^2=r^2$. We can map $C_{old}$ to another circle $C_{new}$ given by $|w-d|^2=s^2$ that is easier to work with using a linear transformation $f(z)=w=az+b$. In other words we'll have $f(C_{old})=C_{new}$ and moreover $f(c)=d$, the center of $C_{old}$ gets mapped by $f$ to the center of $C_{new}$. This should be intuitively clear since the linear transformation translates/scales/rotates $C_{old}$.
One choice of linear map is $f(z)=\frac{z-z_1}{z_2-z_1}$, with $f(z_1)=w_1=0$, $f(z_2)=w_2=1$. Evaluating the defining equation for $C_{new}$ at $w_1=0$, $w_2=1$, $w_3=f(z_3)$ gives \begin{align*} |d|^2&=s^2=d\bar{d},\\ |1-d|^2&=s^2=1-d-\bar{d}+|d|^2,\\ |w_3-d|^2&=s^2=|w_3|^2-w_3\bar{d}-\bar{w}_3d+|d|^2. \end{align*} Simplifying some gives the system of linear equations (in $d$, $\bar{d}$) \begin{align*} d+\bar{d}&=1\\ \bar{w}_3d+w_3\bar{d}&=|w_3|^2. \end{align*} Solving for $d$ gives the center of $C_{new}$ (assuming $w_3\neq\bar{w}_3$) $$ d=\frac{w_3-|w_3|^2}{w_3-\bar{w}_3}, $$ and applying $f^{-1}(w)=(z_2-z_1)w+z_1$ gives the center of $C_{old}$ $$ c=f^{-1}(d)=z_1+(z_2-z_1)d. $$
[Remark: If $w_3=\bar{w}_3$, then $z_3-z_1=w_3(z_2-z_1)$ for the real scalar $w_3$, so that $z_3$ is on the line through $z_1$ and $z_2$.]
Finally, the radius of $C_{old}$ can be found by evaluating $|z_i-c|^2=r^2$ at any of the initial three points.
Here is yet another implementation (Python):
def center_radius(z1, z2, z3):
"""
Return (center, radius), with complex center, of the circle through
complex z1, z2, z3. If z1, z2, z3 are colinear, print "colinear" and
return None.
"""
if z1 == z2:
print("colinear")
return None
def f(z):
return (z-z1)/(z2-z1)
def f_inv(w):
return z1+(z2-z1)*w
w3 = f(z3)
if w3.imag == 0: #may want abs(w3.imag) < eps
print("colinear")
return None
d = (w3-w3*w3.conjugate())/(w3-w3.conjugate())
c = f_inv(d)
r = abs(z1-c)
return (c, r)
For example,
>>> center_radius(0,1,complex(0,1))
((0.5+0.5j), 0.7071067811865476)
I know this question is outdated, but I like to show this easy solution.
First subtract the coordinates of a point from the others'. Now we are looking for the equation of a circle through the origin and two other points, which is of the form
$$(x-x_c)^2+(y-y_c)^2-r^2=x^2-2xx_c+y^2-2yy_c=0,$$ as there is no independent term.
This results in the $2\times2$ linear system of equations
(typo in the right hand side of the first equation of the 2x2 system changed wrong y2^2 to correct y1^2 ) $$\begin{cases}2x_cx_1+2y_cy_1=x_1^2+y_1^2,\\2x_cx_2+2y_cy_2=x_2^2+y_2^2.\end{cases}$$ giving the coordinates of the center and the radius $\sqrt{x_c^2+y_c^2}$. Don't forget to translate back by adding the first point.
Just for the sake of adding an answer:
Take any two points of the given three :$A(x_1,y_1)$ and $B(x_2,y_2)$. We will first determine the circle $\Omega$ having $AB$ as the diameter. This may seem arbitrary, but I assure you it is not.
Let $P(x,y)$ be any point on the circle. Then, $\angle APB$ is an angle in a semicircle and is hence a right angle. By writing $m_{AP}\cdot m_{BP}=-1$, we have $$\left(\frac{y-y_1}{x-x_1}\right)\left(\frac{y-y_2}{x-x_2}\right)=-1$$$$\implies (y-y_1)(y-y_2)+(x-x_1)(x-x_2)=0$$ is the equation of $\Omega$. Now, the straight line, $L$, through $A$ and $B$ is $$\frac{y-y_1}{x-x_1}= \frac{y_2-y_1}{x_2-x_1}$$$$\implies (y-y_1)(x_2-x_1)-(x-x_1)(y_2-y_1)=0$$
Now we will write the equation of a family of circles passing through $A$ and $B$. This can be written as $$\Omega +\lambda L=0$$ and this can be verified to being a circle for all $\lambda$. Now you can understand why we found the equation of $\Omega$: to find a similar (to the required curve) but special curve passing through the points of interest.
Now, put $(x_3,y_3)$ in the equation and solve for $\lambda$. This will give the required equation of the circle.
Equation of a circle with center $(x_0, y_0)$ and radius $r$ is
$$(x-x_0)^2+(y-y_0)^2=r^2$$
which can be expressed as
$$x^2+y^2+ax+by+c=0 \quad \quad \ldots (1)$$
where
$$\begin{array} xx_0 &=& -\frac{a}{2} \\ y_0 &=& -\frac{b}{2} \\ r &=& \sqrt{x_0^2 + y_0^2 - c} \quad \quad \ldots (2) \end{array}$$
Now, if the circle passes through points $(x_1,y_1)$, $(x_2,y_2)$ and $(x_3,y_3)$ (assuming that they are not collinear), then the equation $(1)$ is satisfied by the points, so that we have
$$\begin{array} aax_1 + by_1 + c &=& -x_1^2 -y_1^2 \\ ax_2 + by_2 + c &=& -x_2^2 -y_2^2 \\ ax_3 + by_3 + c &=& -x_3^2 -y_3^2 \end{array}$$
with solution
$$\begin{array} aa &=& \frac{det\left(\begin{bmatrix}-x_1^2-y_1^2 & y_1 & 1 \\ -x_2^2-y_2^2 & y_2 & 1 \\ -x_3^2-y_3^2 & y_3 & 1 \end{bmatrix}\right)}{det\left(\begin{bmatrix}x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{bmatrix}\right)} \\ b &=& \frac{det\left(\begin{bmatrix}x_1 & -x_1^2-y_1^2 & 1 \\ x_2 & -x_2^2-y_2^2 & 1 \\ x_3 & -x_3^2-y_3^2 & 1 \end{bmatrix}\right)}{det\left(\begin{bmatrix}x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{bmatrix}\right)} \\ c &=& \frac{det\left(\begin{bmatrix}x_1 & y_1 & -x_1^2-y_1^2 \\ x_2 & y_2 & -x_2^2-y_2^2 \\ x_3 & y_3 & -x_3^2-y_3^2 \end{bmatrix}\right)}{det\left(\begin{bmatrix}x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{bmatrix}\right)} \end{array}$$
Now, from $(2)$, we can compute the center $(x_0,y_0)$ and radius $r$ of the circle.
We can implement the above with the following r
function find.circle()
, which accepts the arguments $(x_1,y_1)$, $(x_2,y_2)$ and $(x_3,y_3)$ corresponding to the coordinates of the $3$ points the circle passes through, and returns the center $(x_0,y_0)$ and radius $r$ of the circle:
find.circle <- function(x1, y1, x2, y2, x3, y3) {
A <- matrix(c(x1, x2, x3, y1, y2, y3, 1, 1, 1), ncol=3)
b <- c(-x1^2-y1^2, -x2^2-y2^2, -x3^2-y3^2)
res <- solve(A, b)
a <- res[1]
b <- res[2]
c <- res[3]
print(c(a,b,c))
x0 <- -a/2
y0 <- -b/2
r <- sqrt(x0^2 + y0^2 - c)
return(c(x0,y0,r))
}
x1 <- 1
y1 <- 1
x2 <- 2
y2 <- 4
x3 <- 5
y3 <- 3
find.circle(x1, y1, x2, y2, x3, y3)
# [1] -6 -4 8 # (a, b, c)
# [1] 3.000000 2.000000 2.236068 # (x0, y0, r)
The equation of the circle is $x^2 + y^2 - 6x - 4y + 8 = 0$, i.e., $(x-3)^2+(y-2)^2=5$.
res <- find.circle(x1, y1, x2, y2, x3, y3)
x0 <- res[1]
y0 <- res[2]
r <- res[3]
plot(x0, y0, pch=19, xlim=c(x0-r,x0+r), ylim=c(y0-r,y0+r), col='brown')
for (theta in seq(0,2*pi,0.01)) {
points(x0 + r*cos(theta), y0 + r*sin(theta), pch=19, cex=0.3)
}
points(x1, y1, pch=19, xlim=c(x0-r,x0+r), ylim=c(y0-r,y0+r), col='red')
points(x2, y2, pch=19, xlim=c(x0-r,x0+r), ylim=c(y0-r,y0+r), col='green')
points(x3, y3, pch=19, xlim=c(x0-r,x0+r), ylim=c(y0-r,y0+r), col='blue')
grid()