I am trying to fit a set of data to an off-center-circle using python. However, the displacement of the center of the circle (to the center of the axis) is evaluated negative, when it should be positive.

Consider a circle of radius $a$ with center at $(theta_0, r_0)$. These are the equations I developed (I suppose everything is right):

Point $(theta, r) on the circle:

Partial derivative of $r$ with respect to $a$:

Partial derivative of $r$ with respect to $r_0$:

Partial derivative of $r$ with respect to $theta_0$:

This is a MWE:

import numpy as np               #... for mathematical operations and constants
from scipy import optimize       #... for data fitting
import matplotlib.pyplot as plt  #... for plots


dsts = [1.0238822745042868, 1.0179614005977726, 1.0288142236514832, 1.014240255062302, 1.0297548330833965]
angs = [1.2118767129907884, 0.46576418683129606, -0.2945355605442553, 0.7614519513889376, -0.7442328029712463]


#--> DEFINE FUNCTIONS TO FIT THE DATA
#... see: https://scipython.com/book/chapter-8-scipy/examples/non-linear-fitting-to-an-ellipse/

def circle(theta, p): #... defines the function for circle
    a, r0, theta_0 = p
    rc = r0*np.cos(theta - theta_0)
    return ( rc + (rc**2 + a**2 - r0**2)**0.5 )

def residuals(p, r, theta):
    """ Return the observed - calculated residuals using f(theta, p). """
    return r - circle(theta, p)

def jac(p, r, theta):
    """ Calculate and return the Jacobian of residuals. """
    tmp = circle(theta, p)
    a, r0, theta_0 = p
    da  = a/( tmp - r0*np.cos(theta - theta_0) )
    dr  = ( r*np.cos(theta - theta_0) - r0 )/( r - r0*np.cos(theta - theta_0) )
    dt  = tmp*r0*np.sin(theta - theta_0)/( tmp - r0*np.cos(theta - theta_0) )
    
    return -da,  -dr, -dt
    return np.array((-da, -dr, -dt)).T


#--> FIT THE DATA
p0 = (1, 0.05, 0) #... initial guess
plsq = optimize.leastsq(residuals, p0, Dfun=jac, args=(dsts, angs), col_deriv=True)


#--> CREATE FITTING CURVE ARRAYS
t_fit = np.linspace(0, 2*np.pi, 100)
r_fit = circle(t_fit, plsq[0])


#--> PLOT
cm = 1/2.54  #... centimeters-to-inches conversor
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}, figsize=(8.5*cm, 8.5*cm), dpi=300)

#--> Write the title
fit_a = str(round(plsq[0][0], 4))
fit_r = str(round(plsq[0][1], 4))
fit_t = str(round(plsq[0][2]*180/np.pi, 1))
pltTitle = "$a =$" + fit_a + "; $r_0 =$ " + fit_r + "; $\theta_0 =$ " + fit_t + "$^circ$"
ax.set_title(pltTitle, va='bottom', fontsize=8)

#--> Plot the measured data
ax.scatter(angs, dsts)

#--> Plot the fitted orbit
ax.plot(t_fit, r_fit, linestyle='dotted')

#--> Configure the axis
ax.set_rmax(1.2)
ax.set_rticks(np.linspace(0.2, 1, 5))
ax.set_xticklabels([])      #... remove angular tick labels
ax.set_yticklabels([])      #... remove radial tick labels
ax.grid(True)

#--> Display the plot
plt.show()

Any suggestion to get positive $r_0$?