Recently I need to plot the

CodePudding user response：

From the link you provided, it looks like the range over which you are plotting the Cassini ovals change depending on how the ratio b/a compares to 1. For instance, when a<b, the range is whereas it is restricted to when a>=b. Similarly, when a>=b, the curve becomes two disjoint ovals while it is a single one when a<b. Consequently, in order to plot the disjoints ovals you need to plot all different solutions for r (changing signs etc.)

See code below where I plotted the Casini curves when a=2 and the b/a ratio goes from 0.1 to 1.5 (same ratio modulation than the figure in the link you provided)

import numpy as np
import matplotlib.pyplot as plt
from math import radians

a=2    
B= np.arange(0.2,3, 0.3)#b/a ration goes form 0.1 to 1.5
x=[]
y=[]
step=0.00001
lw=1.
colors = plt.cm.viridis(np.linspace(0,1,len(B)))[::-1]

#Function  to compute the first part of the Cassini ovals
def cas_process(rads):
  sin2=np.sin(2*rads)
  sin22=np.power(sin2,2)
  ba4=(b/a)**4
  diff=ba4-sin22
  sq=np.sqrt(diff)
  return sq

for b,i in zip(B,range(len(B))):    
  if a>=b:
    theta_0=0.5*np.arcsin((b/a)**2) 
    rads = np.arange(-theta_0, theta_0, step)[1:-1] #defining appropriate range for the b/a ratio
    tmp2=cas_process(rads)

    #Getting all solutions
    tmp3p=(np.cos(2*rads)) tmp2
    tmp3m=(np.cos(2*rads))-tmp2
    rpp=np.sqrt(tmp3p)*a
    rmm=np.sqrt(tmp3m)*(-a)
    rpm=np.sqrt(tmp3p)*(-a)
    rmp=np.sqrt(tmp3m)*(a)

    xpp=np.multiply(rpp, np.cos(rads))
    ypp= np.multiply(rpp, np.sin(rads))
    xpm=np.multiply(rpm, np.cos(rads))
    ypm= np.multiply(rpm, np.sin(rads))
    xmp=np.multiply(rmp, np.cos(rads))
    ymp= np.multiply(rmp, np.sin(rads))
    xmm=np.multiply(rmm, np.cos(rads))
    ymm= np.multiply(rmm, np.sin(rads))
    
    #Plotting disjoints ovals. Each of the curves below are necessary to get the complete picture
    plt.plot(xpp,ypp,color=colors[i],lw=lw,label='b/a=' str(np.round(b/a,2)))
    plt.plot(xmm,ymm,color=colors[i],lw=lw)
    plt.plot(xpm,ypm,color=colors[i],lw=lw)
    plt.plot(xmp,ymp,color=colors[i],lw=lw)

  elif a<b:
    rads = np.arange(0, (2 * np.pi), step) #defining appropriate range for the b/a ratio
    tmp2=cas_process(rads)
    tmp3=(np.cos(2*rads)) tmp2
    r=np.sqrt((tmp3))* (a)
    x=np.multiply(r, np.cos(rads))
    y= np.multiply(r, np.sin(rads))

    plt.plot(x,y,color=colors[i],lw=lw,label='b/a=' str(np.round(b/a,2)))
plt.legend(bbox_to_anchor=(1.1, 1.0))
plt.show()

and the output gives:

And below is the result you get for the range B= np.arange(1.9,2.1, 0.05) you specified in your question: