##
##   Symmetric Galerkin 3D Laplace
##   Coincident LINEAR 
##   hypersingular
##
##
with(linalg):
readlib(fortran):
##
## define arrays
##
phi   := array(0..2,0..1);
Lterm := array(0..2,0..1):
##
##
## first integration
##
##
rh  := (eps^2 + a^2*rho^2)^(1/2);
JNR := -eps*jp;
##
## phi integral jnJN=jp^2  (JNR)*(jnR) = eps^2*jp^2
##
hyp  :=  jp^2/rh^3 - 3*JNR^2/rh^5;
for j from 0 to 1 do
## powers for Q shape function
 zz := int(rho^(j+1)*hyp,rho=0..QR);
 zz := subs(sqrt(eps^2)=eps,1/sqrt(eps^2)=1/eps,zz);
  if j=1 then
   zz := subs( ln(eps)=loge,eps=0,
               (a^2*QR^2)^(1/2)=a*QR,
               (a^2*QR^2)^(-1/2)=1/(a*QR),
               (a^2*QR^2)^(3/2)=a^3*QR^3,
               (a^2*QR^2)^(-3/2)=1/(a^3*QR^3),
          zz);
## this loge disappears after integrating theta (st,ct coeffs from shapes)
   zz := subs( loge=0,zz);
  fi;
 zz := subs(sqrt(a^2)=a,1/sqrt(a^2)=1/a,zz);
 zz := factor(expand(zz));
##
## change of variables: theta --> t
## polar coords: t-e=Lambda*cp  x=Lambda*sp
##
##  QR = sqrt(x^2+(t-e)^2) = Lambda
##  ct = (t-e)/sqrt(x^2+(t-e)^2) = cp
##  st = -x/sqrt(x^2+(t-e)^2) = -sp
 zz  := subs(QR=Lambda,zz);
## d(theta)/dt = x/(x^2+(t-e)^2) = sp/Lambda    
## 1/Lambda cancels with Lambda d(Lambda)
 zz  :=  zz*sp;
 zz  := expand(zz);
 zz  := normal(zz);
##  new polar coords
##
## second integration
##  powers of rho for shape function product mp_k*mp_j
##
 for k from 0 to 2 do
  phi[k,j]   := int(Lambda^k*zz,Lambda=0..QL);
##
##  Note: 3 subtriangles for QR to divide rectangle
##  0<x<LE    -e<t<1-e     LE = sq3*(1-e) e>0   LE = sq3*(1+e) e<0
##  origin at (0,0)
##
phi[k,j]   := subs(sqrt(eps^2)=eps,1/sqrt(eps^2)=1/eps,
                   sqrt(a^2)=a,1/sqrt(a^2)=1/a,
               phi[k,j]);
phi[k,j]   := subs(ln(eps)=loge,phi[k,j]);
phi[k,j]   := collect(phi[k,j],loge):
Lterm[k,j] := subs(eps=0,coeff(phi[k,j],loge));
phi[k,j]   := coeff(phi[k,j],loge,0);
phi[k,j]   := subs(eps=0,phi[k,j]);
phi[k,j]   := subs((a^2*QL^2)^(1/2)=a*QL,
                   (a^2*QL^2)^(-1/2)=1/(a*QL),
               phi[k,j]);
phi[k,j]   := normal(expand(phi[k,j]));
 od;
od;
## output
file := `c_phi_A`:
fortran(phi,filename=file);
fortran(Lterm,filename=file);