[Python]楕円内での光の反射を考え、その軌跡を描くプログラムを作る

はじめに

今回の楕円内での光の反射の軌跡を描くプログラムを作るにあたって、ただそのプログラムを作成する以上に、いかに数学的に、楕円内で反射する光の軌道を計算し、プログラムするのかを重視している.

Githubから来た人には、work2の改訂版(コードを短くし、分かりやすくした)の紹介も兼ねているので、Githubに公開されているコードよりもこっちの記事を見てほしい.

今回使う公式・数式

• 楕円の公式 $\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$

• 楕円上のある座標(m,n)における接線の方程式 $\frac{mx}{a^2}+\frac{ny}{b^2}=1$
  変形させると $y=-\frac{b^2m}{a^2n}x + \frac{b^2}{n}$

• 楕円と直線の交点の座標

  直線：$y=cx+d$

  $A=(\frac{1}{a^2}+\frac{c^2}{b^2})$
  $B=2(\frac{cd}{b^2})$
  $C=(\frac{d^2}{b^2}-1)$
  $D=B^2-4AC$

  のように、直線と$A,B,C,D$を定義したときに

  $x1=\frac{-B+\sqrt{D}}{2A}$
  $x2=\frac{-B-\sqrt{D}}{2A}$

  交点のx座標である$(x1,x2)$が求まる.

プログラムの作成

使うライブラリはnumpyとmatplotlib

• まず始めにライブラリをimport

import numpy as np
import matplotlib.pyplot as plt

next_position関数を定義

• next_position関数の引数は(今向いている方角(傾き)、今いる座標、$a$、$b$)
• 返り値は、次に到達する座標

def next_position(now_direc ,now_position ,a ,b):
  new_position = [None,None]
  c = now_direc
  d = -1 * now_direc * now_position[0] + now_position[1]

  A = (1/a**2 + c**2/b**2)
  B = 2 * c * d / b**2
  C = (d**2 / b**2) -1
  D = B**2 - 4 * A * C

  x1 = (-B + D**(1/2)) / (2*A) #直線と楕円の交点のx座標(x1,x2)
  x2 = (-B - D**(1/2)) / (2*A)
  #今、x1とx2のどちらの座標にいるか分からないので条件分岐させる
  if round(now_position[0],8) == round(x2,8):
  #コンピューター上では、5が5.000001になったりして正しく計算できないのでroundで少数第8桁まで丸め込む
    new_position[0] = x1
    s = (0 - now_position[1])/(a - now_position[0]) #今いる座標と点(x,y) = (a,0)を通る直線の傾き

    if now_direc >= s:
      new_position[1] = (b**2 - (x1**2 * b**2 / a**2)) ** (1/2)
    else:
      new_position[1] = -(b**2 - (x1**2 * b**2 / a**2)) ** (1/2)
  
  else:
    new_position[0] = x2
    s = (0 - now_position[1])/(-a - now_position[0])#今いる座標と点(x,y) = (-a,0)を通る直線の傾き

    if now_direc >= s:
      new_position[1] = -(b**2 - (x2**2 * b**2 / a**2)) ** (1/2)
    else:
      new_position[1] = (b**2 - (x2**2 * b**2 / a**2)) ** (1/2)

  return new_position#次に到達する座標

next_position関数の説明

if round(now_position[0],8) == round(x2,8):
  new_position[0] = x1

最初の分岐で、今、x1とx2のどちらの座標にいるか判断し、今いない方のx座標(次に到達する座標のx座標)をnew_position[0]に代入.

今回の例では上の条件分岐通り、今いる$x$座標は、$x2$とする$(x1 > x2)$

s = (0 - now_position[1])/(a - now_position[0])
if now_direc >= s:
  new_position[1] = (b**2 - (x1**2 * b**2 / a**2)) ** (1/2)
else:
  new_position[1] = -(b**2 - (x1**2 * b**2 / a**2)) ** (1/2)

円において、$y$の式は、$y = \pm{\sqrt{(b^2-\frac{b^2x^2}{a^2})}}$の2種類あるため、ここでも条件分岐させる必要がある.

下の図では、黄点は$x1$に対応する座標の例、青点は$x2$に対応する座標の例としている(この時点では、まだ黄点の座標は分からないが$x1 > x2$ということは分かっている).

変数$s$:$\frac{0-今いるy座標}{a-今いるx座標}$　　は青点から赤点に直線(赤線)を引いたときの傾きになっている.
そして、現在向いている向き(傾き:now_direc)が青線もしくは緑線となる.

現在いる$x$座標が今回の例のような$x2$だった場合、青線もしくは緑線の傾きが赤線よりが大きい(now_direc$>=s$)と、必ず次の到達点の$y$座標は0以上になる.
よって、$y = \sqrt{(b^2-\frac{b^2x^2}{a^2})}$に次に到達する$x$座標を代入すると、次に到達する$y$座標が出力される.

現在いる座標が$x1$だった場合は、これと逆になる.

---

reflect_direc関数を定義

• 引数は(次に到達する座標、現在向いている方角(傾き)、a、b)
• 返り値は、次に到達する座標に到達したときに反射する向き(傾き)

def reflect_direc(new_position,now_direc,a,b):

  #tan_degは次到達する座標の接線の傾きの角度(-90<=tan_deg<=90)
  if new_position[1] == 0:
    tan_deg = 90
  else:
    deg = -b**2 * new_position[0]/(a**2 * new_position[1])
    tan_deg = np.degrees(np.arctan(deg))

  direc_deg = np.degrees(np.arctan(now_direc))#now_direcを角度に変換(-90<=direc_deg<=90)
  theta = 2*tan_deg - direc_deg

  return_direc = np.tan(np.radians(theta))
  
  return return_direc

reflect_direc関数の説明

 if new_position[0] == 0:
    tan_deg = 90        
 else:
    deg = -b**2 * new_position[0]/(a**2 * new_position[1])
    tan_deg = np.degrees(np.arctan(deg))

  direc_deg = np.degrees(np.arctan(now_direc))

楕円上のある座標における接線の傾きは、上の公式、数式一覧で紹介している通り$-\frac{b^2m}{a^2n}$だが、分母に0を取ることができない.

new_position[1]が$n$にあたるため、new_position[1] == 0の場合、すなわち次に到達する座標が$(a,0)$か$(-a,0)$のとき、あらかじめtan_deg = 90と設定しておく必要がある.

direc_degで現在向ている方角(傾き:now_direc)を角度に変換(-90<=direc_deg<=90)

  theta = 2*tan_deg - direc_deg

  return_direc = np.tan(np.radians(theta))
  
  return return_direc

thetaは反射した後の光の進む方向($x$軸に対する角度)となっている.

注意
角度は-90~90で定義されています.

上の図で、反射板(緑線)に対する入射光(赤線)の角度(direc_deg - tan_deg)は、角度紫(反射板に対する入射角の角度なので、上の例の場合は負)となっていることが分かる.

よって、「反射板(緑線)に対する反射光の反射角度」は、「反射板(緑線)に対する入射光(赤線)の入射角度(角度赤)をプラスマイナス反転したもの」で、tan_deg-direc_degとなる.

その「反射板(緑線)に対する反射光の反射角度(tan_deg-direc_deg)」を角度緑(x軸に対する角度なので、今回の例の場合は負)だけ時計回り方向に回転させる((tan_deg - direc_deg) + tan_deg)と、角度黒(x軸に対する反射光の角度(2*tan_deg - direc_deg))が出てくることが分かる.

つまり上の図のように、まず反射板(緑線)に対する反射光(黒線)の反射角度を求めて、そしてその反射板を任意の角度(橙角度)だけ回転するのを想像すればいい.

• 上の例では、楕円内の第3象限の部分で、光が反射しているときの例を紹介しているが、その場合は反射板に対する反射光の反射角度にtan_deg(橙角度)を足す.

• 第1象限での反射の場合は、180 + tan_degだけ足す(傾き(tan(角度))を求めたいので180は消すことができる).

• 第2象限での反射の場合は、180-tan_degだけ引く(傾き(tan(角度))を求めたいので180は消すことができる).

• 第4象限での反射の場合は、tan_degだけ足す.

よって、反射光の反射角(x軸に対する角度)は2*tan_deg - direc_degとなる.

最後に、return_direcで、角度を傾きに変換してそれを関数の返り値とする.

最後にプロットする

入力値は、最初に光を放つ角度、何回反射させるか、楕円のx方向の大きさ、y方向の大きさ

now_direc = np.tan(np.radians(int(input('打ち出す角度を設定してください(0から180):'))))
h = int(input('何回光を反射させるか決めてください：'))
a = int(input('楕円のx軸方向の長さを設定してください：'))
b = int(input('楕円のy軸方向の長さを設定してください：'))

first_posi = [0,-b]#光を放つ初期位置を設定.楕円の真ん中下から放つ.
next = next_position(now_direc,first_posi,a,b)
direc = reflect_direc(next,now_direc,a,b)#次の場所に向かうための方向を計算

posi = np.array([[0],[-b]])
for _ in range(0,h): #h回反射させます。
  posi = np.append(posi,np.array([[next[0]],[next[1]]]),axis = 1)
  next = next_position(direc,next,a,b)
  direc = reflect_direc(next,direc,a,b)

t = np.linspace(-a,a,a*100)#楕円の作成
y = (b**2 - (t**2 * b**2 / a**2)) ** (1/2)
y1 = -1 * (b**2 - (t**2 * b**2 / a**2)) ** (1/2)

plt.plot(t,y,c = 'b')#楕円をプロット
plt.plot(t,y1,c = 'b')#反射する光の軌跡をプロット

plt.grid()
plt.plot(posi[0,:],posi[1,:])
plt.axis('equal')
plt.show()

結果

80度、60回、$a = 1,b = 1$だと

となる。
$a = 3,b = 2$に変えてみると、上と全然違う形になる.

$a = 4,b = 9$の場合も形が全然違う.

終わりに

放つ角度が同じでも、aとbを変えるだけでここまで出力される軌跡が違うのは、楕円の面白いところですね.

また今度このコードを使って、進化版ビリヤードのようなものや、放つ角度を少しずつ変えていった時に、反射する光の軌跡がどのように変化していくかを見れるプログラムも書きたいと思います.

この記事が、楕円内に放たれた光の反射をどのように考え、プログラムするかを知りたい人の手助けになれれば嬉しいです.