求助渐进多焦点镜片设计问题【附上全部Python代码】

OPO_salt_fish

问题：在定义子午线和轮廓线方程，并给出目标轮廓线对应的屈光度分布。但是通过sag方程（偏微分计算）得到的曲面屈光度分布很诡异，而且随着在点云上细分步数增加，求出的结果会更加奇怪，如下图：
         
从左到右，从上到下，分别为图(a), (b), (c)
其中，图a为步数=20的细分，图b为步数=30的细分，图c为步数=100的细分。

附上全部Python代码，有大佬，能帮我看看问题吗？小弟感激不尽。


import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from scipy.optimize import fsolve
from scipy.interpolate import bisplrep, bisplev
from scipy.interpolate import RectBivariateSpline

import json
import csv
import pandas as pd
import json

from scipy.integrate import quad
from scipy.signal import savgol_filter
from scipy import interpolate

from sympy import symbols, diff, integrate, Eq, solve, IndexedBase, Idx, linsolve, Rational, Integral
from sympy import Sum
from sympy import oo, integrate
from sympy import sin, csc, cos, sec, tan, cot, acsc, acos, asec, atan, acot, atan2, asin
from sympy import lambdify
from sympy import Abs, re

from sympy import init_printing
init_printing()

# ---------------------------------------------------------------------------------------
# 重写坐标系 并调整单位到mm， 20231210
# -----------------

# m, l, h, L = 3, 3, 10, 10
#定义求解需要的方程
#定义由边界条件得出的c_n， h, n， m, l的方程关系
# 视远点镜片度数100度‚视远与视近点之间加光度数为150度；
# 镜片通道长度 h＝21mm；
# 视远点 PA # 与镜片中心 O点的间距 L＝7mm。
# 镜片直径60mm，折射率 n＝1．499。

###########################################
h_fix = 14
L_fix = 7
# u_fix = 0.014 # 应该是约定的h-L
##########################################

solutions = []

for m in range(1,40):
    for l in range(1,40):
        # #第一行
        first_row = [ h_fix**n for n in range(m, m+l)]

        # 将第一行转换为 NumPy 数组（也就是矩阵的一部分）
        first_row_matrix = np.array(first_row).reshape(1,-1)

        # first_row_matrix = np.ones((1, l))

        #下面的l-1 行
        matrix = np.zeros((l-1,l))
        for i in range(1, l):
            for j in range(1, l+1):
                product = (h_fix)**(m+j-1) #和行无关，只和列有关
                for k in range(m+j-i, m+j-1+1):
                    product *= k  
                matrix[i-1][j-1] = product

        # 补充一列 0
        add_col = np.zeros((l, 1))
        add_col[0][0] = 1

        # 构建增广矩阵
        first_row_matrix = first_row_matrix.astype(float)
        matrix = matrix.astype(float)
        add_col = add_col.astype(float)

        Augmented_matrix = np.vstack((first_row_matrix, matrix))
        Augmented_matrix = np.hstack((Augmented_matrix, add_col))


        # 求解
        solution = np.linalg.solve(Augmented_matrix[:,0:-1], Augmented_matrix[:,-1])
        solution = solution.tolist()
        solutions.append(solution)


#定义方程
#定义子午线上关键点的x轴坐标的位置u
u = symbols('u')

#定义视远点PA的曲率半径r_d
r_D = symbols('r_d')

#定义视近点PB的曲率半径r_r
r_R = symbols('r_r')

# 定义需要求解的变量c_n，和对应的连加变量n
c = IndexedBase('c')
n = Idx('n')

#设定corridor高度h和视远点长度L
h, L = symbols('h L')

#定义后续用到的高阶微分的第一个非零的阶数
m, l = symbols('m l')

#设定微分次数d_n
d_n = symbols('d_n')

#输出子午线曲率设计函数，注意不是曲率半径
fc_u = 1/r_D + (1/r_R - 1/r_D) * Sum(c[n] * (u+L)**n, (n, m, m+l-1) )

#输出f_r_u的高阶微分平滑性评价函数
dff_fc_u = diff(fc_u, u, d_n)

#控制输入单位是mm
#定义材料的折射率是1.53
ref_n = symbols('ref_n')

#定义求解需要的方程
#定义由边界条件得出的c_n， h, n， m, l的方程关系
# 视远点镜片度数100度‚视远与视近点之间加光度数为150度；
# 镜片通道长度 h＝21mm；
# 视远点 PA # 与镜片中心 O点的间距 L＝7mm。
# 镜片直径60mm，折射率 n＝1．499。

# 1/f = (n-1)/R
#镜片参数
D_lens = 60
# n_ref = 1.499
n_ref = 1.499
D_D = 1
D_R = 2.5


r_D_0 = (n_ref - 1) *1000 / D_D

r_R_0 = (n_ref - 1) *1000 / D_R

print(f"r_D_0 = {r_D_0:.3f}")
print(f"r_R_0 = {r_R_0:.3f}")

fc_u_practical = fc_u.subs({r_D:r_D_0, r_R:r_R_0, L: L_fix})

fr_u = 1/fc_u_practical

# 带入实际案例测试
########################################
m_0 = 3
l_0 = 3
#######################################
fr_u_0 = fr_u.subs({m:m_0,l:l_0}).doit()
fc_u_0 = fc_u_practical.subs({m:m_0,l:l_0}).doit()


c_n = []
col_length = len(solutions[-1])
current_id = (m_0-1)*(col_length)+l_0 -1

arr=solutions[current_id]

for id, item in enumerate(arr):
    c_n.append(item)

for id, item in enumerate(c_n):
    fr_u_0 = fr_u_0.subs(c[id+m_0],item)
    fc_u_0 = fc_u_0.subs(c[id+m_0],item)


c1 = fc_u_0.subs({u:0})
d2 = fc_u_0.subs({u:L_fix}) * (n_ref - 1) * 1000
d3 = fc_u_0.subs({u:- h_fix + L_fix}) * (n_ref - 1) * 1000
print(f"test, 0点曲率: {c1:.3f}, 视近点D={d2:.3f}, 视远点D={d3:.3f} ")
print(f"当前曲率和u的函数》 {fc_u_0}")
print('\n')
dff_1 = diff(fc_u_0, u, 1).subs(u, h_fix - L_fix)
print(f"在h-L处的1阶导数 = {dff_1.round(5)}")
# ///////////////////////////////////////
# 按照微分几何定义和子午线曲率方程fr_u 求出镜片上个点的矢高，sag equation

#plot 子午线曲率图，Y轴为P，diopter， X轴为u -----------------------------------

# # 创建一个新的图形
# fig, ax = plt.subplots()
fig_1 = plt.figure()
ax = fig_1.add_subplot(111)

ax.set_title(f'm={m_0}, l={l_0}, meridian perf. D, u=x')
ax.set_xlabel("u, or x")
ax.set_ylabel("Diopter")

x_plot = np.linspace( -(L_fix + 1) , h_fix - L_fix + 1, 100)
## 这里的x_plot后续被覆盖


fd_u = (n_ref-1) * 1000 / fr_u_0

fd_u_plot =  lambdify(u, fd_u, "numpy")
y_plot = fd_u_plot(x_plot)
ax.plot(x_plot, y_plot, label='Diapter = (n_ref-1)/ fr_u_0')



point_D = fd_u.subs(u, -L_fix)
ax.plot(-L_fix, point_D, 'o')
# 画出十字线
ax.axhline(y=point_D, color='r', linestyle='--', alpha = 0.2)
ax.axvline(x=-L_fix, color='r', linestyle='--', alpha = 0.2)

# 标注点的坐标
ax.annotate(f'D{-L_fix}, {"%.3f" % point_D})', (-L_fix, point_D), textcoords="offset points", xytext=(0,10), ha='center')

point_D = fd_u.subs(u, h_fix - L_fix)
ax.plot(h_fix - L_fix, point_D, 'o')
# 画出十字线
ax.axhline(y=point_D, color='r', linestyle='--', alpha = 0.2)
ax.axvline(x=h_fix - L_fix, color='r', linestyle='--', alpha = 0.2)

# 标注点的坐标
ax.annotate(f'R{round((h_fix - L_fix), 3)}, {point_D.round(3)})', (h_fix - L_fix, point_D), textcoords="offset points", xytext=(0,10), ha='center')

ax.grid(True)
# 添加图例
ax.legend()
plt.show()

# ------------------------------------------------------------------
# 定义轮廓线方程
# 定义sag，z, z(x, y)
x, y, b = symbols('x y b')
fu_xy = (L**2 + b**2) * (-y) / ( (Abs(x)+b)**2 + y**2 )

########################################################################
b_val = 15
# L_val = 0.007
fu_xy_0 = fu_xy.subs({L: L_fix, b: b_val })
#######################################################################

# -------------------------------------------------------------------------------------
# 单独设定一下3D画图的精细度，步数
steps = 100
#######################################################################

# 创建一个新的图形子图143 u(x,y) performance

fig_2 = plt.figure(figsize=(24,9))
ax1_2 = fig_2.add_subplot(131, projection='3d')

ax1_2.set_title(f'm={m_0}, l={l_0}, u(x,y) performance')
ax1_2.set_xlabel("x")
ax1_2.set_ylabel("y")
ax1_2.set_zlabel("u(x,y)")

fu_xy_0_plot =  lambdify([x, y], fu_xy_0, "numpy")

########################################################
x_plot = np.linspace(-D_lens/2, D_lens/2, steps)
y_plot = np.linspace(-D_lens/2, D_lens/2, steps)
########################################################

x_plot, y_plot = np.meshgrid(x_plot, y_plot)
u_plot = fu_xy_0_plot(x_plot, y_plot)

# 使用 clip 函数将数组中的元素限制在 -0.007 到 0.014 的范围内
u_plot_clip = np.clip(u_plot, -L_fix, h_fix-L_fix)

# 绘制一个三维曲面
surf = ax1_2.plot_surface(x_plot, y_plot, u_plot_clip, cmap='viridis', label = 'u(x,y)', alpha = 0.5)

# 绘制等高线图
conto = ax1_2.contour(x_plot, y_plot, u_plot_clip, 20, zdir='z',  cmap='viridis')

# # 添加颜色条
plt.colorbar(surf, shrink = 0.3, aspect= 10 )

# -------------------------------------------------------
# 尝试用python画一下镜片焦度等值线

D_actual =  fd_u_plot(u_plot_clip)
ax1_3 = fig_2.add_subplot(132, projection='3d')

ax1_3.set_title(f'm={m_0}, l={l_0}, b={b_val}, D(Sphere) perf. of desired lens')
ax1_3.set_xlabel("x")
ax1_3.set_ylabel("y")
ax1_3.set_zlabel("D")

# 绘制一个三维曲面
ax1_3.plot_surface(x_plot, y_plot, D_actual, \
                            cmap='viridis', label = 'sag performance', alpha = 0.5)

# 绘制等高线图
conto = ax1_3.contour(x_plot, y_plot, D_actual, 20, \
                    zdir='z',  cmap='viridis')

# # 添加颜色条
plt.colorbar(conto, shrink = 0.3, aspect= 10 )

# 2D-------------------------------------------------------

ax1_4 = fig_2.add_subplot(133)
# 绘制等高线图
conto = ax1_4.contour(x_plot, y_plot, D_actual, \
                      10, zdir='z',  cmap='viridis')
# # 添加颜色条
plt.colorbar(conto, shrink = 0.3, aspect= 10 )

ax1_4.set_title(f'm={m_0}, l={l_0}, b={b_val}, D(Sphere) perf. actual lens')
ax1_4.set_xlabel("x")
ax1_4.set_ylabel("y")
ax1_4.set_aspect('equal')
ax1_4.grid()



plt.show()



# 定义好子午线，定义好轮廓线之后，再处理sag ——————————————————————————————————————————

def func_fz_sag(u_for_sag, x, y):
    # u_for_sag = 0.005
    # 定义子午线上曲率中心坐标， xi, eta, zeta
    theta = asin( integrate(1/fr_u_0, (u, 0, u )) ).subs(u, u_for_sag)
    ru_x = 0
    ru_y = - (u - fr_u_0 * sin(theta)).subs(u, u_for_sag)
    ru_z = (fr_u_0*cos(theta) + integrate(tan(theta), (u, 0, u))).subs(u, u_for_sag)

    arr = [ru_x, ru_y, ru_z]
    fz = ru_z - (fr_u_0**2 - (x - ru_x)**2 - (y - ru_y)**2 ) ** Rational(1/2)
    z_0 = fz.subs({u:u_for_sag, x:x, y:y})
    return z_0, arr

# ----------------------
# 画出所有的sag, 画出面型

sag_0 = []
# center_0 = []
for i in range(u_plot_clip.shape[0]):
    for j in range(u_plot_clip.shape[1]):
        print(i,j)
        row = []
        z_0, r_center_0 = func_fz_sag(u_plot_clip[i,j], x_plot[i,j], y_plot[i,j]) # 这步计算很慢
        sag_0.append(z_0)

sag_0 = np.array(sag_0).reshape(steps, steps)
np.savetxt("sag_0.csv", sag_0, delimiter=',')
# center_0 = np.array(center_0).reshape(-1,3)

# 创建一个新的图形 sag performance

fig_3 = plt.figure(figsize=(8,8))
ax2_1 = fig_3.add_subplot(111, projection='3d')

ax2_1.set_title(f'm={m_0}, l={l_0}, b={b_val}, sag performance')
ax2_1.set_xlabel("x")
ax2_1.set_ylabel("y")
ax2_1.set_zlabel("sag")

# 绘制一个三维曲面
################ 这里是核心数据 x_plot, y_plot, sag_0 #######################
ax2_1.plot_surface(x_plot, y_plot, sag_0, \
                   cmap='viridis', label = 'sag performance', alpha = 0.8)
############################################################################

# 在画一个保准的边界u对应的度数的镜片 -----------------------
# 取 u = -0.007 即-L_fix ， 即远视力点，画一个基准曲面

x_base = -L_fix
y_base = 0
u_base = fu_xy_0.subs({x:x_base, y:y_base})
r_base = fr_u_0.subs({u:u_base})

z_base_0, r_center_0 = func_fz_sag(u_base, x_base, y_base)

def func_z_base(x,y):
    z = r_center_0[2]-(r_base**2 - (x-r_center_0[0])**2 - (y - r_center_0[1])**2)**0.5
    return z

sag_base=[]
for i in range(u_plot_clip.shape[0]):
    for j in range(u_plot_clip.shape[1]):
        print(i,j)
        row = []
        z_0 = func_z_base(x_plot[i,j], y_plot[i,j])
        sag_base.append(z_0)

sag_base = np.array(sag_base).reshape(steps, steps)
# 绘制一个三维曲面base sag
ax2_1.plot_surface(x_plot, y_plot, sag_base, cmap='viridis', label = 'base sag', alpha = 0.5)

# plt.legend()
plt.show()


# -------------------------------------------------------------------------------
# 从sag推算的光焦度分布-拟合前的图

# 计算 x_plot 和 y_plot 的间距
dx = x_plot[0][1] - x_plot[0][0]
dy = y_plot[1][0] - y_plot[0][0]

# 计算 sag_0 的梯度

# ##拟合前
Zy, Zx = np.gradient(sag_0, dy, dx)
fit_mark = "Before Fit"

# # 拟合后
# Zy, Zx = np.gradient(sag_new, dy, dx)
# fit_mark = "After Fit"

Zxy, Zxx = np.gradient(Zx, dy, dx)
Zyy, _ = np.gradient(Zy, dy, dx)

#换一个求法
# 直接求H，K

E = 1 + Zx**2
F = Zx * Zy
G = 1 + Zy**2
L = Zxx / np.sqrt(1 + Zx**2 + Zy**2)
M = Zxy / np.sqrt(1 + Zx**2 + Zy**2)
N = Zyy / np.sqrt(1 + Zx**2 + Zy**2)

H = (E*N + G*L - 2*F*M) / (2*(E*G - F**2))
K = (L*N - M**2) / (E*G - F**2)
kappa_1 = H - np.sqrt(H**2 - K)
kappa_2 = H + np.sqrt(H**2 - K)

alpha_1 = np.arctan(-(M - kappa_1*F) / (N - kappa_1*G))
alpha_2 = np.arctan(-(M - kappa_2*F) / (N - kappa_2*G))

AOSA = alpha_2

P_1 = (n_ref - 1)*1000 * kappa_1
P_2 = (n_ref - 1)*1000 * kappa_2

SA = np.abs(P_1 - P_2)
SP = (P_1 + P_2)/2
# SP = P_2
print(f"SP.max():{SP.max()}, SP.min():{SP.min()}")

# 出一下图
# 创建一个新的图形窗口

fig_4 = plt.figure(figsize=(16,16))
ax_4_1 = fig_4.add_subplot(221)
# 第一个图，绘制颜色带代表的 SP，再画出轮廓线
cs1 = ax_4_1.contourf(x_plot, y_plot, SP, cmap='jet')
contours_1 = ax_4_1.contour(x_plot, y_plot, SP, colors='black', linewidths=0.5)
# 添加轮廓线标签
plt.clabel(contours_1, inline=True, fontsize=8)

ax_4_1.set_title(f'[{fit_mark}], m={m_0}, l={l_0}, b={b_val}, steps={steps}, SP')
ax_4_1.set_xlabel("x")
ax_4_1.set_ylabel("y")
ax_4_1.set_aspect('equal')

# -------------------------------------------
# 第二个图，绘制颜色带代表的 SA，再画出轮廓线
ax_4_2 = fig_4.add_subplot(222)
cs2 = ax_4_2.contourf(x_plot, y_plot, SA, cmap='jet')
contours_2 = ax_4_2.contour(x_plot, y_plot, SA, colors='black', linewidths=0.5)
# 添加轮廓线标签
plt.clabel(contours_2, inline=True, fontsize=8)


# 为了画出AOSA方向，我们可以创建一个向量场。为了简化图形，我们只在每个5个点上画一个箭头
aosa_u = np.cos(AOSA)[::2, ::2]  # 长度和方向的cos值
aosa_v = np.sin(AOSA)[::2, ::2]  # 长度和方向的sin值
ax_4_2.quiver(x_plot[::2, ::2], y_plot[::2, ::2], aosa_v, aosa_u, color='white')

ax_4_2.set_title(f'[{fit_mark}], m={m_0}, l={l_0}, b={b_val}, steps={steps}, SA')
ax_4_2.set_xlabel("x")
ax_4_2.set_ylabel("y")


ax_4_2.set_aspect('equal')






#-----------------------------------------------------------------------------------
# 拟合处理一下x_plot，y_plot, sag_0

x_plot = x_plot.astype('float64')
y_plot = y_plot.astype('float64')
sag_0 = sag_0.astype('float64')

# 创建一个5阶的RectBivariateSpline对象
degree = 5
spl = RectBivariateSpline(x_plot[0, :], y_plot[:, 0], sag_0, kx=degree, ky=degree)


x_fine = np.linspace(-D_lens/2, D_lens/2, steps*1)
y_fine = np.linspace(-D_lens/2, D_lens/2, steps*1)

x_fine, y_fine = np.meshgrid(x_fine, y_fine)

# 计算新网格上的高度
sag_new = spl(x_fine[0, :], y_fine[:, 0])

# 拟合后的图----------------------------------------------

# 计算 x_fine 和 y_fine 的间距
dx = x_fine[0][1] - x_fine[0][0]
dy = y_fine[1][0] - y_fine[0][0]

# 计算 sag_0 的梯度

# # ##拟合前
# Zy, Zx = np.gradient(sag_0, dy, dx)
# fit_mark = "Before Fit"

# 拟合后
Zy, Zx = np.gradient(sag_new, dy, dx)
fit_mark = "After Fit"

Zxy, Zxx = np.gradient(Zx, dy, dx)
Zyy, _ = np.gradient(Zy, dy, dx)

#换一个求法
# 直接求H，K

E = 1 + Zx**2
F = Zx * Zy
G = 1 + Zy**2
L = Zxx / np.sqrt(1 + Zx**2 + Zy**2)
M = Zxy / np.sqrt(1 + Zx**2 + Zy**2)
N = Zyy / np.sqrt(1 + Zx**2 + Zy**2)

H = (E*N + G*L - 2*F*M) / (2*(E*G - F**2))
K = (L*N - M**2) / (E*G - F**2)
kappa_1 = H - np.sqrt(H**2 - K)
kappa_2 = H + np.sqrt(H**2 - K)

alpha_1 = np.arctan(-(M - kappa_1*F) / (N - kappa_1*G))
alpha_2 = np.arctan(-(M - kappa_2*F) / (N - kappa_2*G))

AOSA = alpha_2

P_1 = (n_ref - 1)*1000 * kappa_1
P_2 = (n_ref - 1)*1000 * kappa_2

SA = np.abs(P_1 - P_2)
SP = (P_1 + P_2)/2
# SP = P_2
print(f"SP.max():{SP.max()}, SP.min():{SP.min()}")

# 出一下图
# 创建一个新的图形窗口


# 第一个图，绘制颜色带代表的 SP，再画出轮廓线
ax_4_3 = fig_4.add_subplot(223)
cs3 = ax_4_3.contourf(x_fine, y_fine, SP, cmap='jet')
contours_3 = ax_4_3.contour(x_fine, y_fine, SP, colors='black', linewidths=0.5)
# 添加轮廓线标签
plt.clabel(contours_3, inline=True, fontsize=8)

ax_4_3.set_title(f'[{fit_mark}], m={m_0}, l={l_0}, b={b_val}, steps={steps}, SP')
ax_4_3.set_xlabel("x")
ax_4_3.set_ylabel("y")

ax_4_3.set_aspect('equal')




# 第二个图，绘制颜色带代表的 SA，再画出轮廓线
ax_4_4 = fig_4.add_subplot(224)
cs4 = ax_4_4.contourf(x_fine, y_fine, SA, cmap='jet')
contours_4 = ax_4_4.contour(x_fine, y_fine, SA, colors='black', linewidths=0.5)
# 添加轮廓线标签
plt.clabel(contours_4, inline=True, fontsize=8)

# 为了画出AOSA方向，我们可以创建一个向量场。为了简化图形，我们只在每个5个点上画一个箭头
aosa_u = np.cos(AOSA)[::2, ::2]  # 长度和方向的cos值
aosa_v = np.sin(AOSA)[::2, ::2]  # 长度和方向的sin值
ax_4_4.quiver(x_fine[::2, ::2], y_fine[::2, ::2], aosa_v, aosa_u, color='white')

ax_4_4.set_title(f'[{fit_mark}], m={m_0}, l={l_0}, b={b_val}, steps={steps}, SA')
ax_4_4.set_xlabel("x")
ax_4_4.set_ylabel("y")


ax_4_4.set_aspect('equal')



#---------------------------------------------------------------

OPO_salt_fish

补充一下子午线和轮廓线的设计逻辑
子午线设计采用的核心方程是： fc_u = 1/r_D + (1/r_R - 1/r_D) * Sum(c[n] * (u+L)**n, (n, m, m+l-1) )
轮廓线的核心方程是： fu_xy = (L**2 + b**2) * (-y) / ( (Abs(x)+b)**2 + y**2 )

附上子午线设计对应的论文解释，请见图片。