我无法在 python 中绘制卡西尼椭圆形

Question

最近需要绘制Cassini oval

但它只适用于 b/a>=1。这是我的简单代码。为了澄清，我使用了许多临时变量。

import numpy as np
import matplotlib.pyplot as plt
from math import cos, sin, sqrt, radians

# setting the values of a , b
a =2    
B= np.arange(1.9,2.1, 0.05)

# containing the radian values
rads = np.arange(0, (2 * np.pi), 0.03)

x=[]
y=[]
# plotting the ellipse
for b in B:    
    sin2=np.sin(2*rads)
    sin22=np.power(sin2,2)
    ba4=(b/a)**4
    tmp1=ba4-sin22
    tmp2=np.sqrt(tmp1)
    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,'g')
plt.show()

结果是

Answer 1

从您提供的 link 来看，您绘制卡西尼椭圆的范围似乎会根据比率 b/a 与 1 的比较而变化。例如，当 a<b时，范围为 whereas it is restricted to 时a>=b。同样，当a>=b时，曲线变成两个不相交的椭圆，而当a<b时，它是一个椭圆。因此，为了绘制不相交的椭圆，您需要绘制 r（改变符号等）

的所有不同解决方案

请参阅下面的代码，其中我在 a=2 和 b/a 比率从 0.1 变为 1.5 时绘制了卡西尼曲线（与您提供的 link 中的数字相同的比率调制）

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

并且输出给出：

下面是您在问题中指定的 B= np.arange(1.9,2.1, 0.05) 范围内得到的结果：

我无法在 python 中绘制卡西尼椭圆形

I can't plot Cassini oval in python

python

matplotlib