import numpy as np
import math
import cmath
import matplotlib.pyplot as plt
im = complex(0.0,1.0)
# Si
ka = 150.0
ca = 710.0
ra = 2330.0
ta = 150.0e-12
va = math.sqrt(ka/(ra*ca*ta))
#Ge
kb = 60.0
cb = 310.0
rb = 5300.0
tb = 200.0e-12
vb = math.sqrt(kb/(rb*cb*tb))
d_1 = 2.00e-9
fa_1 = 0.5
da_1 = fa_1*d_1
db_1 = d_1-da_1
d_2 = 2.35e-9
fa_2 = 0.5
da_2 = fa_2*d_2
db_2 = d_2-da_2
d_3 = 2.75e-9
fa_3 = 0.5
da_3 = fa_3*d_3
db_3 = d_3-da_3
d_4 = 3.25e-9
fa_4 = 0.5
da_4 = fa_4*d_4
db_4 = d_4-da_4
Nf = 1000
f = np.zeros(Nf+1)
n2a = np.zeros(Nf+1,dtype=complex)
na = np.zeros(Nf+1,dtype=complex)
ga = np.zeros(Nf+1,dtype=complex)
Za = np.zeros(Nf+1,dtype=complex)
n2b = np.zeros(Nf+1,dtype=complex)
nb = np.zeros(Nf+1,dtype=complex)
gb = np.zeros(Nf+1,dtype=complex)
Zb = np.zeros(Nf+1,dtype=complex)
r = np.zeros(Nf+1,dtype=complex)
t = np.zeros(Nf+1,dtype=complex)
logaritmo = np.zeros(Nf+1,dtype=complex)
fi = 0.01
ff = 800.0e9
df = (ff-fi)/Nf
for iff in range(0,Nf+1):
f[iff] = fi+df*iff
w = 2.0*math.pi*f[iff]
n2a[iff] = 1+im/(w*ta)
na[iff] = cmath.sqrt(n2a[iff])
ga[iff] = na[iff]*(w/va)
Za[iff] = (ga[iff]*ka)/(w*ta*n2a[iff])
n2b[iff] = 1+im/(w*tb)
nb[iff] = cmath.sqrt(n2b[iff])
gb[iff] = nb[iff]*(w/vb)
Zb[iff] = (gb[iff]*kb)/(w*tb*n2b[iff])
Ma = [[1,1],[Za[iff],-Za[iff]]]
Mai = np.linalg.inv(Ma)
Mb = [[1,1],[Zb[iff],-Zb[iff]]]
Mbi = np.linalg.inv(Mb)
Fa1 = [[cmath.exp(im*ga[iff]*da_1),0],[0,cmath.exp(-im*ga[iff]*da_1)]]
Fb1 = [[cmath.exp(im*gb[iff]*db_1),0],[0,cmath.exp(-im*gb[iff]*db_1)]]
Fa2 = [[cmath.exp(im*ga[iff]*da_2),0],[0,cmath.exp(-im*ga[iff]*da_2)]]
Fb2 = [[cmath.exp(im*gb[iff]*db_2),0],[0,cmath.exp(-im*gb[iff]*db_2)]]
Fa3 = [[cmath.exp(im*ga[iff]*da_3),0],[0,cmath.exp(-im*ga[iff]*da_3)]]
Fb3 = [[cmath.exp(im*gb[iff]*db_3),0],[0,cmath.exp(-im*gb[iff]*db_3)]]
Fa4 = [[cmath.exp(im*ga[iff]*da_4),0],[0,cmath.exp(-im*ga[iff]*da_4)]]
Fb4 = [[cmath.exp(im*gb[iff]*db_4),0],[0,cmath.exp(-im*gb[iff]*db_4)]]
Mi = Mb
Mt = Ma
Mti = np.linalg.inv(Mt)
Maux1_1 = np.matmul(Fa1,Mai)
Maux1_1 = np.matmul(Ma, Maux1_1)
Maux1_1 = np.matmul(Mbi,Maux1_1)
Maux1_1 = np.matmul(Fb1,Maux1_1)
Maux1_1 = np.matmul(Mb, Maux1_1)
Maux1_2 = np.matmul(Fa2,Mai)
Maux1_2 = np.matmul(Ma, Maux1_2)
Maux1_2 = np.matmul(Mbi,Maux1_2)
Maux1_2 = np.matmul(Fb2,Maux1_2)
Maux1_2 = np.matmul(Mb, Maux1_2)
Maux1_3 = np.matmul(Fa3,Mai)
Maux1_3 = np.matmul(Ma, Maux1_3)
Maux1_3 = np.matmul(Mbi,Maux1_3)
Maux1_3 = np.matmul(Fb3,Maux1_3)
Maux1_3 = np.matmul(Mb, Maux1_3)
Maux1_4 = np.matmul(Fa4,Mai)
Maux1_4 = np.matmul(Ma, Maux1_4)
Maux1_4 = np.matmul(Mbi,Maux1_4)
Maux1_4 = np.matmul(Fb4,Maux1_4)
Maux1_4 = np.matmul(Mb, Maux1_4)
# 1 hetero, 3 bicapas
Maux = np.matmul(Maux1_1,Mi)
Maux = np.matmul(Maux1_1,Maux)
Maux = np.matmul(Maux1_1,Maux)
Maux = np.matmul(Maux1_1,Maux)
Maux = np.matmul(Maux1_1,Maux)
Maux = np.matmul(Maux1_1,Maux)
Maux = np.matmul(Maux1_1,Maux)
Maux = np.matmul(Maux1_1,Maux)
# 2 hetero, 3 bicapas
Maux = np.matmul(Maux1_2,Maux)
Maux = np.matmul(Maux1_2,Maux)
Maux = np.matmul(Maux1_2,Maux)
Maux = np.matmul(Maux1_2,Maux)
Maux = np.matmul(Maux1_2,Maux)
Maux = np.matmul(Maux1_2,Maux)
Maux = np.matmul(Maux1_2,Maux)
# 3 hetero, 3 bicapas
# Maux = np.matmul(Maux1_3,Maux)
# Maux = np.matmul(Maux1_3,Maux)
# Maux = np.matmul(Maux1_3,Maux)
# Maux = np.matmul(Maux1_3,Maux)
# Maux = np.matmul(Maux1_3,Maux)
# 4 hetero, 3 bicapas
# Maux = np.matmul(Maux1_4,Maux)
# Maux = np.matmul(Maux1_4,Maux)
# Maux = np.matmul(Maux1_4,Maux)
# Maux = np.matmul(Maux1_4,Maux)
# Maux = np.matmul(Maux1_4,Maux)
# ultima capa
Maux = np.matmul(Mti, Maux)
r[iff] = -Maux[1,0]/Maux[1,1]
t[iff] = Maux[0,0]-(Maux[1,0]*Maux[0,1])/Maux[1,1]
logaritmo[iff]=10*math.log10(abs(t[iff]))