Handling singularities in squircle parametric equations In order to plot a squircle using its parametric equations:

The equations can be found in page 9 of this
this article:

Fong, Chamberlain. &quot;Squircular Calculations.&quot; arXiv preprint arXiv:1604.02174 (2016).

Unfortunately the following singularities are not handled by those equations:

For t = 0, t = pi/2, t = pi, t = 3pi/2, and t=2*pi

Here is the code to plot these equations in Python
&quot;&quot;&quot;Plot squircle shape using its parametric equations&quot;&quot;&quot;
import numpy as np
from numpy import pi, sqrt, cos, sin, sign
import matplotlib.pyplot as plt

# FG stands for Fernandez-Guasti who introduced the algebric equations.
def FG_param(r, s, t):
    x = np.ones_like(t)
    y = np.ones_like(t)
    for i, _t in enumerate(t):
        if np.isclose(_t, 0.0):
            x[i] = r
            y[i] = 0.0
        elif np.isclose(_t, pi/2.0):
            x[i] = 0.0
            y[i] = r
        elif np.isclose(_t, pi):
            x[i] = -r
            y[i] = 0.0
        elif np.isclose(_t, 3*pi/2.0):
            x[i] = 0.0
            y[i] = -r
        elif np.isclose(_t, 2*pi):
            x[i] = r
            y[i] = 0.0
        else:
            x[i] = r*sign(cos(t[i]))/(s*sqrt(2)*np.abs(sin(t[i])))*sqrt(1 - sqrt(1-s**2*sin(2*t[i])**2))
            y[i] = r*sign(sin(t[i]))/(s*sqrt(2)*np.abs(cos(t[i])))*sqrt(1 - sqrt(1-s**2*sin(2*t[i])**2))
    return x, y

s = 0.7
r = 1.0
NTHETA = 90
t = np.linspace(0, 2*pi, NTHETA)

x, y = FG_param(r, s, t)

plt.gca().set_aspect('equal')
plt.scatter(x, y, s=5)
plt.show()

The function FG_param looks very clumsy. I appreciate any help to make it more compact and neat.
Thank you
Here's a different solution.  The singularities occur when sin(_t) or cos(_t) is zero.  So check for those conditions.

def FG_param(r, s, t):
    x = np.ones_like(t)
    y = np.ones_like(t)

    for i, _t in enumerate(t):
        sin_t = sin(_t)
        cos_t = cos(_t)

        if np.isclose(sin_t, 0.0):
            x[i] = -r if np.isclose(_t, pi) else r
            y[i] = 0.0

        elif np.isclose(cos_t, 0.0):
            x[i] = 0.0
            y[i] = r if _t &lt; pi else -r

        else:
            radical = sqrt(1 - sqrt(1 - s**2 * sin(2*_t)**2))

            x[i] = r*sign(cos_t) / (s*sqrt(2)*np.abs(sin_t)) * radical
            y[i] = r*sign(sin_t) / (s*sqrt(2)*np.abs(cos_t)) * radical

    return x, y