# 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. "Squircular Calculations." 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

"""Plot squircle shape using its parametric equations"""
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

• It is possible to simplify via knowing that boolean statements return 0 or 1 if they are coerced to; for example, np.isclose(_t, pi/2.0) returns 1 if true, so that can be inverted (display zero if matched) and then multiplied by the RHS to remove it; the equation can be restructured to be r + xsin(...), where x is the condition. This can be repeated for all remaining conditions and structured accordingly. I may provide a code sample – alexyorke Dec 6 '19 at 4:39
• Something sort of like this, but this is not 100% correct: x[i] = r * \ (np.isclose(_t, 0.0) or np.isclose(_t, 2*pi)) * \ (not np.isclose(_t, 0.0) or np.isclose(_t, 2*pi)) * \ (1 - (not np.isclose(_t, pi))) * \ (not np.isclose(_t, 3*pi/2.0)) – alexyorke Dec 6 '19 at 4:40

## 3 Answers

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 < 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

• The singularities occur also when s=0 – s.ouchene Dec 6 '19 at 12:13
• @Navaro: Exactly. In that case (i.e: s=0) the map will be a circle with radius r. – IamNotaMathematician Dec 7 '19 at 0:31
• @Navaro s==0 wasn't handled in the original code, so I didn't handle it either. It would be an easy case to handle. – RootTwo Dec 7 '19 at 1:08
• @RootTwo: Yes, you could just add something like this at the top of the function: if s == 0: return r*cos(t), r*sin(t) – s.ouchene Dec 7 '19 at 1:12

This solution is an order of magnitude faster than the other ones. It uses numpy functions, so there aren't any explicit loops. np.where chooses values that keep the equations from blowing up where the denominator would go to zero. Case for s == 0 is handled separately.

def FG_param1(r, s, t):

if s == 0:
return r*cos(t), r*sin(t)

sin_t = sin(t)
cos_t = cos(t)

radical = sqrt(1 - sqrt(1 - s**2 * sin(2*t)**2))

x_denom = np.where(np.isclose(sin_t, 0.0), 1.0, s*sqrt(2)*np.abs(sin_t))
x = r * sign(cos_t) * np.where(np.isclose(sin_t, 0.0, 1e-12), 1.0, radical / x_denom)

y_denom = np.where(np.isclose(cos_t, 0.0), 1.0, s*sqrt(2)*np.abs(cos_t))
y = r * sign(sin_t) * np.where(np.isclose(cos_t, 0.0, 1e-12), 1.0, radical / y_denom)

return x, y


Off the top of my head, and untested, how about:

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

for i, _t in enumerate(t):
radical = sqrt(1 - sqrt(1 - s**2 * sin(2*_t)**2))

try:
x[i] = r*sign(cos(_t)) / (s*sqrt(2)*np.abs(sin(_t))) * radical

except ZeroDivisionError:
x[i] = -r if np.isclose(_t, pi) else r

try:
y[i] = r*sign(sin(_t)) / (s*sqrt(2)*np.abs(cos(_t))) * radical

except ZeroDivisionError:
y[i] = r if _t < pi/2.0 else -r

return x, y