ready to publish code

A/4/testplot.py

@@ -1,5 +1,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
+import mpmath as mp
 
 kPitch = 1.7
 kAlpha = kPitch / (2 * np.pi)
@@ -8,7 +9,12 @@ kCriticalRadius = 4.5
 theta_max = (kCriticalRadius) / kAlpha + 2 * 2 * np.pi
 kPlotingRadius = theta_max * kAlpha
 
-kCriticalTheta = 2.86 / ((2 / 3) * kAlpha)
+def f(x):
+  r=(1/3)*kAlpha*mp.sqrt(1+x**2)
+  phi=mp.atan(x)
+  L=mp.mpf("2.86")
+  return (r+3*r*mp.cos(mp.pi-2*phi)-L)**2+(3*r*mp.sin(mp.pi-2*phi))**2-L**2
+kCriticalTheta = float(mp.findroot(f, 15, solver='secant'))
 # 生成角度数组
 theta = np.linspace(kCriticalTheta, theta_max, 1000)
 
@@ -80,7 +86,5 @@ theta_grid = np.arctan2(Y, X)
 valid_points = r_grid <= kPlotingRadius
 ax.scatter(theta_grid[valid_points], r_grid[valid_points], color='grey', s=10)  # 灰色小点
 
-plt.title("The Moving Path")
-
 # 显示图像
 plt.show()