|  | 
|  | 1 | +from numpy import * | 
|  | 2 | +from pylab import * | 
|  | 3 | +import timeit | 
|  | 4 | + | 
|  | 5 | + | 
|  | 6 | +class DBRMirror: | 
|  | 7 | + def __init__(self, ni, nbr): # ni - refractive indices as array, nbr - number of layers (must be even!) | 
|  | 8 | + if nbr % 2 != 0: print("Number must be even! Program will give wrong results") | 
|  | 9 | + self.nbr = nbr | 
|  | 10 | + self.n = array(concatenate(([1], tile(ni, (self.nbr-2)//2), [ni[0]]))) | 
|  | 11 | + | 
|  | 12 | + def reflection_coeff(self, incidentwavelength, braggwavelength): | 
|  | 13 | + nm = (self.n[:-1] / self.n[1:]) | 
|  | 14 | + R = array([]) | 
|  | 15 | + for value in braggwavelength: | 
|  | 16 | + d = array(concatenate(([0], 0.25 * value/self.n[1:-1], [0]))) | 
|  | 17 | + phi = 2j*pi*self.n[:-1]*d[:-1]/incidentwavelength | 
|  | 18 | + transition_matrix = array([[0.5*(1+nm)*exp(-phi), 0.5*(1-nm)*exp(phi)], | 
|  | 19 | + [0.5*(1-nm)*exp(-phi), 0.5*(1+nm)*exp(phi)]]) | 
|  | 20 | + | 
|  | 21 | + general_transition_matrix = 1 | 
|  | 22 | + for k in range(self.nbr - 2,-1,-1): | 
|  | 23 | + general_transition_matrix = dot(general_transition_matrix, transition_matrix[:, :, k]) | 
|  | 24 | + | 
|  | 25 | + # print(general_transition_matrix) | 
|  | 26 | + r = - general_transition_matrix.item(2) / general_transition_matrix.item(3) | 
|  | 27 | + R = append(R, abs(r)**2) | 
|  | 28 | + | 
|  | 29 | + # t = linalg.det(general_transition_matrix) / general_transition_matrix.item(3) | 
|  | 30 | + # T = n.item(N-1)/n.item(0) * abs(t)**2 | 
|  | 31 | + # | 
|  | 32 | + # print(R) | 
|  | 33 | + # print(T) | 
|  | 34 | + # print(R+T) | 
|  | 35 | + return R | 
|  | 36 | + | 
|  | 37 | + | 
|  | 38 | +start = timeit.default_timer() | 
|  | 39 | +wave_length = 0.98 * 10 ** (-6) | 
|  | 40 | +NewMirror = DBRMirror([3.0, 3.5], 42) | 
|  | 41 | +braggwavelength = array(linspace(0.68,1.28,100000)*10**(-6)) | 
|  | 42 | +R = NewMirror.reflection_coeff(wave_length,braggwavelength) | 
|  | 43 | + | 
|  | 44 | +stop = timeit.default_timer() | 
|  | 45 | +print(stop - start) | 
|  | 46 | + | 
|  | 47 | +plot(braggwavelength, R) | 
|  | 48 | +xlabel("$\lambda$") | 
|  | 49 | +ylabel("R") | 
|  | 50 | +show() | 
0 commit comments