From 894ace7b3e1c4acde267a4ae96d95e9411bb92f4 Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Thu, 21 Nov 2019 10:28:19 +0100
Subject: [PATCH] new inc - better
---
DiffractionGratings/grating3D.pro | 72 +++++++++++++++++++++++--------
1 file changed, 54 insertions(+), 18 deletions(-)
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index f82f27f..0851ba9 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -125,9 +125,6 @@ Function{
epsr[PMLbot] = Re[epsr2[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
mur[PMLbot] = TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
- // epsr[PMLbot] = epsr2[];
- // mur[PMLbot] = TensorDiag[1,1,1];
-
// epsr[PMLtop] = TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
// mur[PMLtop] = TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
// epsr[PMLbot] = TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
@@ -147,25 +144,64 @@ Function{
-Sqrt[k0^2*epsr2[]-CompX[k1[]]^2-CompY[k1[]]^2]
];
k1r[] = Vector[CompX[k1[]], CompY[k1[]], -CompZ[k1[]]];
- k1s[] = (Abs[theta0]>1e-12) ? k1[]/\zhat[] / Norm[ k1[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
- k2s[] = (Abs[theta0]>1e-12) ? Re[k2[]]/\zhat[] / Norm[ Re[k2[]]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
- k1rs[] = (Abs[theta0]>1e-12) ? k1r[]/\zhat[] / Norm[ k1r[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
- k1p[] = (Abs[theta0]>1e-12) ? k1[]/\k1s[] / Norm[ k1[]/\k1s[]] : Vector[ Cos[phi0], Sin[phi0], 0.] ;
- k2p[] = (Abs[theta0]>1e-12) ? Re[k2[]]/\k2s[] / Norm[ Re[k2[]]/\k2s[]] : Vector[ Cos[phi0], Sin[phi0], 0.] ;
- k1rp[] = (Abs[theta0]>1e-12) ? k1r[]/\k1rs[] / Norm[ k1r[]/\k1rs[]] : Vector[-Cos[phi0],-Sin[phi0], 0.] ;
rs[] = (CompZ[k1[]]-CompZ[k2[]])/(CompZ[k1[]]+CompZ[k2[]]);
- ts[] = 2.*CompZ[k1[]] /(CompZ[k1[]]+CompZ[k2[]]);
+ ts[] = 2.*CompZ[k1[]] /(CompZ[k1[]]+CompZ[k2[]]);
rp[] = (CompZ[k1[]]*epsr2[]-CompZ[k2[]]*epsr1[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]);
- tp[] = (2.*CompZ[k1[]]*epsr2[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]) * Sqrt[epsr1[]/epsr2[]];
-
- AmplEi[] = Cos[psi0]* k1p[] + Sin[psi0]* k1s[];
- AmplEr[] = rp[]*Cos[psi0]*k1rp[] + rs[]*Sin[psi0]*k1rs[];
- AmplEt[] = tp[]*Cos[psi0]* k2p[] + ts[]*Sin[psi0]* k2s[];
+ tp[] = (2.*CompZ[k1[]]*epsr2[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]);
+
+ Vecs[] = Vector[Sin[phi0],-Cos[phi0],0];
+ AmplEis[] = Vecs[];
+ AmplErs[] = rs[]*Vecs[];
+ AmplEts[] = ts[]*Vecs[];
+ AmplHis[] = Sqrt[ep0*epsr1[]/mu0] *Vecs[];
+ AmplHrs[] = Sqrt[ep0*epsr1[]/mu0]*rp[]*Vecs[];
+ AmplHts[] = Sqrt[ep0*epsr1[]/mu0]*tp[]*Vecs[];
+
+ Eis[] = AmplEis[] * Exp[I[]*k1[] *XYZ[]];
+ Ers[] = AmplErs[] * Exp[I[]*k1r[]*XYZ[]];
+ Ets[] = AmplEts[] * Exp[I[]*k2[] *XYZ[]];
+ His[] = AmplHis[] * Exp[I[]*k1[] *XYZ[]];
+ Hrs[] = AmplHrs[] * Exp[I[]*k1r[]*XYZ[]];
+ Hts[] = AmplHts[] * Exp[I[]*k2[] *XYZ[]];
+ Eip[] = -1/(om0*ep0*epsr1[]) * k1[] /\ His[];
+ Erp[] = -1/(om0*ep0*epsr1[]) * k1r[] /\ Hrs[];
+ Etp[] = -1/(om0*ep0*epsr2[]) * k2[] /\ Hts[];
+ Hip[] = 1/(om0*mu0*mur[]) * k1[] /\ Ei[];
+ Hrp[] = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
+ Htp[] = 1/(om0*mu0*mur[]) * k2[] /\ Et[];
+
+ E1sa[Omega_super] = Eis[]+Ers[];
+ E1sa[Omega_subs] = Ets[];
+ E1pa[Omega_super] = Eip[]+Erp[];
+ E1pa[Omega_subs] = Etp[];
+ H1sa[Omega_super] = His[]+Hrs[];
+ H1sa[Omega_subs] = Hts[];
+ H1pa[Omega_super] = Hip[]+Hrp[];
+ H1pa[Omega_subs] = Htp[];
+
+ E1[] = Cos[psi0]*E1p[]-Sin[psi0]*E1s[];
+ H1[] = Cos[psi0]*H1p[]-Sin[psi0]*H1s[];
+
+ // k1s[] = (Abs[theta0]>1e-12) ? k1[]/\zhat[] / Norm[ k1[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+ // k2s[] = (Abs[theta0]>1e-12) ? Re[k2[]]/\zhat[] / Norm[ Re[k2[]]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+ // k1rs[] = (Abs[theta0]>1e-12) ? k1r[]/\zhat[] / Norm[ k1r[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+ // k1p[] = (Abs[theta0]>1e-12) ? k1[]/\k1s[] / Norm[ k1[]/\k1s[]] : Vector[ Cos[phi0], Sin[phi0], 0.] ;
+ // k2p[] = (Abs[theta0]>1e-12) ? Re[k2[]]/\k2s[] / Norm[ Re[k2[]]/\k2s[]] : Vector[ Cos[phi0], Sin[phi0], 0.] ;
+ // k1rp[] = (Abs[theta0]>1e-12) ? k1r[]/\k1rs[] / Norm[ k1r[]/\k1rs[]] : Vector[-Cos[phi0],-Sin[phi0], 0.] ;
+
+ // rs[] = (CompZ[k1[]]-CompZ[k2[]])/(CompZ[k1[]]+CompZ[k2[]]);
+ // ts[] = 2.*CompZ[k1[]] /(CompZ[k1[]]+CompZ[k2[]]);
+ // rp[] = (CompZ[k1[]]*epsr2[]-CompZ[k2[]]*epsr1[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]);
+ // tp[] = (2.*CompZ[k1[]]*epsr2[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]) * Sqrt[epsr1[]/epsr2[]];
+
+ // AmplEi[] = Cos[psi0]* k1p[] + Sin[psi0]* k1s[];
+ // AmplEr[] = rp[]*Cos[psi0]*k1rp[] + rs[]*Sin[psi0]*k1rs[];
+ // AmplEt[] = tp[]*Cos[psi0]* k2p[] + ts[]*Sin[psi0]* k2s[];
- Ei[] = Ae*AmplEi[] * Exp[I[]*k1[]*XYZ[]];
- Et[] = Ae*AmplEt[] * Exp[I[]*k2[]*XYZ[]];
- Er[] = Ae*AmplEr[] * Exp[I[]*k1r[]*XYZ[]];
+ // Ei[] = Ae*AmplEi[] * Exp[I[]*k1[]*XYZ[]];
+ // Et[] = Ae*AmplEt[] * Exp[I[]*k2[]*XYZ[]];
+ // Er[] = Ae*AmplEr[] * Exp[I[]*k1r[]*XYZ[]];
Hi[] = 1/(om0*mu0*mur[]) * k1[] /\ Ei[];
Hr[] = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
Ht[] = 1/(om0*mu0*mur[]) * k2[] /\ Et[];
--
GitLab