test new inc

DiffractionGratings/grating3D.geo

@@ -205,15 +205,15 @@ pts_LAYER_L1() = PointsOf{ Physical Volume{7}; };
 pts_PMLTOP()   = PointsOf{ Physical Volume{8}; };
 pts_ROD()      = PointsOf{ Physical Volume{9}; };
 
-Characteristic Length{:}              = PML_lc/10;
-// Characteristic Length{:}              = PML_lc;
-// Characteristic Length{pts_LAYER_L1()} = L1_lc/refine_mesh_L1;
-// Characteristic Length{pts_LAYER_L2()} = L2_lc/refine_mesh_L2;
-// Characteristic Length{pts_LAYER_L6()} = L6_lc/refine_mesh_L6;
-// Characteristic Length{pts_LAYER_L5()} = L5_lc/refine_mesh_L5;
-// Characteristic Length{pts_LAYER_L4()} = L4_lc/refine_mesh_L4;
-// Characteristic Length{pts_LAYER_L3()} = L3_lc/refine_mesh_L3;
-// Characteristic Length{pts_ROD()}      = scat_lc;
+// Characteristic Length{:}              = PML_lc/5;
+Characteristic Length{:}              = PML_lc;
+Characteristic Length{pts_LAYER_L1()} = L1_lc/refine_mesh_L1;
+Characteristic Length{pts_LAYER_L2()} = L2_lc/refine_mesh_L2;
+Characteristic Length{pts_LAYER_L6()} = L6_lc/refine_mesh_L6;
+Characteristic Length{pts_LAYER_L5()} = L5_lc/refine_mesh_L5;
+Characteristic Length{pts_LAYER_L4()} = L4_lc/refine_mesh_L4;
+Characteristic Length{pts_LAYER_L3()} = L3_lc/refine_mesh_L3;
+Characteristic Length{pts_ROD()}      = scat_lc;
 If (tag_geom==3) // Split torus weird otherwise
   Mesh.Algorithm = 6;
 EndIf

DiffractionGratings/grating3D.pro

@@ -150,13 +150,13 @@ Function{
   rp[] = (CompZ[k1[]]*epsr2[]-CompZ[k2[]]*epsr1[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]);
   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[];
+  spol[]    = Vector[Sin[phi0],-Cos[phi0],0];
+  AmplEis[] =      spol[];
+  AmplErs[] = rs[]*spol[];
+  AmplEts[] = ts[]*spol[];
+  AmplHis[] = Sqrt[ep0*epsr1[]/mu0]      *spol[];
+  AmplHrs[] = Sqrt[ep0*epsr1[]/mu0]*rp[]*spol[];
+  AmplHts[] = Sqrt[ep0*epsr1[]/mu0]*tp[]*spol[];
 
   Eis[]     = AmplEis[] * Exp[I[]*k1[] *XYZ[]];
   Ers[]     = AmplErs[] * Exp[I[]*k1r[]*XYZ[]];
@@ -167,44 +167,13 @@ Function{
   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[]];
-  Hi[]     =  1/(om0*mu0*mur[]) * k1[]  /\ Ei[];
-  Hr[]     =  1/(om0*mu0*mur[]) * k1r[] /\ Er[];
-  Ht[]     =  1/(om0*mu0*mur[]) * k2[]  /\ Et[];
+
+  Ei[] = Ae*(Cos[psi0]*Eip[]-Sin[psi0]*Eis[]);
+  Er[] = Ae*(Cos[psi0]*Erp[]-Sin[psi0]*Ers[]);
+  Et[] = Ae*(Cos[psi0]*Etp[]-Sin[psi0]*Ets[]);
+  Hi[] =  1/(om0*mu0*mur[]) * k1[]  /\ Ei[];
+  Hr[] =  1/(om0*mu0*mur[]) * k1r[] /\ Er[];
+  Ht[] =  1/(om0*mu0*mur[]) * k2[]  /\ Et[];
 
   E1[Omega_super]  = Ei[]+Er[];
   E1[Omega_subs]   = Et[];