import numpy as np
'''
---L、U计算模块---
包含:
Calc(...), 给定孔径系数、视场系数、物面、镜面、入瞳信息,光线类型,利用迭代递推模块进行中间计算,返回各面L、U、I、PA信息。
'''
def Calc(K_eta, K_W, object, lens, pupil, image, raytype):
"""
作用:计算某一孔径、视场光线在逐个面的(L,U)坐标数据
输入:K_eta为孔径取点系数,K_W为视场取点系数
"""
#初始化(L,U)坐标列表
L0 = []
L1 = []
U = []
I0 = []
I1 = []
PA = []
A = pupil.semi_dia #入瞳半径
eta = K_eta * A #像高
#物在无穷远且光线与光轴平行时的情形
if(object.thickness == None and K_W == 0):
h0 = K_eta * pupil.semi_dia
L0, L1, U, I0, I1, PA = lens_recursion(0, 0, h0, lens, raytype, 1)
#其余情形 -> 转化为轴上点在有限远处
else:
if((object.thickness == None and K_W != 0)):
U = -K_W * object.semi_dia * np.pi/180
L = pupil.thickness + eta/np.tan(U)
elif(object.semi_dia == 0):
L = -object.thickness
U = np.arcsin(K_eta * A/np.sqrt(A**2 + Lp**2))
elif(object.semi_dia != 0):
Lp = pupil.thickness #第一面到入瞳的距离
y = K_W * object.semi_dia #物方线视场
U = np.arctan((y-eta)/(Lp+object.thickness))
L = Lp + eta/((y-eta)/(Lp+object.thickness))
L0, L1, U, I0, I1, PA = lens_recursion(L, U, 0, lens, raytype, 0)
return L0, L1, U, I0, I1, PA
'''
---迭代递推模块---
包含:
lens_recursion(...), 给定入射光线L、U,对各镜面进行光线追迹并返回各面L、U数据
astigmatism_recursion(...), 给定各面L、U、I、PA参数,求子午场曲,弧矢场曲,像散
'''
EPS = 1e-5 #PA校验精度
def lens_recursion(L, U, H, lens, raytype, isParallel):
"""
函数作用: 给定轴上光线入射初值, 对每个镜面进行递推
输入: 初始物方截距(L)、物方倾斜角(U)、镜面数据(lens)、光线类型(raytype = 'd', 'F', 'C')
是否平行于主光轴(isParallel, 1平行, 0倾斜)
输出: 光线经过每面镜面后的L, L', U', I, I', PA数组
"""
snumber = len(surface)
recursion = 0
sinI = 0; sinII = 0
I = 0; II = 0
UU = 0
LL = 0
L_data = [L]
LL_data = []
UU_data = [U]
I_data = []
II_data = []
PA_data = []
if(raytype == 'd'): #d光,折射率取nd
for recursion in range(lensnumber):
if(recursion == 0 and isParallel):
sinI = H / lens[recursion].radius
else:
sinI = (L - lens[recursion].radius) / lens[recursion].radius * np.sin(U)
sinII = lens[recursion].nd0 / lens[recursion].nd1 * sinI
if(np.abs(sinI) < 1 and np.abs(sinII) < 1):
I = np.arcsin(sinI); II = np.arcsin(sinII)
I_data.append(I); II_data.append(II)
UU = U + I - II
LL = lens[recursion].radius + lens[recursion].radius * sinII / np.sin(UU)
LL_data.append(LL); UU_data.append(UU)
#PA校对法
PA1 = L * np.sin(U) / np.cos((I-U)/2)
PA2 = LL * np.sin(UU) / np.cos((II-UU)/2)
if(np.abs(PA1-PA2) > EPS and isParallel == 0):
print("Recursion not right.")
print(f'PA1 = {PA1}, PA2 = {PA2}')
PA_data.append(PA2)
#过渡公式
U = UU
L = LL - lens[recursion].thickness
#舍去最后一个过渡的物方截距
if(recursion != lensnumber - 1):
L_data.append(L)
if(raytype == 'F'): #F光,折射率取nF
for recursion in range(lensnumber):
if(recursion == 0 and isParallel):
sinI = H / lens[recursion].radius
else:
sinI = (L - lens[recursion].radius) / lens[recursion].radius * np.sin(U)
sinII = lens[recursion].nF0 / lens[recursion].nF1 * sinI
if(np.abs(sinI) < 1 and np.abs(sinII) < 1):
I = np.arcsin(sinI); II = np.arcsin(sinII)
I_data.append(I); II_data.append(II)
UU = U + I - II
LL = lens[recursion].radius + lens[recursion].radius * sinII / np.sin(UU)
LL_data.append(LL); UU_data.append(UU)
#PA校对法
PA1 = L * np.sin(U) / np.cos((I-U)/2)
PA2 = LL * np.sin(UU) / np.cos((II-UU)/2)
if(np.abs(PA1-PA2) > EPS and isParallel == 0):
print("Recursion not right.")
print(f'PA1 = {PA1}, PA2 = {PA2}')
PA_data.append(PA2)
#过渡公式
U = UU
L = LL - lens[recursion].thickness
#舍去最后一个过渡的物方截距
if(recursion != lensnumber - 1):
L_data.append(L)
if(raytype == 'C'): #C光,折射率取nC
for recursion in range(lensnumber):
if(recursion == 0 and isParallel):
sinI = H / lens[recursion].radius
else:
sinI = (L - lens[recursion].radius) / lens[recursion].radius * np.sin(U)
sinII = lens[recursion].nC0 / lens[recursion].nC1 * sinI
if(np.abs(sinI) < 1 and np.abs(sinII) < 1):
I = np.arcsin(sinI); II = np.arcsin(sinII)
I_data.append(I); II_data.append(II)
UU = U + I - II
LL = lens[recursion].radius + lens[recursion].radius * sinII / np.sin(UU)
LL_data.append(LL); UU_data.append(UU)
#PA校对法
PA1 = L * np.sin(U) / np.cos((I-U)/2)
PA2 = LL * np.sin(UU) / np.cos((II-UU)/2)
if(np.abs(PA1-PA2) > EPS and isParallel == 0):
print("Recursion not right.")
print(f'PA1 = {PA1}, PA2 = {PA2}')
PA_data.append(PA2)
#过渡公式
U = UU
L = LL - lens[recursion].thickness
#舍去最后一个过渡的物方截距
if(recursion != lensnumber - 1):
L_data.append(L)
#print(f'LEN_DATA\n{len(L_data)}, {len(LL_data)}, {len(UU_data)}, {len(I_data)}, {len(II_data)}, {len(PA_data)}')
return L_data, LL_data, UU_data, I_data, II_data, PA_data
def astigmatism_recursion(L_data, LL_data, UU_data, I_data, II_data, PA_data, lens, raytype, isParallel):
"""
函数作用: 求像散, 场曲
输入: 由lens_recursion算得的各参数
输出: 子午场曲(xt), 弧矢场曲(xs), 像散(delta_x)
"""
lensnumber = len(lens)
for recursion in range(lensnumber):
if(raytype == 'd'):
x = PA_data[recursion] ** 2 / (2 * lens[recursion].radius)
#求t、s初值
if(recursion == 0 and isParallel == 0):
t = s = (L_data[0] - x) / np.cos(UU_data[0])
if(recursion == 1 and isParallel == 1):
t = s = (L_data[1] - x) / np.cos(UU_data[1])
temp_tt = (lens[recursion].nd1 * np.cos(II_data[recursion]) - lens[recursion].nd0 * np.cos(I_data[recursion]))\
/ lens[recursion].radius + lens[recursion].nd0 * np.cos(I_data[recursion]) ** 2 / t
tt = 1 / temp_tt * lens[recursion].nd1 * np.cos(II_data[recursion]) ** 2 #t'
temp_ss = (lens[recursion].nd1 * np.cos(II_data[recursion]) - lens[recursion].nd0 * np.cos(I_data[recursion]))\
/ lens[recursion].radius + lens[recursion].nd0 / s
ss = 1 / temp_ss * lens[recursion].nd1 #s'
if(recursion < lensnumber-1):
#过渡公式
xx = PA_data[recursion + 1] ** 2 / (2 * lens[recursion+1].radius) #x(i+1)
D = (lens[recursion].thickness - x + xx) / np.cos(UU_data[recursion+1])
t = tt - D
s = ss - D
if(recursion == lensnumber-1):
#最后一面
lt = tt * np.cos(UU_data[recursion + 1]) + x
ls = ss * np.cos(UU_data[recursion + 1]) + x
xt = lt - LL_data[-1]
xs = ls - LL_data[-1]
delta_x = xt - xs
if(raytype == 'C'):
x = PA_data[recursion] ** 2 / (2 * lens[recursion].radius)
#求t、s初值
if(recursion == 0 and isParallel == 0):
t = s = (L_data[0] - x) / np.cos(UU_data[0])
if(recursion == 1 and isParallel == 1):
t = s = (L_data[1] - x) / np.cos(UU_data[1])
temp_tt = (lens[recursion].nC1 * np.cos(II_data[recursion]) - lens[recursion].nC0 * np.cos(I_data[recursion]))\
/ lens[recursion].radius + lens[recursion].nC0 * np.cos(I_data[recursion]) ** 2 / t
tt = 1 / temp_tt * lens[recursion].nC1 * np.cos(II_data[recursion]) ** 2 #t'
temp_ss = (lens[recursion].nC1 * np.cos(II_data[recursion]) - lens[recursion].nC0 * np.cos(I_data[recursion]))\
/ lens[recursion].radius + lens[recursion].nC0 / s
ss = 1 / temp_ss * lens[recursion].nC1 #s'
if(recursion < lensnumber-1):
#过渡公式
xx = PA_data[recursion + 1] ** 2 / (2 * lens[recursion+1].radius) #x(i+1)
D = (lens[recursion].thickness - x + xx) / np.cos(UU_data[recursion+1])
t = tt - D
s = ss - D
if(recursion == lensnumber-1):
#最后一面
lt = tt * np.cos(UU_data[recursion + 1]) + x
ls = ss * np.cos(UU_data[recursion + 1]) + x
xt = lt - LL_data[recursion]
xs = ls - LL_data[recursion]
delta_x = xt - xs
if(raytype == 'F'):
x = PA_data[recursion] ** 2 / (2 * lens[recursion].radius)
#求t、s初值
if(recursion == 0 and isParallel == 0):
t = s = (L_data[0] - x) / np.cos(UU_data[0])
if(recursion == 1 and isParallel == 1):
t = s = (L_data[1] - x) / np.cos(UU_data[1])
temp_tt = (lens[recursion].nF1 * np.cos(II_data[recursion]) - lens[recursion].nF0 * np.cos(I_data[recursion]))\
/ lens[recursion].radius + lens[recursion].nF0 * np.cos(I_data[recursion]) ** 2 / t
tt = 1 / temp_tt * lens[recursion].nF1 * np.cos(II_data[recursion]) ** 2 #t'
temp_ss = (lens[recursion].nF1 * np.cos(II_data[recursion]) - lens[recursion].nF0 * np.cos(I_data[recursion]))\
/ lens[recursion].radius + lens[recursion].nF0 / s
ss = 1 / temp_ss * lens[recursion].nF1 #s'
if(recursion < lensnumber-1):
#过渡公式
xx = PA_data[recursion + 1] ** 2 / (2 * lens[recursion+1].radius) #x(i+1)
D = (lens[recursion].thickness - x + xx) / np.cos(UU_data[recursion+1])
t = tt - D
s = ss - D
if(recursion == lensnumber-1):
#最后一面
lt = tt * np.cos(UU_data[recursion + 1]) + x
ls = ss * np.cos(UU_data[recursion + 1]) + x
xt = lt - LL_data[recursion]
xs = ls - LL_data[recursion]
delta_x = xt - xs
return xt, xs, delta_x
'''
---数据处理模块---
包含:
data_processing(...), 总输出函数,给定物面、镜面、入瞳、像面信息,向主函数返回所有待求量
Calc_ff(...), 计算计算像方焦距、焦点位置、主点位置
Calc_ll(...), 计算像距、像高、像散、场曲
Calc_lp(...), 计算出瞳距
Calc_sp(...), 计算全孔径和0.7孔径时的实际像距、球差、位置色差
Calc_coma(...), 计算全视场、0.7视场、正负全孔径、正负0.7孔径的彗差、两个视场dFC光实际像高,绝对畸变、相对畸变、倍率色差
'''
def data_processing(object, lens, pupil, image):
'''
共输出38个数据
data1:f', l', lH', lp', y0', y0, 0.707', xt', xs', delta_x', F光近轴像位置, C光近轴像位置;共12个数据
data2:轴上点全孔径的 d光像点, F光像点, C光像点, 球差, 位置色差;共5个数据
data3:轴上点0.707孔径的 d光像点, F光像点, C光像点, 球差, 位置色差;共5个数据
data4:轴上点0孔径位置色差;共1个数据
data5:全视场的 全孔径子午彗差, 0.7孔径子午彗差, d光实际像高, F光实际像高, C光实际像高, 绝对畸变, 相对畸变, 倍率色差;共8个数据
data6:0.707视场的 全孔径子午彗差, 0.7孔径子午彗差, d光实际像高, F光实际像高, C光实际像高, 绝对畸变, 相对畸变, 倍率色差;共8个数据
'''
oo = 1e-20 #定义近轴小量
ff, l_H = Calc_ff(object, lens, pupil, image, oo)
ll, ll_F, ll_C, y0, y0_707, xt, xs, delta_x = Calc_ll(object, lens, pupil, image, oo, l_H, ff)
l_p = Calc_lp(lens, pupil, image, oo)
data1 = [ff, ll, l_H, l_p, y0, y0_707, xt, xs, delta_x, ll_F, ll_C]
data2, data3, data4 = Calc_sp(object, lens, pupil, image, ll, oo)
data5, data6 = Calc_coma(object, lens, pupil, image, oo, y0, y0_707)
return data1, data2, data3, data4, data5, data6
#1、计算像方焦距、焦点位置、主点位置
def Calc_ff(object, lens, pupil, image, oo):
#定义一个近轴无穷远物面
object_ = Surface(0, None, object.nd1, object.nF1, object.nC1, None)
#计算平行入射的近轴光线
L0 = Calc(oo, 0, object_, lens, pupil, image, 'd')[0]
L1 = Calc(oo, 0, object_, lens, pupil, image, 'd')[1]
#计算L0,L1的乘积
L0_product = L1_product = L0F_product = L1F_product = L0C_product = L1C_product = 1
for i in range(len(L0)-1):
L0_product = L0_product * L0[i+1]
for i in range(len(L1)):
L1_product = L1_product * L1[i]
#利用正切计算法计算焦距
ff = L1_product/L0_product
#计算像方焦点和主点位置
l_F = L1[-1]
l_H = ff - l_F
return ff, l_H
这个是正确计算的代码,但是写的比较混乱,下方代码计算结果不正确import tkinter as tk
from tkinter import ttk, messagebox, filedialog
import json
import os
import math
import numpy as np
class OpticalSystem:
def __init__(self):
self.surfaces = [] # 存储光学表面对象
self.entrance_pupil_diameter = None # 入瞳直径 (mm)
self.entrance_pupil_position = None # 入瞳位置 (相对第一个面顶点) (mm)
self.object_infinite = True # 物在无穷远
self.field_angle = None # 半视场角 (度)
self.object_distance = None # 物距 (有限远) (mm)
self.object_height = None # 物高 (有限远) (mm)
self.aperture_angle = None # 孔径角 (有限远) (度)
self.light_type = "d" # 色光类型,默认d光
self.ray_paths = [] # 存储光线路径数据
# 计算结果存储
self.results = {
"focal_length": None, # 焦距 f'
"ideal_image_distance": None, # 理想像距 l'
"actual_image_position": None, # 实际像位置
"image_principal_plane": None, # 像方主面位置 lH'
"exit_pupil_distance": None, # 出瞳距 lp'
"ideal_image_height": None, # 理想像高 y0'
"spherical_aberration": None, # 球差
"longitudinal_chromatic": None, # 位置色差
"tangential_field_curvature": None, # 子午场曲 xt'
"sagittal_field_curvature": None, # 弧矢场曲 xs'
"astigmatism": None, # 像散 Δxts'
"actual_image_height": None, # 实际像高
"relative_distortion": None, # 相对畸变
"absolute_distortion": None, # 绝对畸变
"lateral_chromatic": None, # 倍率色差
"tangential_coma": None # 子午慧差
}
class Surface:
def __init__(self, r=float('inf'), d=0.0, nd=1.0, vd=0.0):
self.r = r # 曲率半径 (mm)
self.d = d # 厚度/间隔 (mm)
self.nd = nd # d光折射率
self.vd = vd # 阿贝数
def to_dict(self):
"""将表面数据转换为字典"""
return {
"r": self.r,
"d": self.d,
"nd": self.nd,
"vd": self.vd
}
@classmethod
def from_dict(cls, data):
"""从字典创建表面对象"""
return cls(
r=data.get('r', float('inf')),
d=data.get('d', 0.0),
nd=data.get('nd', 1.0),
vd=data.get('vd', 0.0)
)
class Calculate:
EPS = 1e-5 # PA校验精度
@staticmethod
def recursion(L, U, H, surfaces, P, light_type='d'):
"""
光线追迹递归计算
:param L: 物方截距
:param U: 物方孔径角
:param H: 光线高度
:param surfaces: 光学表面列表
:param P: PA校验标志
:param light_type: 色光类型
:return: 光线路径数据
"""
# 数据初始化
snumber = len(surfaces)
sinI = 0; sinII = 0
I = 0; II = 0
UU = 0
LL = 0
L_data = [L]
LL_data = []
UU_data = [U]
I_data = []
II_data = []
PA_data = []
# 前一个面的折射率(初始为空气)
n_prev = 1.0
for i, surf in enumerate(surfaces):
# 根据色光类型选择折射率
if light_type == 'd':
n_current = surf.nd
elif light_type == 'f':
# F光折射率计算 (nd + (nf - nc)/(vd * (nf - nc)/ (nd - 1)))
n_current = surf.nd + 1/(surf.vd * (1.0/(surf.nd - 1.0)))
elif light_type == 'c':
# C光折射率计算
n_current = surf.nd - 1/(surf.vd * (1.0/(surf.nd - 1.0)))
else:
n_current = surf.nd # 默认为d光
# 计算入射角
if i == 0 and P:
# 第一面且平行光入射
sinI = H / surf.r if surf.r != 0 else 0
sinII = (n_prev / n_current) * sinI
else:
# 普通情况
sinI = (L_data[-1] - surf.r) / surf.r * np.sin(U) if surf.r != 0 else 0
sinII = (n_prev / n_current) * sinI
# 确保sin值在合理范围内
sinI = np.clip(sinI, -1.0, 1.0)
sinII = np.clip(sinII, -1.0, 1.0)
I = np.arcsin(sinI)
II = np.arcsin(sinII)
I_data.append(I)
II_data.append(II)
# 计算折射后角度
UU = U + I - II
UU_data.append(UU)
# 计算折射后截距
if abs(sinII) > 1e-10: # 避免除以0
LL = surf.r + surf.r * sinII / np.sin(UU)
else:
LL = surf.r # 垂直入射
LL_data.append(LL)
# PA校对法
if abs(U) > 1e-5 and abs(UU) > 1e-5: # 避免除以0
PA1 = L_data[-1] * np.sin(U) / np.cos((I - U)/2)
PA2 = LL * np.sin(UU) / np.cos((II - UU)/2)
if abs(PA1 - PA2) > Calculate.EPS and P == 0:
print(f"Warning: PA校验不通过! 表面 {i+1}: PA1={PA1:.6f}, PA2={PA2:.6f}")
PA_data.append(PA2)
else:
PA_data.append(0.0) # 无法计算PA值
# 过渡到下一个面
U = UU
if i < snumber - 1: # 最后一个面不需要过渡
L = LL - surf.d
L_data.append(L)
# 更新折射率为当前面的折射率
n_prev = n_current
return L_data, LL_data, UU_data, I_data, II_data, PA_data
@staticmethod
def calc(system, K_eta=0.0, K_W=0.0):
"""
计算某一孔径、视场光线在逐个面的(L,U)坐标数据
:param system: OpticalSystem对象
:param K_eta: 孔径取点系数 (0-1)
:param K_W: 视场取点系数 (0-1)
:return: 光线路径数据
"""
# 初始化(L,U)坐标列表
L0 = []
U_data = []
I0 = []
I1 = []
PA = []
A = system.entrance_pupil_diameter / 2 # 入瞳半径
eta = K_eta * A # 实际像高系数
# 物在无穷远且光线与光轴平行时的情形
if system.object_infinite and K_W == 0:
h0 = K_eta * A
L0, LL, U_data, I0, I1, PA = Calculate.recursion(
0, 0, h0, system.surfaces, True, system.light_type
)
else:
# 物在无穷远但视场不为0
if system.object_infinite:
W = K_W * system.field_angle * np.pi / 180 # 转换为弧度
U_angle = -W
L_start = system.entrance_pupil_position + eta / np.tan(U_angle)
else:
# 有限远物距
Lp = system.entrance_pupil_position # 第一面到入瞳的距离
y = K_W * system.object_height # 物方线视场
if K_W == 0: # 轴上点
dist = np.sqrt(Lp**2 + A**2)
U_angle = np.arcsin(K_eta * A / dist)
L_start = -system.object_distance
else: # 轴外点
U_angle = np.arctan((y - eta) / (Lp + system.object_distance))
L_start = Lp + eta / np.tan(U_angle)
L0, LL, U_data, I0, I1, PA = Calculate.recursion(
L_start, U_angle, 0, system.surfaces, False, system.light_type
)
return L0, LL, U_data, I0, I1, PA
@staticmethod
def calculate_focal_length(system):
"""计算焦距"""
# 计算近轴光线
_, LL, U_data, _, _, _ = Calculate.calc(system, K_eta=0.0001, K_W=0.0)
last_U = U_data[-1]
last_L = LL[-1]
# 焦距 f' = h / tan(U')
if abs(last_U) > 1e-5:
focal_length = system.entrance_pupil_diameter / 2 / np.tan(last_U)
return focal_length
return float('inf') # 平行光出射
@staticmethod
def calculate_ideal_image_distance(system):
"""计算理想像距"""
# 计算主光线
L0, LL, U_data, _, _, _ = Calculate.calc(system, K_eta=0.0, K_W=0.0)
last_L = LL[-1]
return last_L
@staticmethod
def calculate_spherical_aberration(system):
"""计算球差"""
# 计算边缘光线
_, LL_edge, U_edge, _, _, _ = Calculate.calc(system, K_eta=1.0, K_W=0.0)
# 计算近轴光线
_, LL_paraxial, U_paraxial, _, _, _ = Calculate.calc(system, K_eta=0.0001, K_W=0.0)
# 球差 δL' = L'边缘 - L'近轴
return LL_edge[-1] - LL_paraxial[-1]请说说第二段代码运算逻辑哪里有问题
最新发布