コードの書き方
import math
import cmath
import numpy as np
import matplotlib.pyplot as plt
def multi_reflection(n1, n2, n3, theta1_deg, d, wavelength):
#calculation
theta1_rad = math.radians(theta1_deg)
#屈折角の算出
theta2_rad = math.asin(n1/n2*math.sin(theta1_rad))
theta2_deg = math.degrees(theta2_rad)
theta3_rad = math.asin(n2/n3*math.sin(theta2_rad))
theta3_deg = math.degrees(theta2_rad)
#位相差
delta = 4 * cmath.pi * n2 * d * math.cos(theta2_rad)/ wavelength
r12s = (n1 * math.cos(theta1_rad) - n2 * math.cos(theta2_rad)) / (n1 * math.cos(theta1_rad) + n2 * math.cos(theta2_rad))
#t12s = 2 * n1 * math.cos(theta1_rad) / (n1 * math.cos(theta1_rad) + n2 * math.cos(theta2_rad))
r23s = (n2 * math.cos(theta2_rad) - n3 * math.cos(theta3_rad)) / (n2 * math.cos(theta2_rad) + n3 * math.cos(theta3_rad))
r12p = (n2 * math.cos(theta1_rad) - n1 * math.cos(theta2_rad)) / (n2 * math.cos(theta1_rad) + n1 * math.cos(theta2_rad))
#t12p = 2 * n1 * math.cos(theta1_rad) / (n2 * math.cos(theta1_rad) + n1 * math.cos(theta2_rad))
r23p = (n3 * math.cos(theta2_rad) - n2 * math.cos(theta3_rad)) / (n3 * math.cos(theta2_rad) + n2 * math.cos(theta3_rad))
delta_j = complex(0, delta)
Rs = (r12s + r23s * cmath.exp(delta_j))/(1 + r23s * r12s * cmath.exp(delta_j))
Rs = abs(Rs) * abs(Rs)
Rp = (r12p + r23p * cmath.exp(delta_j))/(1 + r23p * r12p * cmath.exp(delta_j))
Rp = abs(Rp) * abs(Rp)
if wavelength == 1000:
print("-------input parameters-------")
print("theta1={0} deg".format(theta1_deg))
#print("n1={0}, n2={2}".format(n1, n2))
print("-------output@1000nm-------")
print("theta2_deg={0} deg".format(theta2_deg)) #計算チェック完了
print("Rp={0}".format(Rp)) #水の反射率は再現
print("Rs={0}".format(Rs)) #水の反射率は再現
return [Rp, Rs]
#parameters
#wavelength = 1000 #nm
wavelength = np.arange(900, 1500, 1)
n1 = 1
theta1_deg = 0
n2 = 1.5
d = 210 #nm
n3 = 3.6
Rs = []
Rp = []
for s in wavelength:
Rs.append(multi_reflection(n1, n2, n3, theta1_deg, d, s)[1])
Rp.append(multi_reflection(n1, n2, n3, theta1_deg, d, s)[0])
plt.plot(wavelength, Rp)
plt.ylabel('R')
plt.ylim(0, 1)
plt.xlabel('wavelength nm')
plt.show()
