ファラデー効果の偏光回転をPythonで再現

import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy import pi,sin,cos,tan,arcsin,linspace
from matplotlib.pyplot import plot,show,xlabel,ylabel,title,legend,grid,axis,tight_layout

n1=1 #媒質1の屈折率
n2=1.5 #媒質2の屈折率

t1_deg=linspace(0,90,91)#入射角t1(deg)
t1=t1_deg/180*pi#入射角t1
t2=arcsin((n1/n2)*sin(t1))#屈折角t2

tp=2*n1*cos(t1)/(n2*cos(t1)+n1*cos(t2))#P偏光透過係数
rp=(n2*cos(t1)-n1*cos(t2))/(n2*cos(t1)+n1*cos(t2))#P偏光反射係数
ts=2*n1*cos(t1)/(n1*cos(t1)+n2*cos(t2))#S偏光透過係数
rs=(n1*cos(t1)-n2*cos(t2))/(n1*cos(t1)+n2*cos(t2))#S偏光反射係数

Rp=rp**2#P偏光反射率
Rs=rs**2#S偏光反射率

plt.rcParams['xtick.direction'] = 'in' 
plt.rcParams['ytick.direction'] = 'in' 
plt.figure(figsize=(6,6))
plot(t1Deg,Rp,label=r"$r_{12}^{\rm{p}}$",linewidth=3.0,color='red')
plot(t1Deg,Rs,label=r"$r_{12}^{\rm{s}}$",linewidth=3.0,color='black')
xlabel(r"$\theta_1$ (deg.)",fontsize=20)
ylabel("Reflectivity",fontsize=20)
grid(True)
axis([0.0,90,0,1])
plt.tick_params(labelsize=20)
tight_layout()
show()

#直線偏光
from sympy import *
from scipy import pi,sin,cos,tan,arcsin,linspace
import numpy as np

plt.rcParams['xtick.direction'] = 'in' 
plt.rcParams['ytick.direction'] = 'in' 

#変数設定
list_incident=[(w, 2*pi*5.6e+14), (c, 3e+8),(z,1e-3),(Ex,1),(Ey,0),(n,1)]#パラメータを変える場合はここを変更
x,y,z,Ex,Ey,Ez,w,t = symbols("x y z Ex Ey Ez w t")
r=Matrix([x, y, z]) 
E0 = Matrix([Ex, Ey, Ez]) 
K= Matrix([0, 0, w*n/c])#光が進む方向のベクトル(波数ベクトル)#w:角振動数、n:物質中の屈折率、c:真空中での光速

#直線偏光を設定
E_incident=E0*exp(-1j*w*t)*exp(1j*(K.T*r)[0,0])
E_incident=E_incident.subs(list_incident)
#print(E_incident)

#直線偏光を描画
t3_deg=linspace(0,180,181)#角度を0-180度の範囲で振る
#t3=t3_deg/180*pi

plt.figure(figsize=(7,7)) 
for t_ in t3:
    plt.scatter(re(E_incident[0].subs(t,t_/180*pi)),re(E_incident[1].subs(t,t_/180*pi)),c="red")

plt.xlim(-1,1)
plt.ylim(-0.2,0.2)
xlabel(r"$Ex$",fontsize=20)
ylabel(r"$Ey$",fontsize=20)
plt.show()

#z軸方向に磁化をもつ磁性体中の誘電率テンソル
e_xx,e_xy,e_yy,e_zz = symbols("e_xx e_xy e_yy e_zz ",real=True)
epsilon=Matrix([[e_xx, e_xy, 0],
                  [-e_xy, e_xx, 0],
                  [0, 0, e_zz]])
#入射偏光を定義 波数ベクトルK=N_hat*w/c
c,w = symbols("c w",real=True)
N = symbols("N")
N_hat=Matrix([0, 0, N])
x,y,z,t = symbols("x y z t",real=True)
Ex,Ey,Ez= symbols("Ex Ey Ez")
r=Matrix([x, y, z]) 
E0 = Matrix([Ex, Ey, Ez]) 
E=E0*exp(-1j*w*(t-(N_hat.T*r)[0,0]/c))

#マクスウェル方程式から導かれる、平面波が満たすべき式(固有方程式)
eq=(N**2)*E-(E.T*N_hat)[0,0]*N_hat-epsilon*E
#display(eq)

#固有値を導出
Coef=Matrix([[eq[0].coeff(Ex, 1),eq[0].coeff(Ey, 1)],[eq[1].coeff(Ex, 1),eq[1].coeff(Ey, 1)]])#Ex,Eyの項の係数(2×2行列)抽出
Coef=simplify((Coef)/(exp(-1j*w*(t-(N_hat.T*r)[0,0]/c))))#全ての項に共通の部分を除去。simplifyで式整理
display(Coef)
Det=Coef.det()
#print(Det)
s=solve(Det,N**2)#固有値の2乗の値を導出
display(s)

Ey1=solve((s[0]-e_xx)*Ex-e_xy*Ey,Ey)#固有値を代入して固有ベクトルのy成分を導出(Ey1,Ey2)
Ey2=solve((s[1]-e_xx)*Ex-e_xy*Ey,Ey)
print('Ey1:',solve((s[0]-e_xx)*Ex-e_xy*Ey,Ey))
print('Ey2:',solve((s[1]-e_xx)*Ex-e_xy*Ey,Ey))
#print(Ey_1)
N1,N2 = symbols("N1 N2")

E1=Matrix([Ex, Ey1[0], Ez])*exp(-1j*w*(t-N1/c*z))#固有ベクトルその１
E2=Matrix([Ex, Ey2[0], Ez])*exp(-1j*w*(t-N2/c*z))#固有ベクトルその２

#固有ベクトルの1つに着目→円偏光であることを確認
print("固有ベクトルの1つ:",E1)
#print(E2[1])
print("x成分:",simplify(re(E1[0].coeff(Ex, 1))))#変数の状態で、固有ベクトルの1つのx,y成分の係数を導出
print("y成分:",simplify(re(E1[1].coeff(Ex, 1))))



list_epsilon=[(e_xx, 16),(e_xy,0.3j)]
n1=sqrt(s[0]).subs(list_epsilon)#sが2つ解を持ち、固有値は±√|s|の計4つ。±√sの±それぞれでの2つの固有円偏光の結果は変わらないため、プラス側の2つを以下では考える。)
n2=sqrt(s[1]).subs(list_epsilon)
#パラメータの数値代入
list_all=[(w, 2*pi*5.6e+14), (c, 3e+8),(z,1e-3),(N1,n1),(N2,n2),(Ex,1)]

t3_deg=linspace(0,180,181)
t3=t3_deg/180*pi#入射角t1 
#1つの固有ベクトルを描画→円偏光であることを確認
f0=re(E1[0].coeff(Ex, 1))#固有ベクトル1のx成分
f1=re(E1[1].coeff(Ex, 1))#固有ベクトル1のy成分
g0=f0.subs(list_all)#固有ベクトル1のx成分にパラメータ代入
g1=f1.subs(list_all)#固有ベクトル1のy成分にパラメータ代入
#print('f1:',f1)
#print('g1:',g1)
print(simplify(re(E1[0].coeff(Ex, 1))))
plt.figure(figsize=(7,7)) 
for t_ in t3:
    plt.scatter(g0.subs(t,t_),g1.subs(t,t_),c="red")
#plt.xlim(-1,1)
#plt.ylim(-0.2,0.2)
xlabel(r"$Ex$",fontsize=20)
ylabel(r"$Ey$",fontsize=20)
plt.show()



#透過光(２つの固有ベクトルの和)を導出
Transmitted_light=simplify(E1+E2)

print("Transmitted light↓")
display(Transmitted_light)
#print(E1[0])
Transmitted_x=simplify(Transmitted_light[0].subs(list_all))#透過光のx成分にパラメータ代入
Transmitted_y=simplify(Transmitted_light[1].subs(list_all))#透過光のy成分にパラメータ代入
print('Transmitted_x:',Transmitted_x)
print('Transmitted_y:',Transmitted_y)
#透過光描画
plt.figure(figsize=(7,7)) 
for t_ in t3:
    plt.scatter(re(Transmitted_x.subs(t,t_)),re(Transmitted_y.subs(t,t_)),c="red")
plt.xlim(-1,1)
plt.ylim(-0.2,0.2)
xlabel(r"$Ex$",fontsize=20)
ylabel(r"$Ey$",fontsize=20)
plt.show()