# モルフォチョウの多層膜を再現して各波長に対する反射率を計算するプログラム # # モルフォチョウの翅を 11層の多層膜をと見立て、 # クチクラ層 (屈折率1.6, 膜幅65nm), 空気層 (屈折率1.0, 膜幅130nm) として計算 # # 参考: 構造色研究会「エクセルで多層膜干渉!! 」 # http://mph.fbs.osaka-u.ac.jp/~ssc/excelde/excelde.html # # プログラム作者: @tsujimotter

multilayer.rb

require 'complex'

# モルフォチョウの多層膜を再現して各波長に対する反射率を計算するプログラム
# 
# モルフォチョウの翅を 11層の多層膜をと見立て、
# クチクラ層 (屈折率1.6, 膜幅65nm), 空気層 (屈折率1.0, 膜幅130nm) として計算
# 
# 参考: 構造色研究会「エクセルで多層膜干渉!! 」 
# http://mph.fbs.osaka-u.ac.jp/~ssc/excelde/excelde.html
# 
# プログラム作者: @tsujimotter
# 

# S偏光の反射係数
def rs(n1, n2, cos1, cos2)
  (n1*cos1 - n2*cos2) / (n1*cos1 + n2*cos2)
end
# S偏光の透過係数
def ts(n1, n2, cos1, cos2)
  2*n1*cos1 / (n1*cos1 + n2*cos2)
end

# P偏光の反射係数
def rp(n1, n2, cos1, cos2)
  (n2*cos1 - n1*cos2) / (n2*cos1 + n1*cos2)
end
# P偏光の透過係数
def tp(n1, n2, cos1, cos2)
  2*n1*cos1 / (n2*cos1 + n1*cos2)
end



# 入射角 (cos = 1.0, 垂直入射)
cos_in = 1.0

# 各層の屈折率
# 入射層を0番目 (入射層) とし、最後が透過した光が出ていく層 
n   = [1.0, 1.6, 1.0, 1.6, 1.0, 1.6, 1.0, 1.6, 1.0, 1.6, 1.0, 1.6, 1.0]

# 層の数 (最初と最後は空気層)
size = n.length - 2

# 各層の厚さ
# 入射層と透過光が出ていく層は 0 
d = [0, 65e-9, 130e-9, 65e-9, 130e-9, 65e-9, 130e-9, 65e-9, 130e-9, 65e-9, 130e-9, 65e-9, 0]

# スネルの法則に基づいて各層の余弦を計算する関数
def cosine(n0, cos0, n)
  c = 0.0
  if n0 == n
    c = cos0
  else
    c = Math::sqrt(1-(1-cos0*cos0)*(n0*n0)/(n*n))
  end
  c
end
  
# 300 nm から 900 nm まで計算
for l in 300..900

  # 波長 (単位メートル)
  lambda = 1e-9 * l
  # 波数を計算
  k = 2*Math::PI/lambda

  # 最終層の余弦を計算  
  cos_c = cosine(n[0], cos_in, n[size])
  # 最終層から出ていく余弦を計算  
  cos_out = cosine(n[0], cos_in, n[size+1])

  # 最終から戻っていく反射係数を計算
  rs = rs(n[size], n[size+1], cos_c, cos_out)
  rp = rp(n[size], n[size+1], cos_c, cos_out)
  
  for i in 1..size
    # 現在の層の上層の余弦を計算
    cos_up = cosine(n[0], cos_in, n[i-1])
    # 現在の層の余弦を計算
    cos_c = cosine(n[0], cos_in, n[i])
        
    # S偏光の反射係数を計算
    rs_c = rs(n[i-1], n[i], cos_up, cos_c)
    # S偏光の透過係数を計算
    ts_c = ts(n[i-1], n[i], cos_up, cos_c)  
    # P偏光の反射係数を計算
    rp_c = rp(n[i-1], n[i], cos_up, cos_c)
    # P偏光の透過係数を計算
    tp_c = tp(n[i-1], n[i], cos_up, cos_c)  
  
    #　位相差を計算
    delta = 2*n[i]*k*d[i]*cos_c
  
    # 現在の層から出ていく反射係数を計算
    rs = (rs_c + rs * Math::exp(Complex::I*delta)) / (1 + rs*rs_c * Math::exp(Complex::I*delta))
    rp = (rp_c + rp * Math::exp(Complex::I*delta)) / (1 + rp*rp_c * Math::exp(Complex::I*delta))
  
  end

  # 計算結果を出力 (波長[ナノメートル] S偏光反射率 P偏光反射率)
  puts (lambda*1e9).to_s + " " +(rs.abs * rs.abs).to_s + " " +(rp.abs * rp.abs).to_s

end