1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

from math import *
from sympy import *
import matplotlib.pyplot as plt
import numpy as np

#理想光线追迹
def idealight(r, d, n1, n2, l1=-1000):
    u11 = 0.0000001

    i11 = (l1 - r[1]) * u11 / r[1]
    i12 = 1 * i11 / n1
    u12 = u11 + i11 - i12
    l12 = r[1] * (1 + i12 / u12)

    l21 = l12 - d[1]
    u21 = u12

    i21 = (l21 - r[2]) * u21 / r[2]
    i22 = n1 * i21 / n2
    u22 = u21 + i21 - i22
    l22 = r[2] * (u22 + i22) / u22

    l31 = l22 - d[2]
    u31 = u22

    i31 = (l31 - r[3]) * u31 / r[3]
    i32 = n2 * i31
    u32 = u31 + i31 - i32
    l32 = r[3] * (1 + i32 / u32)
    return l32

#实际光线追迹
def reallight(u11, r, d, n1, n2, l1=-1000):
    sini11 = (l1 - r[1]) * sin(u11) / r[1]
    sini12 = 1 * sini11 / n1
    u12 = u11 + asin(sini11) - asin(sini12)
    l12 = r[1] * (1 + sini12 / sin(u12))

    l21 = l12 - d[1]
    u21 = u12

    sini21 = (l21 - r[2]) * sin(u21) / r[2]

    sini22 = n1 * sini21 / n2
    u22 = u21 + asin(sini21) - asin(sini22)

    l22 = r[2] * (sin(u22) + sini22) / sin(u22)

    l31 = l22 - d[2]
    u31 = u22

    sini31 = (l31 - r[3]) * sin(u31) / r[3]
    sini32 = n2 * sini31
    u32 = u31 + asin(sini31) - asin(sini32)
    l32 = r[3] * (1 + sini32 / sin(u32))
    return l32

#初始条件：物距，透镜半径，透镜距离,对应
#初始0是因为python从0索引
l = -1000
r = [0, 150, -120, -200]
d = [0, 20, 20]
#折射率 C,d,F
n1 = [1.4,1.55,1.6]
n2 = [1.9, 2.1, 2.2]

#求解光阑参数：以第一个透镜为光阑
h = Symbol('h')
maxh = solve( -(r[1]**2-h**2)**0.5+r[1]-r[2]-(r[2]**2-h**2)**0.5-d[1],h)
maxu = int(10000*atan(maxh[1] / (-l + d[1] - (r[1] - d[1] - (r[1] ** 2 - maxh[1] ** 2) ** 0.5))))


#初始化序列
divi=500
uin = np.linspace(0, maxu/10000, divi)


#初始化光线
deltac = [0]*divi
deltad = [0]*divi
deltaf = [0]*divi
#列表
for i in range(divi):
    deltac[i] = [idealight(r, d, n1[0], n2[0], l) - reallight(uin[i], r, d, n1[0], n2[0], l)]
    deltad[i] = [idealight(r, d, n1[1], n2[1], l) - reallight(uin[i], r, d, n1[1], n2[1], l)]
    deltaf[i] = [idealight(r, d, n1[2], n2[2], l) - reallight(uin[i], r, d, n1[2], n2[2], l)]
#画图
pc,=plt.plot(uin, deltac, label='C',color='r')
pd,=plt.plot(uin, deltad, label='d',color='gold')
pf,=plt.plot(uin, deltaf, label='F',color='b')
plt.xlabel('U/rad')
plt.ylabel('delta/mm')
plt.legend(handles = [pc,pd,pf],loc=3)
plt.show()