Plotting polar function using matplotlib

I'm trying to plot this function using matplotlib.

As you can see in the Desmos app, the equation correctly plot the function as circle, but when I try to port it to Python, I got this instead:

import numpy as np
import matplotlib.pyplot as plt

def fungsi_r4(theta, theta0, r0, a):
  return r0 * np.cos(theta - theta0) + np.sqrt((a ** 2) - (r0 ** 2) * (np.sin(theta - theta0) ** 2))

theta = np.linspace(0, 2 * np.pi, 100)
r = fungsi_r4(theta, 2.4, 5.1, 2.6)

ax = plt.subplot(projection='polar')
ax.plot(theta, r)

My feeling tells me it has something to do with the negative values retuned from the function but I don't know what to do with it.

The difference is in how the two programs handle negative radii: Desmos flips them back through the origin, while matplotlib extends the radial scale into the negatives.

Here's some code to modify the data points like Desmos does:

def flip_negative_radii(theta, r):
    flip_mask = r < 0
    r[flip_mask] *= -1
    theta[flip_mask] = (theta[flip_mask] - np.pi) % (2*np.pi)
    return theta, r

Example usage:

import numpy as np
import matplotlib.pyplot as plt

def fungsi_r4(theta, theta0, r0, a):
    return r0 * np.cos(theta - theta0) + np.sqrt((a ** 2) - (r0 ** 2) * (np.sin(theta - theta0) ** 2))

theta = np.linspace(0, 2 * np.pi, 100)
r = fungsi_r4(theta, 2.4, 5.1, 2.6)

ax = plt.subplot(projection='polar')
ax.plot(*flip_negative_radii(theta, r))
plt.show()

The gaps are due to the term in the square root going negative and producing NaNs. If possible, I would try to come up with a parametric representation of the function, which would avoid both of these problems.