Include "PlaneWaveBiperiodic3D_data.geo";
Group {
  // SubDomains
  PML_bot = Region[1];
  Subs    = Region[2];
  Super   = Region[3];
  PML_top = Region[4];

  // Boundaries
  SurfBlochXm    = Region[101];
  SurfBlochXp    = Region[102];
  SurfBlochYm    = Region[201];
  SurfBlochYp    = Region[202];
  SurfIntTop     = Region[301];
  SurfIntBot     = Region[302];

  // Grouping
  SurfYOZ  = Region[{SurfBlochXm,SurfBlochXp}];
  SurfXOZ  = Region[{SurfBlochYm,SurfBlochYp}];
  Omega = Region[{PML_bot,Subs,Super,PML_top}];
  Domain = Region[{Subs,Super}];

  // Optional Flag to exlude facets from Bloch linking
  // Since Webb elements are hierarcical, we just remove
  // facet basis functions from periodic faces, so Bloch ties
  // only relies on edge basis functions
  If (FLAG_EXCLUDE_FACETS_FROM_LINK==1)
    SurfExcludeFacets = Region[{SurfXOZ,SurfYOZ}];
  Else
    SurfExcludeFacets = Region[{}];
  EndIf
}


Function{
  I[]    = Complex[0,1];
  k0     = 2*Pi / lambda0;
  om0    = cel * k0;

  // TODO
  epsr1[] = Complex[epsr_super_re, 0];
  epsr2[] = Complex[epsr_subs_re, epsr_subs_im];
  n1[]  = Sqrt[epsr1[]];
  n2[]  = Sqrt[epsr2[]];
  period_x = d;
  period_y = d;

  // handling possibly skewed lattices
  xsi = 0;

  // plane wave amplitude dans power entering the cell
  Ae   = 1;
  Pinc = 0.5*Ae^2*Sqrt[epsr_super_re*ep0/mu0] * Cos[theta0];

  // simple PML strecth along Z
  sx   = 1.;
	sy   = 1.;
	sz[] = Complex[1,1];

  // Permittivities
  epsr[Super]  = epsr1[] * TensorDiag[1,1,1];
  epsr[Subs]   = epsr2[] * TensorDiag[1,1,1];
  epsr[PML_top] = epsr_super_re * TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
  epsr[PML_bot] = epsr_subs_re * TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
  
  // Permeabilities
  mur[Super]  = TensorDiag[1,1,1];
  mur[Subs]   = TensorDiag[1,1,1];
	mur[PML_top] = TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
	mur[PML_bot] = TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
  mur[SurfIntTop]  = TensorDiag[1,1,1];

  // Analytical solution of the simple diopter
  // = Fresnel coefficients for a 3D plane wave
  // we project the 3D problem on the s/p (TE/TM) basis 
  // to use scalar well know fresnel coefficients

  // incident, transmitted & reflected wave vectors
  k1x[] = -k0*n1[]*Sin[theta0]*Cos[phi0];
  k1y[] = -k0*n1[]*Sin[theta0]*Sin[phi0];
  k1z[] = -k0*n1[]*Cos[theta0] ;
  k2x[] =  k1x[];
  k2y[] =  k1y[];
  k2z[] = -Sqrt[k0^2*epsr2[]-k1x[]^2-k1y[]^2];
  k1[]  = Vector[k1x[],k1y[], k1z[]];
  k2[]  = Vector[k2x[],k2y[], k2z[]];
  k1r[] = Vector[k1x[],k1y[],-k1z[]];

  // spol is a unit vector perpendicular to the plane of incidence (Oz,k1)
  spol[] = Vector[Sin[phi0], -Cos[phi0],0];
  // ppol is in the plane of incidence, orthogonal to k_{1,2} and spol
  ppol_r[] = (k1r[]/Norm[k1r[]]) /\ spol[];
  ppol_t[] = (k2[] /Norm[k2[]] ) /\ spol[];

  // r/t Fresnel coefficients in 2D
  // p "parallel"  polar : E is parallel      to the plane of incidence (Oz,k1)
  // s "senkrecht" polar : E is perpendicular to the plane of incidence (Oz,k1)
  rs[] = (k1z[]-k2z[])/(k1z[]+k2z[]);
  ts[] =     2.*k1z[] /(k1z[]+k2z[]);
  rp[] = (k1z[]*epsr2[]-k2z[]*epsr1[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
  tp[] =            (2.*k1z[]*epsr2[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
  
  // that's where scalar Fresnel coefficients are used
  // E in TE (E is orthogonal to the plane of incidence)
  AmplEis[] =      spol[];
  AmplErs[] = rs[]*spol[];
  AmplEts[] = ts[]*spol[];
  // H in TM (H is orthogonal to the plane of incidence)
  AmplHis[] = Sqrt[ep0*epsr1[]/mu0]     *spol[];
  AmplHrs[] = Sqrt[ep0*epsr1[]/mu0]*rp[]*spol[];
  AmplHts[] = Sqrt[ep0*epsr1[]/mu0]*tp[]*spol[];

  // just multiply by the appropriate phasor to get the fields
  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[]];
  
  // We need to deduce the E-field corresponding the (TM) H-field
  Eip[]     = -1/(om0*ep0*epsr1[]) * k1[]  /\ His[];
  Erp[]     = -1/(om0*ep0*epsr1[]) * k1r[] /\ Hrs[];
  Etp[]     = -1/(om0*ep0*epsr2[]) * k2[]  /\ Hts[];

  // Now that we have our basis, the final field with ad hoc polarization
  // is a linear combination of Es and Ep
  If (FLAG_CIRCPOL==1)
    // NOTE : Sign[k1z[]] convention for left/right handed polar
    Ei[] = Ae/Sqrt[2]*(Sign[k1z[]]*Eip[]+FLAG_LRCIRC*I[]*Eis[]);
    Er[] = Ae/Sqrt[2]*(Sign[k1z[]]*Erp[]+FLAG_LRCIRC*I[]*Ers[]);
    Et[] = Ae/Sqrt[2]*(Sign[k1z[]]*Etp[]+FLAG_LRCIRC*I[]*Ets[]);
  Else
    Ei[] = Ae*(Cos[psi0]*Eip[]-Sin[psi0]*Eis[]);
    Er[] = Ae*(Cos[psi0]*Erp[]-Sin[psi0]*Ers[]);
    Et[] = Ae*(Cos[psi0]*Etp[]-Sin[psi0]*Ets[]);
  EndIf

  // corresponding H-field
  Hi[] =  1/(om0*mu0*1) * k1[]  /\ Ei[];
  Hr[] =  1/(om0*mu0*1) * k1r[] /\ Er[];
  Ht[] =  1/(om0*mu0*1) * k2[]  /\ Et[];

  // define the analytical E-field Eana by parts
  Eana[#{SurfIntTop,Super}] = Ei[]+Er[];
  Eana[#{SurfIntBot,Subs}]  = Et[];
  Eana[#{PML_bot,PML_top}]  = 0;

  // define the corresponding analytical H-field Hana by parts
  Hana[#{SurfIntTop,Super}] = Hi[]+Hr[];
  Hana[#{SurfIntBot,Subs}]  = Ht[];
  Hana[#{PML_bot,PML_top}]  = 0;
  // Now we have the analytical solution

  
  // virtual antenna at (PML_top/Super) interface 
  source_surf_tot[] = -2*I[]*om0*mu0*Vector[0,0,1] /\ Hi[];

  // Bloch phase shifts
  dephX[] = Exp[I[]*k1x[]*d];
  dephY[] = Exp[I[]*k1y[]*d];

  // To compute the energy balance, we just expand the solution as
  // a Rayleigh sum and take order 0 (specular).
  // We could use the Poynting vector, but we would need to rely on the curl
  // of the numerical solution... Here, we only rely on
  // tangential components of the E-field (readily available).
  // we will use this for both analytical and numerical solutions.
  Nb_ordre = 2*Nmax+1;
  For i In {0:Nb_ordre-1}
    For j In {0:Nb_ordre-1}
      alpha~{i}~{j}[] = k1x[] + 2*Pi/period_x*(i-Nmax);
      beta~{i}~{j}[]  = k1y[] + 2*Pi/period_y*(j-Nmax)/Cos[xsi] - 2*Pi/period_x*(i-Nmax)*Tan[xsi];
      expialphaxy~{i}~{j}[] = Exp[-I[]*(alpha~{i}~{j}[]*X[]+beta~{i}~{j}[]*Y[])];
    EndFor
  EndFor
  For i In {0:Nb_ordre-1}
    For j In {0:Nb_ordre-1}
      gammar~{i}~{j}[] = Sqrt[k0^2*epsr1[] - alpha~{i}~{j}[]^2 - beta~{i}~{j}[]^2];
      gammat~{i}~{j}[] = Sqrt[k0^2*epsr2[] - alpha~{i}~{j}[]^2 - beta~{i}~{j}[]^2];
    EndFor
  EndFor

}

// Apply Bloch link
Constraint {
  { Name BlochX;
    Case {
    { Region SurfBlochXp; Type LinkCplx ; RegionRef SurfBlochXm;
      Coefficient dephX[]; Function Vector[$X-d,$Y,$Z] ; }
    }
  }
  { Name BlochY;
    Case {
      { Region SurfBlochYp; Type LinkCplx ; RegionRef SurfBlochYm;
      Coefficient dephY[]; Function Vector[$X,$Y-d,$Z] ; }
    }
  }
}

Jacobian {
  { Name JVol ; Case {{ Region All ; Jacobian Vol ; }}}
  { Name JSur ; Case {{ Region All ; Jacobian Sur ; }}}
  { Name JLin ; Case {{ Region All ; Jacobian Lin ; }}}
}

// we should adapt the number of Gauss points to flag_interp...
Integration {
  { Name I1 ;
    Case {
      { Type Gauss ;
        Case {
          { GeoElement Point       ; NumberOfPoints   1 ; }
          { GeoElement Line        ; NumberOfPoints   4 ; }
          { GeoElement Triangle    ; NumberOfPoints  12 ; }
          { GeoElement Triangle2   ; NumberOfPoints  12 ; }
          { GeoElement Tetrahedron ; NumberOfPoints  16 ; }
          { GeoElement Tetrahedron2; NumberOfPoints  16 ; }
        }
      }
    }
  }
}

// 4 discrete spaces "H^1(\Omega,Curl,bloch)"
// 1- Webb order 1 (same as Nedelec) = one DOF per edge
// 2- two DOFs per edge (this is an intermediate order)
// 3- two DOFs per edge, one per facet (this is an intermediate order)
// 4- Webb order 2 = two DOFs per edge, two DOFs per facets
FunctionSpace {
  { Name Hcurl; Type Form1;
    BasisFunction {
      If (flag_interp == 1)
        { Name sn;  NameOfCoef un;  Function BF_Edge;      Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
      EndIf
      If (flag_interp == 2)
        { Name sn;  NameOfCoef un;  Function BF_Edge;      Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn2; NameOfCoef un2; Function BF_Edge_2E;   Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
      EndIf
      If (flag_interp == 3)
        { Name sn;  NameOfCoef un;  Function BF_Edge;      Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn2; NameOfCoef un2; Function BF_Edge_2E;   Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }
      EndIf
      If (flag_interp == 4)
        { Name sn;  NameOfCoef un;  Function BF_Edge;      Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn2; NameOfCoef un2; Function BF_Edge_2E;   Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }
        { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }          
      EndIf
      If (flag_interp == 5)
        { Name sn;  NameOfCoef un;  Function BF_Edge;      Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn2; NameOfCoef un2; Function BF_Edge_2E;   Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }
        { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }          
        { Name sn5; NameOfCoef un5; Function BF_Edge_3F_c; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }          
      EndIf
      If (flag_interp == 6)
        { Name sn;  NameOfCoef un;  Function BF_Edge;      Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn2; NameOfCoef un2; Function BF_Edge_2E;   Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
        { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }
        { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }
        { Name sn5; NameOfCoef un5; Function BF_Edge_3F_c; Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity FacetsOf[All, Not SurfExcludeFacets]; }
        { Name sn6; NameOfCoef un6; Function BF_Edge_4E;   Support Region[{Omega,SurfIntTop,SurfIntBot}]; Entity EdgesOf[All]; }
      EndIf
    }
    Constraint {
      If (flag_interp == 1)
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochY; }
      EndIf
      If (flag_interp == 2)
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
      EndIf
      If (flag_interp == 3)
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochY; }
      EndIf
      If (flag_interp == 4)
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochY; }
        { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochY; }
      EndIf
      If (flag_interp == 5)
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochY; }
        { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochY; }
        { NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint BlochY; }
      EndIf
      If (flag_interp == 6)
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochY; }
        { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochY; }
        { NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint BlochX; }
        { NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint BlochY; }
        { NameOfCoef un6; EntityType EdgesOf ; NameOfConstraint BlochX; }
        { NameOfCoef un6; EntityType EdgesOf ; NameOfConstraint BlochY; }
      EndIf
    }
  }
}

Formulation {
  { Name form_helmholtz_vector; Type FemEquation;
    Quantity {
      { Name u  ; Type Local; NameOfSpace Hcurl; }
    }
    Equation{
      Galerkin { [1/mur[]  * Dof{Curl u} , {Curl u}]; In Omega      ; Jacobian JVol; Integration I1; }
      Galerkin { [-k0^2*epsr[] * Dof{u}  ,      {u}]; In Omega      ; Jacobian JVol; Integration I1; }
      Galerkin { [-source_surf_tot[]     ,      {u}]; In SurfIntTop ; Jacobian JVol; Integration I1; }
    }
  }
}

Resolution {
  { Name res_helmholtz_vector;
    System {
      { Name S; NameOfFormulation form_helmholtz_vector; Type ComplexValue;}
    }
    Operation {
      Generate[S];
      Solve[S];
    }
  }
}

PostProcessing {
  { Name postpro_helmholtz_vector; NameOfFormulation form_helmholtz_vector; NameOfSystem S;
    Quantity {
      { Name epsr_xx ; Value { Local { [ CompXX[epsr[]] ]; In Omega; Jacobian JVol; } } }
      { Name E ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
      { Name H ; Value { Local { [ (-I[]/(mur[]*mu0*om0)*{Curl u}) ]; In Omega; Jacobian JVol; } } }      
      { Name Eana ; Value { Local { [ Eana[] ]; In Omega; Jacobian JVol; } } }
      { Name Hana ; Value { Local { [ Hana[] ]; In Omega; Jacobian JVol; } } }
      { Name diff_E ; Value { Local { [ {u}-Eana[] ]; In Domain; Jacobian JVol; } } }
      { Name diff_H ; Value { Local { [ (-I[]/(mur[]*mu0*om0)*{Curl u})-Hana[] ]; In Domain; Jacobian JVol; } } }
      { Name diff_E_L2norm ; Value { Integral { [Norm[{u}-Eana[]] / (d*d*(h_subs)) ] ; In Subs ; Integration I1 ; Jacobian JVol ; } } }
      { Name diff_H_L2norm ; Value { Integral { [Norm[(-I[]/(mur[]*mu0*om0)*{Curl u})-Hana[]] / (d*d*(h_subs)) ] ; In Subs ; Integration I1 ; Jacobian JVol ; } } }
      { Name Ei ; Value { Local { [ Ei[] ]; In Omega; Jacobian JVol; } } }
      For i In {0:Nb_ordre-1}
        For j In {0:Nb_ordre-1}
          { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u}     ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
          { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}     ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
          { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
          { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
          { Name eff_t1~{i}~{j}   ; Value { Term{ Type Global; [
                1/(Ae^2*gammat~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * (
                    (gammat~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_t~{i}~{j}]
                   +(gammat~{i}~{j}[]^2+ beta~{j}~{j}[]^2)*SquNorm[$int_y_t~{i}~{j}]
                   +2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_t~{i}~{j}*Conj[$int_y_t~{i}~{j}]] 
                  ) ] ; In SurfIntBot ; } } }
          { Name eff_r1~{i}~{j}   ; Value { Term{ Type Global; [
                1/(Ae^2*gammar~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * (
                   (gammar~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_r~{i}~{j}]
                  +(gammar~{i}~{j}[]^2+ beta~{i}~{j}[]^2)*SquNorm[$int_y_r~{i}~{j}]
                  +2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_r~{i}~{j}*Conj[$int_y_r~{i}~{j}]]
                ) ] ; In SurfIntTop ; } } }
          { Name int_x_t_ana~{i}~{j} ; Value { Integral { [ CompX[Eana[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
          { Name int_y_t_ana~{i}~{j} ; Value { Integral { [ CompY[Eana[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
          { Name int_x_r_ana~{i}~{j} ; Value { Integral { [ CompX[Eana[]-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
          { Name int_y_r_ana~{i}~{j} ; Value { Integral { [ CompY[Eana[]-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
          { Name eff_t1_ana~{i}~{j}  ; Value { Term{ Type Global; [
                1/(Ae^2*gammat~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * (
                   (gammat~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_t_ana~{i}~{j}]
                  +(gammat~{i}~{j}[]^2+ beta~{j}~{j}[]^2)*SquNorm[$int_y_t_ana~{i}~{j}]
                  +2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_t_ana~{i}~{j}*Conj[$int_y_t_ana~{i}~{j}]] 
                ) ] ; In SurfIntBot ; } } }
          { Name eff_r1_ana~{i}~{j}   ; Value { Term{ Type Global; [
                1/(Ae^2*gammar~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * (
                   (gammar~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_r_ana~{i}~{j}]
                  +(gammar~{i}~{j}[]^2+ beta~{i}~{j}[]^2)*SquNorm[$int_y_r_ana~{i}~{j}]
                  +2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_r_ana~{i}~{j}*Conj[$int_y_r_ana~{i}~{j}]]
                ) ] ; In SurfIntTop ; } } }
        EndFor
      EndFor
    }
  }
}

PostOperation {
  { Name postop_helmholtz_vector; NameOfPostProcessing postpro_helmholtz_vector ;
    Operation {
      // Print [ E , OnElementsOf Omega, File "E.pos"];
      // Print [ Eana , OnElementsOf Omega, File "Eana.pos"];
      // Print [ H , OnElementsOf Omega, File "H.pos"];
      // Print [ Hana , OnElementsOf Omega, File "Hana.pos"];
      // Print [ diff_E , OnElementsOf Domain, File "diff_E.pos"];
      // Print [ diff_H , OnElementsOf Domain, File "diff_H.pos"];
      Print[ diff_E_L2norm[Subs], OnGlobal, File "diff_E_L2.out", Format Table ];
      Print[ diff_H_L2norm[Subs], OnGlobal, File "diff_H_L2.out", Format Table ];
      For i In {0:Nb_ordre-1}
        For j In {0:Nb_ordre-1}
          Print[ int_x_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_x_t~{i}~{j}, Format Table];
          Print[ int_y_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_y_t~{i}~{j}, Format Table];
          Print[ int_x_r~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_x_r~{i}~{j}, Format Table];
          Print[ int_y_r~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_y_r~{i}~{j}, Format Table];
          Print[ int_x_t_ana~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_x_t_ana~{i}~{j}, Format Table];
          Print[ int_y_t_ana~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_y_t_ana~{i}~{j}, Format Table];
          Print[ int_x_r_ana~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_x_r_ana~{i}~{j}, Format Table];
          Print[ int_y_r_ana~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_y_r_ana~{i}~{j}, Format Table];
        EndFor
      EndFor

      For i In {0:Nb_ordre-1}
        For j In {0:Nb_ordre-1}
          Print[ eff_t1~{i}~{j}[SurfIntBot]    , OnRegion SurfIntBot, File "T.out"    , Format Table, SendToServer "5Results/T numerical"  ];
          Print[ eff_r1~{i}~{j}[SurfIntTop]    , OnRegion SurfIntTop, File "R.out"    , Format Table, SendToServer "5Results/R numerical"  ];
          Print[ eff_t1_ana~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File "T_ana.out", Format Table, SendToServer "5Results/T analytical" ];
          Print[ eff_r1_ana~{i}~{j}[SurfIntTop], OnRegion SurfIntTop, File "R_ana.out", Format Table, SendToServer "5Results/R analytical" ];
        EndFor
      EndFor
    }
  }
}

DefineConstant[
  R_ = {"res_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
  C_ = {"-solve -pos -petsc_prealloc 200 ", Name "GetDP/9ComputeCommand", Visible 1},
  P_ = {"postop_helmholtz_vector", Name "GetDP/2PostOperationChoices", Visible 1}
];
