虽然可能但是matlab、mathmatica更好用,少掉点头发还是用了python。
一维也算是比较简单了,我看很多哥哥们都喜欢收费,我想了想肯定是免费。
这里设置了一个函数,Film(),在前面自己改一下波长lam还有一些个initial condition
import numpy as np
from matplotlib import pyplot as plt
class Film():
def __init__(self, n_list, d_list, labels=None, nd=True, lambda0=550):
self.n_list = n_list
self.nd = nd
self.lambda0 = lambda0
self.d_list = d_list
self.labels = labels
def __part_fit(self, rect_list, theta_cos, n_list, j):
rect_list = np.array(rect_list, dtype=np.complex128)
rect_list = rect_list.transpose(3, 0, 1, 2)
n_theta_cos = n_list[j][-1] * theta_cos
Y = []
for i in range(rect_list.shape[0]):
rect = np.eye(2)
for k in range(rect_list[i].shape[0]):
rect = np.dot(rect, rect_list[i, k])
rect = np.dot(rect, np.array([[1], [n_theta_cos[i]]]))
Y.append(rect)
Y = np.array(Y)
Y = Y[:, 1, 0] / Y[:, 0, 0]
return Y
def fit(self, wavelength, theta=0):
if type(wavelength) == float or type(wavelength) == int:
wavelength = [wavelength]
self.wavelength = np.array(wavelength)
n_list = []
n_lambda0 = []
for i in self.n_list:
n_list.append([])
n_lambda0.append([])
for j in i:
if type(j).__name__ != 'function':
n_list[-1].append(np.ones_like(wavelength) * j)
n_lambda0[-1].append(j)
else:
n_list[-1].append(j(wavelength))
n_lambda0[-1].append(j(np.array([self.lambda0]))[0])
d_list = []
if self.nd:
for i in range(len(n_list)):
d_list.append([])
for j in range(len(self.d_list[i])):
d_list[-1].append(self.d_list[i][j] * n_list[i][j+1] / n_lambda0[i][j+1])
else:
for i in range(len(n_list)):
d_list.append([])
for j in range(len(self.d_list[i])):
d_list[-1].append(self.d_list[i][j] * n_list[i][j+1])
self.R_p = []
self.R_s = []
theta_cos_0 = np.cos(theta)
for j in range(len(n_list)):
rect_list = [[], []]
theta_copy = theta
for i in range(1, len(d_list[j])+1):
# 折射定律
theta_copy = np.arcsin(n_list[j][i-1] * np.sin(theta_copy) / n_list[j][i])
theta_cos = np.cos(theta_copy)
delta = 2 * np.pi * d_list[j][i-1] / self.wavelength * theta_cos
# p光
rect_single_p = np.array([[np.cos(delta), 1j * np.sin(delta) * theta_cos / n_list[j][i]],
[1j * n_list[j][i] * np.sin(delta) / theta_cos, np.cos(delta)]])
rect_list[0].append(rect_single_p)
# s光
rect_single_s = np.array([[np.cos(delta), 1j * np.sin(delta) / (n_list[j][i] * theta_cos)],
[1j * n_list[j][i] * np.sin(delta) * theta_cos, np.cos(delta)]])
rect_list[1].append(rect_single_s)
theta_copy = np.arcsin(n_list[j][-2] * np.sin(theta_copy) / n_list[j][-1])
theta_cos = np.cos(theta_copy)
Y_p = self.__part_fit(rect_list[0], 1 / theta_cos, n_list, j)
Y_s = self.__part_fit(rect_list[1], theta_cos, n_list, j)
self.R_p.append(abs((n_list[j][0] / theta_cos_0 - Y_p) / (n_list[j][0] / theta_cos_0 + Y_p))**2)
self.R_s.append(abs((n_list[j][0] * theta_cos_0 - Y_s) / (n_list[j][0] * theta_cos_0 + Y_s))**2)
self.R_p = np.array(self.R_p)
self.T_p = 1 - self.R_p
self.R_s = np.array(self.R_s)
self.T_s = 1 - self.R_s
self.R = 0.5 * (self.R_p + self.R_s)
self.T = 1 - self.R
def __part_process(self, data, title):
if not self.labels:
for i in data:
plt.plot(self.wavelength, i)
else:
for i in range(len(data)):
plt.plot(self.wavelength, data[i], label=self.labels[i])
plt.legend()
plt.title(title)
def __single_image(self, percentage, img_type, R, T):
if percentage:
R *= 100
T *= 100
if img_type == 'r' or img_type == 'R':
self.__part_process(R, 'R')
elif img_type == 't' or img_type == 'T':
self.__part_process(T, 'T')
def __multi_image(self, percentage, img_type, R, T):
R = np.array(R)
T = np.array(T)
light = ['ave', 'p', 's']
if percentage:
R *= 100
T *= 100
if img_type == 'r' or img_type == 'R':
for i in range(len(R)):
plt.plot(self.wavelength, R[i][0], label=light[i])
plt.title('R')
elif img_type == 't' or img_type == 'T':
for i in range(len(T)):
plt.plot(self.wavelength, T[i][0], label=light[i])
plt.title('T')
plt.legend()
def __display(self, percentage, img_type, light_type):
if light_type == 'ave':
self.__single_image(percentage, img_type, self.R, self.T)
elif light_type == 'p':
self.__single_image(percentage, img_type, self.R_p, self.T_p)
elif light_type == 's':
self.__single_image(percentage, img_type, self.R_s, self.T_s)
elif light_type == 'all':
self.__multi_image(percentage, img_type,
[self.R, self.R_p, self.R_s], [self.T, self.T_p, self.T_s])
def show(self, percentage=False, img_type='R', light_type='ave'):
self.__display(percentage, img_type, light_type)
plt.show()
def savefig(self, filename, percentage=False, img_type='R', light_type='ave'):
self.__display(percentage, img_type, light_type)
plt.savefig(filename)
具体咋用来俩栗子:
单层
#计算单层膜
wavelength = np.arange(400,900,0.1)#arange函数类比于python的内置函数,通过开始值、终值和步长来创建等差数列。注意不包含终值。
# arange()语法是括号内四个值分别是:开始值start、终值stop、步长step、数组类型dtype。
n_list = [[1, 1.38, 1.52], [1, 1.38, 1.75], [1, 1.38, 1.9]]
d_list = [[130], [130], [130]]
labels = ['1.52', '1.75', '1.90']
film = Film(n_list, d_list, labels)
film.fit(wavelength)
film.show()
双层/多层
#双层薄膜
wavelength = np.arange(400, 800, 0.1)
n_list = [[1, 1.38, 1.7, 1.52], [1, 1.38, 2.3, 1.52]]
d_list = [[135, 135], [174.6, 28.3]]
labels = ['a','b']
film = Film(n_list, d_list, labels)
film.fit(wavelength)
film.show()
注意区分光的偏振
#绘制平均光、p光和s光
wavelength = np.arange(400, 800, 0.1)
n_list = [[1, 1.38, 1.7, 1.52]]
d_list = [[137.5, 275]]
film = Film(n_list, d_list)
film.fit(wavelength, theta=np.pi/4)
film.show(light_type='all')
之前认为折射率恒定,但是像银Ag折射率随着入射波的波长变化,这个就得修改一下了
#折射率为变量,是函数
def F(x):
return 1 + np.sin(0.002 * x)
wavelength = np.arange(380, 780, 0.1)
n_list = [[1, F, 1.52]]
d_list = [[137.5]]
film = Film(n_list, d_list)
film.fit(wavelength)
film.show(percentage=True)
在这个基础上能玩好多花,但是一维的就到此结束吧。二维的过两天出一个。