コードの書き方
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()
お礼
ありがとうございます。自分でも計算はしたのですが不安だったので質問させていただきました。