##
##   Symmetric Galerkin 3D Laplace
##   Edge Adjacent LINEAR 
##   Log Term
##
with(linalg):
readlib(fortran):
##
## define arrays
##
mq   := array(1..3):  mp  := array(1..3):
Q    := array(1..3):  P   := array(1..3): R  := array(1..3):
JN   := array(1..3):  jn  := array(1..3):
N    := array(1..3): 
jc   := array(1..2);
LL   := array(1..1);
cnum := array(0..3);
CFS  := array(0..2,0..2);
rc   := array(1..3,0..1);
chks := array(0..3,1..3);
#
sq3 := sqrt(3);
##
## 3 noded triangle shape functions
##
mq[1] :=  (sq3*(1-ee) - xx)/(2*sq3);
mq[2] :=  (sq3*(1+ee) - xx)/(2*sq3);
mq[3] :=   xx/sq3;
#
mp[1] :=  (sq3*(1-e) - x)/(2*sq3);
mp[2] :=  (sq3*(1+e) - x)/(2*sq3);
mp[3] :=   x/sq3;
##
## WLOG
##
xp[1] := 0; yp[1] := 0; zp[1] := 0;
            yp[2] := 0; zp[2] := 0;
                        zp[3] := 0;
xq[2] := 0; yq[2] := 0; zq[2] := 0;
xq[1] := xp[2]; yq[1] := 0; zq[1] := 0;
##
xp[2] := 1;
#zq[3] := 0;
##
##
Q[1] := 0; Q[2] := 0; Q[3] := 0;
P[1] := 0; P[2] := 0; P[3] := 0;
for j from 1 to 3 do
 Q[1] := Q[1] + xq[j]*mq[j];
 Q[2] := Q[2] + yq[j]*mq[j];
 Q[3] := Q[3] + zq[j]*mq[j];
#
 P[1] := P[1] + xp[j]*mp[j];
 P[2] := P[2] + yp[j]*mp[j];
 P[3] := P[3] + zp[j]*mp[j];
od;
##  normals
n11 := diff(P[1],e); n12 := diff(P[2],e); n13 := diff(P[3],e);
n21 := diff(P[1],x); n22 := diff(P[2],x); n23 := diff(P[3],x);
JN1 := n12*n23-n22*n13;
JN2 := -(n11*n23-n13*n21);
JN3 := n11*n22-n12*n21;
JN[1] := JN1; JN[2] := JN2; JN[3] :=  JN3;
N[1] := 0; N[2] := 0; N[3] :=  1;
##
n11 := diff(Q[1],ee); n12 := diff(Q[2],ee); n13 := diff(Q[3],ee);
n21 := diff(Q[1],xx); n22 := diff(Q[2],xx); n23 := diff(Q[3],xx);
jn1 := n12*n23-n22*n13;
jn2 := -(n11*n23-n13*n21);
jn3 := n11*n22-n12*n21;
##
jn[1] := jn1; jn[2] := jn2; jn[3] := jn3;
ndN := JN[3]*jn[3];
##
P := add(P,N,1,eps);
R := add(Q,P,1,-1);
for j from 1 to 3 do
 R[j] := subs(ee=rho*ct-e,xx=rho*st,R[j]);
 R[j] := subs(rho=Lambda*cp,x=Lambda*sp,R[j]);
 R[j] := expand(R[j]);
 R[j] := collect(R[j],Lambda);
 for ll from 0 to 1 do
  rc[j,ll] := coeff(R[j],Lambda,ll);
  rc[j,ll] := collect(rc[j,ll],[ct,st]);
 od;
od;
##
##  r**2 = eps^2 + eps*b1*Lambda + b2*Lambda**2
##  b[1] := 2*(rc[1,0]*rc[1,1] + rc[2,0]*rc[2,1] + rc[3,0]*rc[3,1] )/eps;
##  b[2] := rc[1,1]^2 + rc[2,1]^2 + rc[3,1]^2;
##  jnr := jq1*Lambda-eps*ndN/jp
##  jNr := jp1*Lambda-jp*eps
##  jq1 := jn dot rc[*,1]   jp1 := jN dot rc[*,1]
##
jc[1] := jn[1]*rc[1,1] + jn[2]*rc[2,1] + jn[3]*rc[3,1];
jc[2] := JN[3]*rc[3,1];
jss   := expand(jc[1]*jc[2]);
##
## Log term
##
b2 := rc[1,1]^2 + rc[2,1]^2 + rc[3,1]^2;
b2 := expand(b2):
b2 := collect(b2,[cp,sp]);
cc := coeff(b2,cp,2);
cs := coeff(b2,cp,1)/sp;
ss := coeff(b2,cp,0)/sp^2;
## c2=cc c1=cs c0=ss
check := expand(b2 - (cc*cp^2+cs*cp*sp+ss*sp^2));
b2n  := c2*cp^2+c1*cp*sp+c0*sp^2;
##
##
#Lterm := cp*( - ndN/b2n^(3/2));
#Lterm := cp*( 3*jss/b2n^(5/2) );
Lterm := cp*( 3*jss/b2n^(5/2) - ndN/b2n^(3/2) );
Lterm := expand(Lterm);
#Lterm := subs(sq3^2=3,Lterm);
##
## q = cotan(p)   dq = -1/sp^2 * dp   0<psi<pi/2 => 0<q<infinity
##
Lterm := subs(
              (c2*cp^2+c1*cp*sp+c0*sp^2)^(-3/2)=
                -(c2*q^2+c1*q+c0)^(-3/2)*sp^(-1),
              (c2*cp^2+c1*cp*sp+c0*sp^2)^(-5/2)=
                -(c2*q^2+c1*q+c0)^(-5/2)*sp^(-3),
          Lterm);
Lterm := subs(cp=q*sp,Lterm);
Lterm := int(Lterm,q=0..infinity);
##
##
## second integration
##
##
## q = cotan(t)   dq = -1/st^2 * dt 
##
CC := subs(ct=q*st,cc);
CC := collect(CC,q);
CS := subs(ct=q*st,cs);
CS := collect(CS,q);
SS := subs(ct=q*st,ss);
SS := collect(SS,q);
c2 := st^2*(c22*q^2+c21*q+c20);
c1 := st*(c11*q+c10);
c0 := c00;
##
##  coefficients  ckj = CFS[k,j]
##
CFS[2,2] := expand(coeff(CC,q,2)/st^2);
CFS[2,1] := expand(coeff(CC,q,1)/st^2);
CFS[2,0] := expand(coeff(CC,q,0)/st^2);
CFS[1,1] := expand(coeff(CS,q,1)/st);
CFS[1,0] := expand(coeff(CS,q,0)/st);
CFS[0,0] := expand(coeff(SS,q,0));
#
Lterm := -st^2*Lterm;
Lterm := subs((st^2*(c22*q^2+c21*q+c20))^(-1/2)=
              st^(-1)*(c22*q^2+c21*q+c20)^(-1/2),
              (st^2*(c22*q^2+c21*q+c20))^(1/2)=
              st*(c22*q^2+c21*q+c20)^(1/2),
              (st^2*c22*q^2+st^2*c21*q+st^2*c20)^(1/2)=
              st*(c22*q^2+c21*q+c20)^(1/2),
          Lterm);
Lterm := expand(Lterm);
Lterm := subs(
          (4*st^2*c00*c22*q^2+4*st^2*c00*c21*q+4*st^2*c00*c20
           -st^2*c11^2*q^2-2*st^2*c11*q*c10-st^2*c10^2)^(-2)=
          st^(-4)*(4*c00*c22*q^2+4*c00*c21*q+4*c00*c20
           -c11^2*q^2-2*c11*q*c10-c10^2)^(-2),
          Lterm);
Lterm := subs(
              1/sqrt(c22*q^2+c21*q+c20)=1/FF,
          Lterm);
Lterm :=  collect(Lterm,FF);
LTF0  := coeff(Lterm,FF,0);
CHECK := expand(LTF0 + coeff(Lterm,FF,-1)/FF - Lterm);
LTF0  :=  int(LTF0,q=-infinity..infinity);
LTF0  := subs(
             d00=4*c00*c20-c10^2,
             d11=4*c00*c21-2*c11*c10,
             d22=4*c00*c22-c11^2,
             c00=CFS[0,0],
             c10=CFS[1,0],
             c11=CFS[1,1],
             c20=CFS[2,0],
             c21=CFS[2,1],
             c22=CFS[2,2],
            LTF0);
 LTF0 := expand(LTF0);
#
#LTF1  := ( (c22*q^2+c21*q+c20)*coeff(Lterm,FF,1) + coeff(Lterm,FF,-1) )
#           /sqrt(c22*q^2+c21*q+c20);
LTF1  := coeff(Lterm,FF,-1)/sqrt(c22*q^2+c21*q+c20);
LTF1  :=  subs(
           (4*c00*c22*q^2+4*c00*c21*q+4*c00*c20
           -c11^2*q^2-2*c11*q*c10-c10^2)^(-2)=
           (d00 + d11*q + d22*q^2)^(-2),
## d00 = 4*c00*c20-c10^2
## d11 = 4*c00*c21-2*c11*c10
## d22 = 4*c00*c22-c11^2
           LTF1);
LTF1  :=  normal(expand(LTF1));
den  := denom(LTF1);
num := numer(LTF1); num := collect(expand(num),q);
for jj from 0 to 3 do
 cnum[jj] := coeff(num,q,jj);
od;
chk3 := expand(num - (cnum[0]+cnum[1]*q+cnum[2]*q^2+cnum[3]*q^3));
##
##
  LTsum := 0;
for jj from 0 to 3 do
 LTj := q^jj/den;
 LTj := int(LTj,q=-infinity..infinity):
 LTj := convert(LTj,ln,arctanh):
 LTj := expand(LTj):
##
##  d11^2-4*d22*d00 < 0
##
LTj := subs(
           sqrt(d11^2-4*d22*d00)=DD,
           1/sqrt(d11^2-4*d22*d00)=1/DD,
           1/(d11^2-4*d22*d00)^(3/2)=1/DD^3,
           1/(-d11^2+4*d22*d00)=-1/DD^2,
         LTj):
LTj := subs(
(2*c22*d11^2/d22^2-4*c22*d00/d22-2*c22*d11*DD/d22^2-
 2*c21*d11/d22+2*c21*DD/d22+4*c20)^(-1/2)=
2^(-1/2)*(c22*d11^2/d22^2-2*c22*d00/d22-c22*d11*DD/d22^2-
 c21*d11/d22+c21*DD/d22+2*c20)^(-1/2),
(2*c22*d11^2/d22^2-4*c22*d00/d22+2*c22*d11*DD/d22^2-
 2*c21*d11/d22-2*c21*DD/d22+4*c20)^(-1/2)=
2^(-1/2)*(c22*d11^2/d22^2-2*c22*d00/d22+c22*d11*DD/d22^2-
 c21*d11/d22-c21*DD/d22+2*c20)^(-1/2),
         LTj):
LTj := subs(
(c22*d11^2/d22^2-2*c22*d00/d22-c22*d11*DD/d22^2-
 c21*d11/d22+c21*DD/d22+2*c20)^(-1/2)=FFa,
(c22*d11^2/d22^2-2*c22*d00/d22+c22*d11*DD/d22^2-
 c21*d11/d22-c21*DD/d22+2*c20)^(-1/2)=FFb,
         LTj):
LTj := subs(
(c22*d11^2-2*c22*d00*d22+c22*d11*DD-
 c21*d11*d22-c21*DD*d22+2*c20*d22^2)^(-1)=FFb^2/d22^2,
(-c22*d11^2+2*c22*d00*d22-c22*d11*DD+
 c21*d11*d22+c21*DD*d22-2*c20*d22^2)^(-1)=-FFb^2/d22^2,
(-c22*d11^2+2*c22*d00*d22+c22*d11*DD+
 c21*d11*d22-c21*DD*d22-2*c20*d22^2)^(-1)=-FFa^2/d22^2,
(-c22*d11^2-c21*DD*d22+2*c22*d00*d22+c22*d11*DD+
 c21*d11*d22-2*c20*d22^2)^(-1)=-FFa^2/d22^2,
         LTj):
LTj := subs(
ln(2+sqrt(2)*FFb*c21/sqrt(c22)-sqrt(2)*FFb*sqrt(c22)*d11/d22-sqrt(2)*FFb*sqrt(c22)*DD/d22)=
  ln(2) + LQ1,
ln(2-sqrt(2)*FFb*c21/sqrt(c22)+sqrt(2)*FFb*sqrt(c22)*d11/d22+sqrt(2)*FFb*sqrt(c22)*DD/d22)=
  ln(2) + LQ2,
ln(1+sqrt(2)*FFa*c21/sqrt(c22)/2-sqrt(2)*FFa*sqrt(c22)*d11/d22/2+
   sqrt(2)*FFa*sqrt(c22)*DD/d22/2)=LQ1c,
ln(1-sqrt(2)*FFa*c21/sqrt(c22)/2+sqrt(2)*FFa*sqrt(c22)*d11/d22/2-
   sqrt(2)*FFa*sqrt(c22)*DD/d22/2)= LQ2c,
         LTj):
LTj := expand(LTj):
LTj := subs(
(c22*d11^2-2*c22*d00*d22+c22*d11*DD-
 c21*d11*d22-c21*DD*d22+2*c20*d22^2)^(-1)=FFb^2/d22^2,
(-c22*d11^2+2*c22*d00*d22+c22*d11*DD+
 c21*d11*d22-c21*DD*d22-2*c20*d22^2)^(-1)=-FFa^2/d22^2,
(-c22*d11^2-c21*DD*d22+2*c22*d00*d22+c22*d11*DD+
 c21*d11*d22-2*c20*d22^2)^(-1)=-FFa^2/d22^2,
(-c22*d11^2+2*c22*d00*d22-c22*d11*DD+
 c21*d11*d22+c21*DD*d22-2*c20*d22^2)^(-1)=-FFb^2/d22^2,
         LTj):
LTj := subs(
ln(2+sqrt(2)*FFb*c21/sqrt(c22)-sqrt(2)*FFb*sqrt(c22)*d11/d22-sqrt(2)*FFb*sqrt(c22)*DD/d22)=
  ln(2) + LQ1,
ln(2-sqrt(2)*FFb*c21/sqrt(c22)+sqrt(2)*FFb*sqrt(c22)*d11/d22+sqrt(2)*FFb*sqrt(c22)*DD/d22)=
  ln(2) + LQ2,
ln(1+sqrt(2)*FFa*c21/sqrt(c22)/2-sqrt(2)*FFa*sqrt(c22)*d11/d22/2+
   sqrt(2)*FFa*sqrt(c22)*DD/d22/2)=LQ1c,
ln(1-sqrt(2)*FFa*c21/sqrt(c22)/2+sqrt(2)*FFa*sqrt(c22)*d11/d22/2-
   sqrt(2)*FFa*sqrt(c22)*DD/d22/2)= LQ2c,
         LTj):
LTj := collect(LTj,[LQ1,LQ2,LQ1c,LQ2c]):
cf1  := coeff(LTj,LQ1,1):
cf1c := coeff(LTj,LQ1c,1):
cf2  := coeff(LTj,LQ2,1):
cf2c := coeff(LTj,LQ2c,1):
chks[jj,1] := expand(cf1 - subs(FFa=ffb,FFb=FFa,ffb=FFb,DD=-DD,cf1c)):
chks[jj,2] := expand(cf2 - subs(FFa=ffb,FFb=FFa,ffb=FFb,DD=-DD,cf2c)):
chks[jj,3] := expand(cf1+cf2):
cf0 := subs(LQ1=0,LQ2=0,LQ1c=0,LQ2c=0,LTj):
LTsum := LTsum + cnum[jj]*(cf0 + cf1*LQ12):
## LQ12 = LQ1 - LQ2
od:
LTsum := subs(c22=1/4,4^(1/2)=2,LTsum);
##########
LTsum := expand(LTsum):
LTsum := collect(LTsum,LQ12):
ltc0  := normal(coeff(LTsum,LQ12,0));
ltc0  := normal(expand(ltc0));
ltc0 := ltc0*d22^2*DD^2/FFa^2/FFb^2;
##  d22^2*DD^2/(FFa^2*FFb^2)
FFp := expand(
(c22*d11^2/d22^2-2*c22*d00/d22-c22*d11*DD/d22^2-
 c21*d11/d22+c21*DD/d22+2*c20)*
(c22*d11^2/d22^2-2*c22*d00/d22+c22*d11*DD/d22^2-
 c21*d11/d22-c21*DD/d22+2*c20)):
FFp := d22^2*DD^2*FFp:
FFp := expand(subs(DD^2=d11^2-4*d22*d00,FFp)):
FFp := subs(
                d00 = 4*c00*c20-c10^2,
                d11 = 4*c00*c21-2*c11*c10,
                d22 = 4*c00*c22-c11^2,
                 c00=CFS[0,0],
                 c10=CFS[1,0],
                 c11=CFS[1,1],
                 c20=CFS[2,0],
                 c21=CFS[2,1],
                 c22=CFS[2,2],
        FFp):
FFp := factor(expand(FFp));
##
##
ltc0  := subs(
               d00 = 4*c00*c20-c10^2,
               d11 = 4*c00*c21-2*c11*c10,
               d22 = 4*c00*c22-c11^2,
                 c00=CFS[0,0],
                 c10=CFS[1,0],
                 c11=CFS[1,1],
                 c20=CFS[2,0],
                 c21=CFS[2,1],
                 c22=CFS[2,2],
       ltc0):
ltc0  := factor(expand(ltc0/FFp));
##########
ltcL  := normal(coeff(LTsum,LQ12,1));
ltcL  := (ltcL/FFa^2/FFb^3/2^(1/2)/yp[3])*(144*d22^3*DD^5);
ltcL  := subs(
               d00 = 4*c00*c20-c10^2,
               d11 = 4*c00*c21-2*c11*c10,
               d22 = 4*c00*c22-c11^2,
       ltcL):
ltcL  := factor(expand(ltcL));
ltcL  := subs(c22=1/4,4^(1/2)=2,ltcL);
ltcL  := collect(expand(ltcL),DD);
z3    := coeff(ltcL,DD,3);
z0    := coeff(ltcL,DD,0):
check5 := expand(ltcL - z0 -z3*DD^3);
z3    := subs(
                 c00=CFS[0,0],
                 c10=CFS[1,0],
                 c11=CFS[1,1],
                 c20=CFS[2,0],
                 c21=CFS[2,1],
                 c22=CFS[2,2],
       z3):
z3    := factor(expand(z3));
## ==> z3 = 0
z0    := factor(expand(z0));
z0    := z0/c00^3/32;
z0    := z0/(c10-2*c21*c11);
z0    := z0/(c10^2-4*c21*c11*c10+4*c20*c11^2)^2;
z0    := subs(
                 c00=CFS[0,0],
                 c10=CFS[1,0],
                 c11=CFS[1,1],
                 c20=CFS[2,0],
                 c21=CFS[2,1],
                 c22=CFS[2,2],
       z0):
z0    := factor(expand(z0));
## ==> z3 = 0