Include "scat2D_data.geo";
cel      = 1.0/Sqrt[epsilon0 * mu0];
myDir = "output/";

Group {
  // Domains
	PMLxy      = Region[1001];
	PMLx       = Region[1002];
	PMLy       = Region[1003];
	SrcPt      = Region[1004];
	ScatShape  = Region[1005];
	BackGnd    = Region[{1006,SrcPt}];
  
  Omega_pml    = Region[{PMLxy,PMLx,PMLy}];
	Omega_nopml  = Region[{ScatShape,BackGnd}];
  Omega        = Region[{Omega_nopml,Omega_pml}];
  // we do not plot exactly at the source point to avoid singularity
  If (Flag_PW==0)
    Omega_plot       = Region[{Omega,-SrcPt}];
  Else
    Omega_plot       = Region[{Omega}];
  EndIf
  Omega_plot_nopml = Region[{Omega_plot,-Omega_pml}];

}

Function{
  I[] = Complex[0.0,1.0];

  epsr_shape[] = Complex[epsr_shape_re,epsr_shape_im];
  epsr_bg[]    = Complex[epsr_bg_re,epsr_bg_im];

  mur_shape[] = Complex[mur_shape_re,mur_shape_im];
  mur_bg[]    = 1;
	
  Freq   = cel/lambda0;
  omega0 = 2.*Pi*cel/lambda0;
  k0     = 2.*Pi/lambda0;
  n_bg[] = (epsr_bg[])^(0.5);
  k_bg[] = k0*n_bg[];
	
  a_pml =    1.;
  b_pml = -siwt;
	
	sx[PMLxy]       = Complex[a_pml,b_pml];
	sx[PMLx]        = Complex[a_pml,b_pml];
	sx[PMLy]        = 1;
	sx[Omega_nopml] = 1;

	sy[PMLxy]       = Complex[a_pml,b_pml];
	sy[PMLx]        = 1;
	sy[PMLy]        = Complex[a_pml,b_pml];
	sy[Omega_nopml] = 1;

	sz = 1;
  
	PML_Tensor[] = TensorDiag[sz*sy[]/sx[],sx[]*sz/sy[],sx[]*sy[]/sz];

  epsilonr[ScatShape] = epsr_shape[] * TensorDiag[1,1,1];
  epsilonr[BackGnd]   = epsr_bg[]    * TensorDiag[1,1,1];
  epsilonr[Omega_pml] = epsr_bg[]    * PML_Tensor[];
  
  epsilonr_annex[ScatShape] = epsr_bg[] * TensorDiag[1,1,1];
  epsilonr_annex[BackGnd]   = epsr_bg[] * TensorDiag[1,1,1];
  epsilonr_annex[Omega_pml] = epsr_bg[] * PML_Tensor[];

  mur[ScatShape] = mur_shape[] * TensorDiag[1,1,1];
  mur[BackGnd]   = mur_bg[]    * TensorDiag[1,1,1];
  mur[Omega_pml] = mur_bg[]    * PML_Tensor[];

  mur_annex[ScatShape] = TensorDiag[1,1,1];
  mur_annex[BackGnd]   = TensorDiag[1,1,1];
  mur_annex[Omega_pml] = PML_Tensor[];
	
	
	// source : Plane Wave
	A                  = 1.;
  alpha[]            =  k_bg[]*Sin[theta_PW_deg*deg2rad];
  beta_bg[]          =  k_bg[]*Cos[theta_PW_deg*deg2rad];
  ui_PW[Omega_nopml] = A*Exp[ siwt*I[]*(alpha[]*X[]+beta_bg[]*Y[]) ];
  ui_PW[Omega_pml]   = 0.;
	Grad_ui_PW[]       = siwt*Vector[alpha[],beta_bg[],0.] * I[] * ui_PW[];
	
	// 2D Green function  Δe+k^2 e = -δ
  // Hankel function of the 2nd kind of order 0
	XYZ_src[] = Vector[x_src,y_src,0];
  R_src[]   = Norm[XYZ[] - XYZ_src[]]; 
  ui_Green2D[Omega_nopml] = I[]/4 * Complex[Jn[0,k_bg[]*R_src[]],-siwt*Yn[0,k_bg[]*R_src[]]];
  ui_Green2D[Omega_pml]   = 0.;
  // Gradient of the Green function (for TM case)
  Grad_ui_Green2D[] = -I[]/4 * k_bg[] 
								* Complex[Jn[1,k_bg[]*R_src[]],-siwt*Yn[1,k_bg[]*R_src[]]]
								* (XYZ[] - XYZ_src[])/R_src[];

  If(Flag_PW==0)
    ui[]      = ui_Green2D[];
    Grad_ui[] = Grad_ui_Green2D[];
  Else
    ui[]      = ui_PW[];
    Grad_ui[] = Grad_ui_PW[];
  EndIf

  // Compute the tensors xsi and chi to handle both E// and H// cases
  If (flag_Hparallel==1)
	  detepst[]       = CompXX[epsilonr[]]*CompYY[epsilonr[]]-CompXY[epsilonr[]]*Conj[CompYX[epsilonr[]]];
	  detepst_annex[] = CompXX[epsilonr_annex[]]*CompYY[epsilonr_annex[]]-CompXY[epsilonr_annex[]]*Conj[CompYX[epsilonr_annex[]]];
	  xsi[]           = Transpose[epsilonr[]]/detepst[];
	  xsi_annex[]     = Transpose[epsilonr_annex[]]/detepst_annex[];
	  chi[]           = CompZZ[mur[]];
		chi_annex[]     = CompZZ[mur_annex[]];		
  Else
    detmut[]       = CompXX[mur[]]*CompYY[mur[]]-CompXY[mur[]]*Conj[CompYX[mur[]]];
  	detmut_annex[] = CompXX[mur_annex[]]*CompYY[mur_annex[]]-CompXY[mur_annex[]]*Conj[CompYX[mur_annex[]]];
	  xsi[]          = Transpose[mur[]]/detmut[];
	  xsi_annex[]    = Transpose[mur_annex[]]/detmut_annex[];
	  chi[]          = CompZZ[epsilonr[]];
	  chi_annex[]    = CompZZ[epsilonr_annex[]];
  EndIf

	// -1. comes from the integration by part
  source_xsi[] = (xsi[]-xsi_annex[]) * -1. * Grad_ui[];
  source_chi[] = k0^2 * (chi[]-chi_annex[]) * ui[];
}

Constraint {
}

Jacobian {
  { Name JVol ;
    Case { 
      { Region All ; Jacobian Vol ; }
    }
  }
}

Integration {
  { Name Int_1 ;
    Case { 
      { Type Gauss ;
        Case { 
          { GeoElement Point    ; NumberOfPoints  1 ; }
          { GeoElement Line     ; NumberOfPoints  4 ; }
          { GeoElement Triangle ; NumberOfPoints  12 ; }
        }
      }
    }
  }
}

FunctionSpace {
  { Name Hgrad; Type Form0;
    BasisFunction {
      If (FE_lagrange_order==1)
        { Name sn;  NameOfCoef un;  Function BF_Node;    Support Region[Omega]; Entity NodesOf[All]; }
      EndIf
      If (FE_lagrange_order==2)
        { Name sn;  NameOfCoef un;  Function BF_Node;    Support Region[Omega]; Entity NodesOf[All]; }
        { Name sn2; NameOfCoef un2; Function BF_Node_2E; Support Region[Omega]; Entity EdgesOf[All]; }
      EndIf
      If (FE_lagrange_order==3)
        { Name sn;  NameOfCoef un;  Function BF_Node;    Support Region[Omega]; Entity NodesOf[All]; }
        { Name sn2; NameOfCoef un2; Function BF_Node_2E; Support Region[Omega]; Entity EdgesOf[All]; }
        { Name sn3; NameOfCoef un3; Function BF_Node_2F; Support Region[Omega]; Entity FacetsOf[All]; }
        { Name sn4; NameOfCoef un4; Function BF_Node_3E; Support Region[Omega]; Entity EdgesOf[All]; }
      EndIf
    }
    Constraint {
    }
  }
}

Formulation {
  {Name form_helmholtz_scalar; Type FemEquation;
    Quantity {
      { Name us; Type Local; NameOfSpace Hgrad;}}
    Equation {
      Galerkin { [ k0^2 * chi[] * Dof{us} , {us}]     ; In Omega; Jacobian JVol; Integration Int_1; }
      Galerkin { [ -xsi[] *  Dof{Grad us} , {Grad us}]; In Omega; Jacobian JVol; Integration Int_1; }
      Galerkin { [ source_xsi[] , {Grad us}]; In Omega; Jacobian JVol; Integration Int_1; }
      Galerkin { [ source_chi[] , {us}     ]; In Omega; Jacobian JVol; Integration Int_1; }
    }
  }
  {Name form_helmholtz_scalar_tot; Type FemEquation;
    Quantity {
      { Name u; Type Local; NameOfSpace Hgrad;}}
    Equation {
      Galerkin { [ k0^2 * chi[] * Dof{u} , {u}]    ; In Omega; Jacobian JVol; Integration Int_1; }
      Galerkin { [ -xsi[] * Dof{Grad u} , {Grad u}]; In Omega; Jacobian JVol; Integration Int_1; }
      // tiny disk source
      Galerkin { [1/(Pi*r_src^2) , {u}]; In SrcPt; Jacobian JVol; Integration Int_1;  }
    }
  }
}

Resolution {
  { Name res_helmholtz_scalar;
    System {
      { Name S; NameOfFormulation form_helmholtz_scalar; Type ComplexValue; Frequency Freq;}
    }
    Operation {
      CreateDir[Str[myDir]];
      Printf["%f",lambda0/nm];
      Generate[S] ;
      Solve[S] ;
    }
  }
  { Name res_helmholtz_scalar_tot;
    System {
      { Name S; NameOfFormulation form_helmholtz_scalar_tot; Type ComplexValue; Frequency Freq;}
    }
    Operation {
      CreateDir[Str[myDir]];
      Printf["%f",lambda0/nm];
      Generate[S] ;
      Solve[S] ;
    }
  }
}

PostProcessing {
  { Name postpro_helmholtz_scalar; NameOfFormulation form_helmholtz_scalar;
    Quantity {
      { Name xsi  ; Value { Local { [ CompXX[xsi[]]]; In Omega_plot; Jacobian JVol; } } }
      { Name chi  ; Value { Local { [ CompXX[chi[]]]; In Omega_plot; Jacobian JVol; } } }
      If (flag_Hparallel==1)
	      { Name Hz_inc  ; Value { Local { [ ui[]     ]; In Omega; Jacobian JVol; } } }
	      { Name Hz_scat ; Value { Local { [ {us}     ]; In Omega; Jacobian JVol; } } }
				{ Name Hz_tot  ; Value { Local { [ {us}+ui[]]; In Omega; Jacobian JVol; } } }
      EndIf
			If (flag_Hparallel==0)
	      { Name Ez_inc  ; Value { Local { [ ui[]     ]; In Omega; Jacobian JVol; } } }
	      { Name Ez_scat ; Value { Local { [ {us}     ]; In Omega; Jacobian JVol; } } }
				{ Name Ez_tot  ; Value { Local { [ {us}+ui[]]; In Omega; Jacobian JVol; } } }
      EndIf
	  }
  }


  { Name postpro_helmholtz_scalar_tot; NameOfFormulation form_helmholtz_scalar_tot;
    Quantity {
      { Name xsi  ; Value { Local { [ CompXX[xsi[]]]; In Omega_plot; Jacobian JVol; } } }
      { Name chi  ; Value { Local { [ CompXX[chi[]]]; In Omega_plot; Jacobian JVol; } } }
      If (flag_Hparallel==1)
        // { Name Grad_Hz_tot ; Value { Local { [ {Grad u} ]; In Omega_plot; Jacobian JVol; } } }
        // { Name Grad_ui_Green2D ; Value { Local { [ Grad_ui_Green2D[] ]; In Omega_plot; Jacobian JVol; } } }
	      { Name Hz_inc  ; Value { Local { [ ui[] ]; In Omega; Jacobian JVol; } } }
	      { Name Hz_scat ; Value { Local { [ {u}-ui[] ]; In Omega; Jacobian JVol; } } }
				{ Name Hz_tot  ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }

      EndIf
			If (flag_Hparallel==0)
        // { Name Grad_Ez_tot ; Value { Local { [ {Grad u} ]; In Omega_plot; Jacobian JVol; } } }
        // { Name Grad_ui_Green2D ; Value { Local { [ Grad_ui_Green2D[] ]; In Omega_plot; Jacobian JVol; } } }
	      { Name Ez_inc  ; Value { Local { [ ui[] ]; In Omega; Jacobian JVol; } } }
	      { Name Ez_scat ; Value { Local { [ {u}-ui[] ]; In Omega; Jacobian JVol; } } }
				{ Name Ez_tot  ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
      EndIf
	  }
  }
}

PostOperation {
  { Name postop_helmholtz_scalar; NameOfPostProcessing postpro_helmholtz_scalar ;
    Operation {
      // Print[ xsi  , OnElementsOf Omega, File StrCat[myDir, "xsi.pos"]  ];
      // Print[ chi  , OnElementsOf Omega, File StrCat[myDir, "chi.pos"]  ];
      If (flag_Hparallel==1)
				Print[ Hz_inc    , OnElementsOf Omega_plot_nopml, File StrCat[myDir, "Hz_inc.pos"]  ];
				Print[ Hz_tot    , OnElementsOf Omega_plot_nopml, File StrCat[myDir, "Hz_tot.pos"]  ];
				Print[ Hz_scat   , OnElementsOf Omega_plot, File StrCat[myDir, "Hz_scat.pos"] ];
			EndIf
			If (flag_Hparallel==0)
				Print[ Ez_inc    , OnElementsOf Omega_plot_nopml, File StrCat[myDir, "Ez_inc.pos"]  ];
				Print[ Ez_tot    , OnElementsOf Omega_plot_nopml, File StrCat[myDir, "Ez_tot.pos"]  ];
				Print[ Ez_scat   , OnElementsOf Omega_plot, File StrCat[myDir, "Ez_scat.pos"] ];
      EndIf
		}
  }

  { Name postop_helmholtz_scalar_tot; NameOfPostProcessing postpro_helmholtz_scalar_tot ;
    Operation {
      // Print[ xsi  , OnElementsOf Omega, File StrCat[myDir, "xsi.pos"]];
      // Print[ chi  , OnElementsOf Omega, File StrCat[myDir, "chi.pos"]];
      If (flag_Hparallel==1)
			  Print[ Hz_inc  , OnElementsOf Omega_plot_nopml, File StrCat[myDir, "Hz_inc.pos"]];
				Print[ Hz_tot  , OnElementsOf Omega_plot, File StrCat[myDir, "Hz_tot.pos"]];
			  Print[ Hz_scat , OnElementsOf Omega_plot_nopml, File StrCat[myDir, "Hz_scat.pos"]];
			EndIf
			If (flag_Hparallel==0)
			  Print[ Ez_inc  , OnElementsOf Omega_plot_nopml, File StrCat[myDir, "Ez_inc.pos"]];
				Print[ Ez_tot  , OnElementsOf Omega_plot, File StrCat[myDir, "Ez_tot.pos"]];
			  Print[ Ez_scat , OnElementsOf Omega_plot_nopml, File StrCat[myDir, "Ez_scat.pos"]];
      EndIf
		}
  }
}

If (Flag_PW==0)
  If (FLAG_TOT==1)
    DefineConstant[
      R_ = {"res_helmholtz_scalar_tot", Name "GetDP/1ResolutionChoices", Visible 1},
      P_ = {"postop_helmholtz_scalar_tot", Name "GetDP/2PostOperationChoices", Visible 1}
    ];
  EndIf
EndIf