function [I_DE] = DIRECT_WS_ST_RWG(r1,r2,r3,Np_1D)
%% Main body of the DIRECT EVALUATION method for the
% evaluation of the coincident 4-D weakly singular integrals over planar
% triangular elements.
% Licensing: This code is distributed under the GNU LGPL license.
% Modified: 19 October 2011
% Author: Athanasios Polimeridis
% References
% 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
% r1,r2,r3 = point vectors of the triangular element's vertices
% Outer triangle P:(rp1,rp2,rp3)=(r1,r2,r3)
% Inner triangle Q:(rq1,rq2,rq3)=(r1,r2,r3)
% Np_1D = order of the Gauss-Legendre quadrature rule
% OUTPUT DATA
% I_DE(1) = I_f1_f1
% I_DE(2) = I_f1_f2
% I_DE(3) = I_f1_f3
% I_DE(4) = I_f2_f1
% I_DE(5) = I_f2_f2
% I_DE(6) = I_f2_f3
% I_DE(7) = I_f3_f1
% I_DE(8) = I_f3_f2
% I_DE(9) = I_f3_f3
%
%%
%%%%%%%%%%%%%%%%%% Triangle Area %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ap = (1/2)*norm(cross(r2-r1,r3-r1));
%%%%%%%%%%%%%%%%%% Jacobian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Jp = Ap/sqrt(3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Asing0 = norm(r1-r2)^2/4;
Asing1 = dot(r1-r2,r1+r2-2*r3)/(2*sqrt(3));
Asing2 = norm(r1+r2-2*r3)^2/12;
%
Asing3 = (1/4) *Asing0 -(sqrt(3)/4) *Asing1 +(3/4) *Asing2;
Asing4 = (sqrt(3)/2) *Asing0 -(1/2) *Asing1 -(sqrt(3)/2) *Asing2;
Asing5 = (3/4) *Asing0 +(sqrt(3)/4) *Asing1 +(1/4) *Asing2;
%
Asing6 = (1/4) *Asing0 +(sqrt(3)/4) *Asing1 +(3/4) *Asing2;
Asing7 = -(sqrt(3)/2) *Asing0 -(1/2) *Asing1 +(sqrt(3)/2) *Asing2;
Asing8 = (3/4) *Asing0 -(sqrt(3)/4) *Asing1 +(1/4) *Asing2;
%
Asing = [Asing0 Asing3 Asing6
Asing1 Asing4 Asing7
Asing2 Asing5 Asing8];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1D Gauss-Legendre Quadrature Rule %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[w,z] = Gauss_1D(Np_1D);
Isub = zeros(3,3,3);
%---------------------------------------------
for m = 1:3
%%%%%%%%% Define the coefs of the appropriate subtriangle
acc = Asing(1,m);
acs = Asing(2,m);
ass = Asing(3,m);
%%%%%%%%%%%%%%%%%%%%% Gauss quadrature %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Int1_1_1=0;Int2_1_1=0;Int3_1_1=0;
for kk = 1:Np_1D
%%%%%%%%%%%%%%% Int1, 0 =< PSI <= pi/3 %%%%%%%%%%%%%%%%%%%%%%%%%%%
PSI_1a = 0;
PSI_1b = pi/3;
PSI_1k = ((PSI_1b-PSI_1a)*z(kk)+(PSI_1b+PSI_1a))/2;
a_PSI_1k = sqrt(acc*cos(PSI_1k)^2-acs*cos(PSI_1k)*sin(PSI_1k)+ass*sin(PSI_1k)^2);
% Int1_1_1
PHI_a_1_1 = Phi_1_1(1,PSI_1k,a_PSI_1k);
PHI_e_1_1 = Phi_1_1(5,PSI_1k,a_PSI_1k);
PHI_f_1_1 = Phi_1_1(6,PSI_1k,a_PSI_1k);
F1_1_1 = PHI_a_1_1+PHI_e_1_1+PHI_f_1_1;
Int1_1_1 = Int1_1_1+w(kk)*F1_1_1;
%%%%%%%%%%%%%%% Int2, pi/3 =< PSI <= 2pi/3 %%%%%%%%%%%%%%%%%%%%%%%
PSI_2a = pi/3;
PSI_2b = 2*pi/3;
PSI_2k = ((PSI_2b-PSI_2a)*z(kk)+(PSI_2b+PSI_2a))/2;
a_PSI_2k = sqrt(acc*cos(PSI_2k)^2-acs*cos(PSI_2k)*sin(PSI_2k)+ass*sin(PSI_2k)^2);
% Int2_1_1
PHI_b_1_1 = Phi_1_1(2,PSI_2k,a_PSI_2k);
PHI_g_1_1 = Phi_1_1(7,PSI_2k,a_PSI_2k);
F2_1_1 = PHI_b_1_1+PHI_g_1_1;
Int2_1_1 = Int2_1_1+w(kk)*F2_1_1;
%%%%%%%%%%%%%%% Int3, 2pi/3 =< PSI <= pi %%%%%%%%%%%%%%%%%%%%%%%%%
PSI_3a = 2*pi/3;
PSI_3b = pi;
PSI_3k = ((PSI_3b-PSI_3a)*z(kk)+(PSI_3b+PSI_3a))/2;
a_PSI_3k = sqrt(acc*cos(PSI_3k)^2-acs*cos(PSI_3k)*sin(PSI_3k)+ass*sin(PSI_3k)^2);
% Int3_1_1
PHI_c_1_1 = Phi_1_1(3,PSI_3k,a_PSI_3k);
PHI_d_1_1 = Phi_1_1(4,PSI_3k,a_PSI_3k);
PHI_h_1_1 = Phi_1_1(8,PSI_3k,a_PSI_3k);
F3_1_1 = PHI_c_1_1+PHI_d_1_1+PHI_h_1_1;
Int3_1_1 = Int3_1_1+w(kk)*F3_1_1;
end % for kk=1:Npp
% Isub_1_1
Int1_1_1 = ((PSI_1b-PSI_1a)/2)*Int1_1_1;
Int2_1_1 = ((PSI_2b-PSI_2a)/2)*Int2_1_1;
Int3_1_1 = ((PSI_3b-PSI_3a)/2)*Int3_1_1;
Isub(1,1,m) = Int1_1_1+Int2_1_1+Int3_1_1;
end % for m = 1:3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I_DE = Jp^2 * (Isub(1,1,1)+Isub(1,1,2)+Isub(1,1,3));
function [PHI] = Phi_1_1(argument,psi,Apsi)
%% Phi_1_1 function
% Licensing: This code is distributed under the GNU LGPL license.
% Modified: 20 September 2011
% Author: Athanasios Polimeridis
% References
% 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
% argument = 1-8 -> a-h
% psi = Variable Psi
% Apsi = alpha(Psi)
% OUTPUT DATA
% Phi_1_1
%
global ko
j = sqrt(-1);
A = j*ko*Apsi;
B = cos(psi);
C = sin(psi)/sqrt(3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Int_1_1 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = sin(psi)/((j*ko)^2*Apsi^3);
switch argument
case 1
%%%%%%%%%%%%%%%%%%%%%% PHI_a %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSI_a1 = A/(2*B);
PSI_a2 = -1;
PSI_a3 = (B/A)*(1-exp(-A/B));
PSI_a = PSI_a1+PSI_a2+PSI_a3;
PHI = D*PSI_a;
case 2
%%%%%%%%%%%%%%%%%%%%%% PHI_b %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSI_b1 = A/(2*C);
PSI_b2 = -1;
PSI_b3 = (C/A)*(1-exp(-A/C));
PSI_b = PSI_b1+PSI_b2+PSI_b3;
PHI = D*PSI_b;
case 3
%%%%%%%%%%%%%%%%%%%%%% PHI_c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta = tan(pi-psi)/sqrt(3);
gamma = (1-beta)/(1+beta);
PSI_c1 = (A/C)*(1/2-(gamma-gamma^2/2));
PSI_c2 = -(1-gamma);
PSI_c3 = (C/A)*(1-exp(-A*(1-gamma)/C));
PSI_c = PSI_c1+PSI_c2+PSI_c3;
PHI = D*PSI_c;
case 4
%%%%%%%%%%%%%%%%%%%%%% PHI_d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta = tan(pi-psi)/sqrt(3);
gamma = (1-beta)/(1+beta);
PSI_d1 = -A*(gamma+gamma^2/2)/B;
PSI_d2 = -gamma;
PSI_d3 = (B/A)*(exp(A*(1+gamma)/B)-exp(A/B));
PSI_d = PSI_d1+PSI_d2+PSI_d3;
PHI = D*PSI_d;
case 5
%%%%%%%%%%%%%%%%%%%%%% PHI_e %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsilon = tan(psi)/sqrt(3);
delta = -(1-epsilon)/(1+epsilon);
PSI_e1 = A*(delta^2/2-delta)/B;
PSI_e2 = delta;
PSI_e3 = (B/A)*(exp(-A/B)-exp(-A*(1-delta)/B));
PSI_e = PSI_e1+PSI_e2+PSI_e3;
PHI = D*PSI_e;
case 6
%%%%%%%%%%%%%%%%%%%%%% PHI_f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsilon = tan(psi)/sqrt(3);
delta = -(1-epsilon)/(1+epsilon);
PSI_f1 = (A/C)*(delta+delta^2/2+1/2);
PSI_f2 = -(1+delta);
PSI_f3 = (C/A)*(1-exp(-A*(1+delta)/C));
PSI_f = PSI_f1+PSI_f2+PSI_f3;
PHI = D*PSI_f;
case 7
%%%%%%%%%%%%%%%%%%%%%% PHI_g %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSI_g1 = A/(2*C);
PSI_g2 = -1;
PSI_g3 = (C/A)*(1-exp(-A/C));
PSI_g = PSI_g1+PSI_g2+PSI_g3;
PHI = D*PSI_g;
case 8
%%%%%%%%%%%%%%%%%%%%%% PHI_h %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSI_h1 = -A/(2*B);
PSI_h2 = -1;
PSI_h3 = (B/A)*(exp(A/B)-1);
PSI_h = PSI_h1+PSI_h2+PSI_h3;
PHI = D*PSI_h;
end %switch argument