%% Example for 4D singular integrals arising in IE formulations

%  Licensing: This code is distributed under the GNU LGPL license. 

%  Modified:  8 February 2013

%  Author:    Athanasios Polimeridis

% References

% A. G. Polimeridis, ?Direct evaluation method in computational electromagnetics (DEMCEM) package,? 2011.
% [Online]. Available: http://lema.epfl.ch/Members/thanos/Software.html

% A. G. Polimeridis and T. V. Yioultsis, “On the direct evaluation of weakly singular
% integrals in Galerkin mixed potential integral equation formulations,” IEEE Trans.
% Antennas Propag., vol. 56, no. 9, pp. 3011-3019, Sep. 2008.

% A. G. Polimeridis and J. R. Mosig, “Complete semi-analytical treatment of weakly
% singular integrals on planar triangles via the direct evaluation method,” Int. J.
% Numerical Methods Eng., vol. 83, pp. 1625-1650, 2010.

% A. G. Polimeridis, J. M. Tamayo, J. M. Rius and J. R. Mosig, “Fast and accurate
% computation of hyper-singular integrals in Galerkin surface integral equation
% formulations via the direct evaluation method,” IEEE Trans.
% Antennas Propag., vol. 59, no. 6, pp. 2329-2340, Jun. 2011.

% A. G. Polimeridis and J. R. Mosig, “On the direct evaluation of surface integral
% equation impedance matrix elements involving point singularities,” IEEE Antennas
% Wireless Propag. Lett., vol. 10, pp. 599-602, 2011.

% INPUT DATA
% f      = working frequency
% dx, dy = dimensions of rectangle
% Np_1D  = order of the Gauss-Legendre quadrature rule

% OUTPUT DATA
%             I_ST  = Singular Galerkin inner product for coincident
%             elements, where the kernel is just the free-space Green's
%             function (or Helmholtz fundamental solution), i.e.,
%             G = exp(-j*k*r) / (4*pi*r)
%
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               INPUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = 300e6; % working frequency

dx = 1e-6;        % dx,dy = edge lengths for the rectangle 
dy = dx;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          EM basic parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu = 4*pi*1e-7;
co = 299792458;
eo = 1/co^2/mu;
%
omega   = 2 * pi * f;
lambda  = co/f;
%
k0 = 2*pi/lambda % wavenumber
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          Compute Self-Term
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[I_ST] = DEMCEM_ST(k0,dx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      PRINT RESULTS                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('I_ST: (%4.16e, %4.16e) \n', real(I_ST),imag(I_ST))