Skip to content

Commit e767723

Browse files
author
Jon Schlipf
committed
implement + verify
1 parent b77e180 commit e767723

File tree

2 files changed

+171
-0
lines changed

2 files changed

+171
-0
lines changed

‎src/RigorousCoupledWaveAnalysis.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,8 @@ include("SRCWA/SRCWA.jl")
1919
include("ETM/ETM.jl")
2020
include("BasicMaterials/BasicMaterials.jl")
2121

22+
include("Taylor/Taylor.jl")
23+
export taylor_reftra
2224

2325
using .Common
2426
using .SRCWA

‎src/Taylor/Taylor.jl‎

Lines changed: 169 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,169 @@
1+
using LinearAlgebra
2+
using LinearAlgebra,CUDA
3+
using ..Common
4+
5+
export squarescale_reftra
6+
7+
function taylorx(dnx,dny,Kx,Ky,λ,l::PatternedLayer)
8+
k0=2π/real(λ)
9+
#get the base permittivity
10+
εxx=get_permittivity(l.materials[1],λ,1)*I
11+
#add the permittivity for all inclusions
12+
if minimum([typeof(m)<:Common.Isotropic for m in l.materials])
13+
#all isotropic
14+
εxx=get_permittivity(l.materials[1],λ)*I
15+
for ct=1:length(l.geometries)
16+
rec=reciprocal(l.geometries[ct],dnx,dny)
17+
εxx+=rec*(get_permittivity(l.materials[ct+1],λ)-get_permittivity(l.materials[ct],λ))
18+
end
19+
εzz=εyy=εxx
20+
εxy=εyx=0I
21+
else
22+
#anisotropic
23+
εxx=get_permittivity(l.materials[1],λ,1)*I
24+
εxy=get_permittivity(l.materials[1],λ,2)*I
25+
εyx=get_permittivity(l.materials[1],λ,3)*I
26+
εyy=get_permittivity(l.materials[1],λ,4)*I
27+
εzz=get_permittivity(l.materials[1],λ,5)*I
28+
for ct=1:length(l.geometries)
29+
rec=reciprocal(l.geometries[ct],dnx,dny)
30+
εxx+=rec*(get_permittivity(l.materials[ct+1],λ,1)-get_permittivity(l.materials[ct],λ,1))
31+
εxy+=rec*(get_permittivity(l.materials[ct+1],λ,2)-get_permittivity(l.materials[ct],λ,2))
32+
εyx+=rec*(get_permittivity(l.materials[ct+1],λ,3)-get_permittivity(l.materials[ct],λ,3))
33+
εyy+=rec*(get_permittivity(l.materials[ct+1],λ,4)-get_permittivity(l.materials[ct],λ,4))
34+
εzz+=rec*(get_permittivity(l.materials[ct+1],λ,5)-get_permittivity(l.materials[ct],λ,5))
35+
end
36+
end
37+
a01=0
38+
a11=-0.10036558103014462001
39+
a21=-0.00802924648241156960
40+
a31=-0.00089213849804572995
41+
42+
b01=0
43+
b11=0.39784974949964507614
44+
b21=1.36783778460411719922
45+
b31=0.49828962252538267755
46+
b41=-0.00063789819459472330
47+
48+
b02=-10.9676396052962062593
49+
b12=1.68015813878906197182
50+
b22=0.05717798464788655127
51+
b32=-0.00698210122488052084
52+
b42=0.00003349750170860705
53+
54+
b03=-0.09043168323908105619
55+
b13=-0.06764045190713819075
56+
b23=0.06759613017704596460
57+
b33=0.02955525704293155274
58+
b43=-0.00001391802575160607
59+
60+
b04=0
61+
b14=0
62+
b24=-0.09233646193671185927
63+
b34=-.01693649390020817171
64+
b44=-0.00001400867981820361
65+
η=inv(εzz)
66+
P=[Kx*η*Ky I-Kx*η*Kx;Ky*η*Ky-I -Ky*η*Kx]
67+
Q=[Kx*Ky+εyx εyy-Kx*Kx;Ky*Ky-εxx -εxy-Ky*Kx]
68+
A0=[0I P;Q 0I]*k0*l.thickness
69+
nrm=maximum(sum(abs.(A0),dims=1))
70+
m=Int(ceil(log2(nrm)))
71+
A=A0*2.0^-m
72+
PQ=P*Q
73+
PQP=PQ*P
74+
QP=Q*P
75+
QPQ=QP*Q
76+
A2=[PQ 0I;0I QP]*(k0*l.thickness*2.0^-m)^2
77+
A3=[0I PQP;QP*Q 0I]*(k0*l.thickness*2.0^-m)^3
78+
A6=[PQP*QPQ 0I;0I QPQ*PQP]*(k0*l.thickness*2.0^-m)^6
79+
#A2=A*A
80+
#A3=A2*A
81+
#A6=A3*A3
82+
B1=a01*I+a11*A+a21*A2+a31*A3
83+
B2=b01*I+b11*A+b21*A2+b31*A3+b41*A6
84+
B3=b02*I+b12*A+b22*A2+b32*A3+b42*A6
85+
B4=b03*I+b13*A+b23*A2+b33*A3+b43*A6
86+
B5=b04*I+b14*A+b24*A2+b34*A3+b44*A6
87+
A9=B1*B5+B4
88+
X=B2+(B3+A9)*A9
89+
return X^(2^m)
90+
end
91+
function squarescalex(dnx,dny,Kx,Ky,λ,l::PatternedLayer)
92+
k0=2π/real(λ)
93+
#get the base permittivity
94+
εxx=get_permittivity(l.materials[1],λ,1)*I
95+
#add the permittivity for all inclusions
96+
if minimum([typeof(m)<:Common.Isotropic for m in l.materials])
97+
#all isotropic
98+
εxx=get_permittivity(l.materials[1],λ)*I
99+
for ct=1:length(l.geometries)
100+
rec=reciprocal(l.geometries[ct],dnx,dny)
101+
εxx+=rec*(get_permittivity(l.materials[ct+1],λ)-get_permittivity(l.materials[ct],λ))
102+
end
103+
εzz=εyy=εxx
104+
εxy=εyx=0I
105+
else
106+
#anisotropic
107+
εxx=get_permittivity(l.materials[1],λ,1)*I
108+
εxy=get_permittivity(l.materials[1],λ,2)*I
109+
εyx=get_permittivity(l.materials[1],λ,3)*I
110+
εyy=get_permittivity(l.materials[1],λ,4)*I
111+
εzz=get_permittivity(l.materials[1],λ,5)*I
112+
for ct=1:length(l.geometries)
113+
rec=reciprocal(l.geometries[ct],dnx,dny)
114+
εxx+=rec*(get_permittivity(l.materials[ct+1],λ,1)-get_permittivity(l.materials[ct],λ,1))
115+
εxy+=rec*(get_permittivity(l.materials[ct+1],λ,2)-get_permittivity(l.materials[ct],λ,2))
116+
εyx+=rec*(get_permittivity(l.materials[ct+1],λ,3)-get_permittivity(l.materials[ct],λ,3))
117+
εyy+=rec*(get_permittivity(l.materials[ct+1],λ,4)-get_permittivity(l.materials[ct],λ,4))
118+
εzz+=rec*(get_permittivity(l.materials[ct+1],λ,5)-get_permittivity(l.materials[ct],λ,5))
119+
end
120+
end
121+
η=inv(εzz)
122+
P=[Kx*η*Ky I-Kx*η*Kx;Ky*η*Ky-I -Ky*η*Kx]
123+
Q=[Kx*Ky+εyx εyy-Kx*Kx;Ky*Ky-εxx -εxy-Ky*Kx]
124+
A0=[0I P;Q 0I]*k0*l.thickness
125+
nrm=maximum(sum(abs.(A0),dims=1))
126+
m=Int(ceil(log2(nrm)))
127+
m=0
128+
A=A0*2.0^-m
129+
X=exp(A)
130+
return X^(2^m)
131+
end
132+
function taylor_reftra(ψin,m::RCWAModel,grd::RCWAGrid,λ)
133+
X=[taylorx(grd.dnx,grd.dny,grd.Kx,grd.Ky,λ,l) for l in m.layers]
134+
Xp=I
135+
for X in X
136+
Xp*=X
137+
end
138+
ref=halfspace(grd.Kx,grd.Ky,m.εsup,λ) #superstrate and substrate
139+
tra=halfspace(grd.Kx,grd.Ky,m.εsub,λ)
140+
Y=Xp*[I;-tra.V]
141+
S=[Y [-I;-ref.V]]\[I;-ref.V]*ψin
142+
to,ro=slicehalf(S)
143+
144+
kzin=grd.k0[3]#*real(sqrt(get_permittivity(m.εsup,λ)))
145+
R=a2p(0ro,ro,ref.V,I,kzin) #compute amplitudes to powers
146+
T=-a2p(to,0to,tra.V,I,kzin)
147+
148+
return R,T
149+
150+
end
151+
function squarescale_reftra(ψin,m::RCWAModel,grd::RCWAGrid,λ)
152+
X=[squarescalex(grd.dnx,grd.dny,grd.Kx,grd.Ky,λ,l) for l in m.layers]
153+
Xp=I
154+
for X in X
155+
Xp*=X
156+
end
157+
ref=halfspace(grd.Kx,grd.Ky,m.εsup,λ) #superstrate and substrate
158+
tra=halfspace(grd.Kx,grd.Ky,m.εsub,λ)
159+
Y=Xp*[I;-tra.V]
160+
S=[Y [-I;-ref.V]]\[I;-ref.V]*ψin
161+
to,ro=slicehalf(S)
162+
163+
kzin=grd.k0[3]#*real(sqrt(get_permittivity(m.εsup,λ)))
164+
R=a2p(0ro,ro,ref.V,I,kzin) #compute amplitudes to powers
165+
T=-a2p(to,0to,tra.V,I,kzin)
166+
167+
return R,T
168+
169+
end

0 commit comments

Comments
 (0)