<jupyter_start><jupyter_text>$\left(x-a\right)^2+\left(y-b\right)^2=r^2$ If $i$ and $j$ represent the horizontal and vertical distances from the starting point $P$ to the center of the circle, then $c_1 = p_1+i$ and $c_2=p_2+j$. That was convenenient. If, on the other hand, all we have is the radius of the circle... $$\left(p_1-a\right)^2+\left(p_2-b\right)^2=r^2$$ $$\left(q_1-a\right)^2+\left(q_2-b\right)^2=r^2$$ So, with numbers obtained from the `gcode` block in question capitalized, we need to solve the following system of equations for $a$ and $b$. $$\left(X_0-a\right)^2+\left(Y_0-b\right)^2=R^2$$ $$\left(X-a\right)^2+\left(Y-b\right)^2=R^2$$ $X_0$ and $Y_0$ represent the current tool position in the equation above, and are therefore known quantities at the time of computation.<jupyter_code>x, y, z, X, X0, Y, Y0, Z, Z0, R = symbols("x y z X X0 Y Y0 Z Z0 R")
e1 = (x-X0)**2 + (y-Y0)**2
e2 = (x-X)**2 + (y-Y)**2<jupyter_output><empty_output><jupyter_text>[Here](https://math.stackexchange.com/questions/1781438/finding-the-center-of-a-circle-given-two-points-and-a-radius-algebraically) is the answer from Stack Exchange. That question / answer also references [this](http://mathforum.org/library/drmath/view/54490.html) newsgroup message from who knows when?<jupyter_code>init_printing()
help(init_printing)
from enum import IntEnum

class CircularPlane(IntEnum):
    XY = 17
    ZX = 18
    YZ = 19
    

xa = (X - X0) / 2
ya = (Y - Y0) / 2
x0, y0 = X0 + xa, Y0 + ya
a = sqrt(xa**2 + ya**2)
b = sqrt(R**2 - a**2)
x3 = x0 + b*ya / a
x4 = x0 - b*ya / a
y3 = y0 - b*xa / a
y4 = y0 + b*xa / a
x3
y3
x3.simplify()
y3.simplify()
params = dict(zip([X0, Y0, X, Y, R], [1, 4, 5, 1, 13/2]))
x3 = x3.subs(params)
y3 = y3.subs(params)
x4 = x4.subs(params)
y4 = y4.subs(params)
T1, T2 = symbols("T1, T2")
C1 = Point3D(x3.subs(params), y3.subs(params))
x3, y3, x4, y4
from IPython.display import Image
Image("png/circles0000.png")<jupyter_output><empty_output><jupyter_text>$$AC_1=0$$<jupyter_code>args1 = [-x3, -y3]
args2 = [1/params[R]] * 3
Q0 = Point3D(params[X0], params[Y0], 0).translate(*args1).scale(*args2).evalf()
Q = Point3D(params[X], params[Y], 0).translate(*args1).scale(*args2).evalf()
print("1/R:", 1/params[R])
print("Q0:", Q0)
print("Q:", Q)<jupyter_output>1/R: 0.15384615384615385
Q0: Point3D(0.246153846153846, 0.969230769230769, 0)
Q: Point3D(0.861538461538462, 0.507692307692308, 0)<jupyter_text>Get the angles...<jupyter_code>theta1 = atan(Q0.y / Q0.x)
theta1
def rad2deg(t):
    return ((t % pi) * 180 / pi).evalf()

def deg2rad(t):
    return (pi * t / 180 % (pi / 2)).evalf()
from anglr import Angle
dir(Angle)
help(Angle.__init__)
P = Point3D()
Q
P.translate(*-Q)
rad2deg(theta1)
theta2 = atan(Q.y / Q.x)
rad2deg(theta2)
help(Angle.angle_between)
t1 = Angle(theta1)
t2 = Angle(theta2)
t3 = t1.angle_between_clockwise(t2)
from math import pi as PI
PI
t3.radians > PI
import matplotlib.path
matplotlib.path.Path.arc(t1.degrees, t2.degrees).vertices.transpose()[0]
help(map)
help(numpy.full)
import numpy
public(numpy)
public(numpy.ndarray)
help(numpy.ndarray.transpose)
public(matplotlib.path.Path)
from enum import Enum
class Quadrants(Enum):
    I = "I"
    II = "II"
    III = "III"
    IV = "IV"
    
class Directions(Enum):
    CW = 0
    CCW = 1

def quadrant(P):
    if P.x >= 0:
        if P.y > 0:
            return Quadrants.I
        else: return Quadrants.IV
    else:
        if P.y > 0:
            return Quadrants.II
        else:
            return Quadrants.III
        
def get_angle(P):
    return atan(P.y / P.x)

def rotate(*args, radians=(pi/2).evalf(), degrees=90):
    pass
        
def direction(P1, P2):
    """ Return the direction travelled to reach P2 from P1 via the shortest path.
    """
    t1 = get_angle(P1)
    t2 = get_angle(P2)
    if not quadrant(P1) == Quadrants.I:
        pass
    if quadrant(P1) == Quadrants.I:
        if quadrant(P2) == Quadrants.I:
            if t1 - t2 > 0:
                return Directions.CW
            else:
                return Directions.CCW
        if quadrant(P2) == Quadrants.II:
            return CCW
        if quadrant(P2) == Quadrants.III:
            if t2 < t1 + 180:
                return Directions.CCW
            else:
                return Directions.CCW
        if quadrant(P2) == Quadrants.IV:
            return Directions.CW
rotate
T1 = Matrix([[1/params[R], 0, -x3], [0, 1/params[R], -y3]])<jupyter_output><empty_output><jupyter_text>$\begin{matrix}  -1 & 3 \\  2 & -4 \end{matrix}$<jupyter_code>v3 = Matrix([x3, y3, 1])
v4 = Matrix([x4, y4, 1])
v3
T1 * v3
T1 * v4
import matplotlib.pyplot as plt
plt.plot([v3[0], v3[1]])
v3
pwd
from os import chdir
chdir('../../')
from py.startup import *
import mpl_toolkits.mplot3d
plt.ion()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')<jupyter_output><empty_output>