物理类-一维光子晶体反射率以及透射率

虽然可能但是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)

在这个基础上能玩好多花，但是一维的就到此结束吧。二维的过两天出一个。