##
##   Symmetric Galerkin 3D Laplace
##   Edge Adjacent LINEAR 
##
##
with(linalg):
readlib(fortran):
##
## define arrays
##
phi   := array(0..2);
sng   := array(0..2);
##
rh  := (eps^2 + eps*a1*Lambda + a2*Lambda^2)^(1/2);
##
## kernel
##  Note: cos(psi) factor from jacobian product omitted
##
hyp :=  Lambda^2*(ndN/rh^3 - 3*(j1p*Lambda-jp*eps)*(j1q*Lambda-eps*ndN/jp)/rh^5);
## powers from shape function product
for ll from 0 to 2 do
 phi[ll] := Lambda^ll*hyp;
 phi[ll] := int(phi[ll],Lambda=0..QL):
 phi[ll] := expand(phi[ll]);
 phi[ll] := subs( (eps^2)^(1/2)=eps,
                  (eps^2)^(3/2)=eps^3,
                  (eps^2)^(-1/2)=1/eps,
                  (eps^2)^(-3/2)=1/eps^3,
                  ln((a2)^(1/2)*eps*a1+2*eps*a2)=loge+logaa,
##  logaa = log(2*a2+a2s*a1)
            phi[ll]);
 phi[ll] := subs( eps=0,phi[ll]);
 phi[ll] := subs( (a2*QL^2)^(1/2)=a2s*QL,
                  (a2*QL^2)^(3/2)=a2s^3*QL^3,
                  (a2*QL^2)^(-1/2)=a2s^(-1)*QL^(-1),
                  (a2*QL^2)^(-3/2)=a2s^(-3)*QL^(-3),
                  ln((a2)^(3/2)*QL+a2s*QL*a2)=logQ,
                  ln(a2s^3*QL+a2s*QL*a2)=logQ,
            phi[ll]);
##  logQ = log(2*a2s^3*QL)
 phi[ll] := expand(phi[ll]):
 phi[ll] := collect(phi[ll],loge);
 sng[ll] := normal(coeff(phi[ll],loge,1));
 phi[ll] := normal(coeff(phi[ll],loge,0));
 phi[ll] := subs(logQ=logS-ln(2)+logaa,phi[ll]):
##  logS = log(4*a2s^3*QL/(2*a2+a2s*a1))
##         log(4*a2*QL/(2*a2s+a1))
 phi[ll] := subs(ln(QL)=logS-2*ln(2)-3*ln(a2s)+logaa,phi[ll]):
 phi[ll] := subs(ln(a2s^3*QL+a2s*QL*a2)=logS-ln(2)+logaa,phi[ll]):
 phi[ll] := subs(a2^(9/2)=a2s^9,
                 a2^(7/2)=a2s^7,
                 a2^(5/2)=a2s^5,
                 a2^(3/2)=a2s^3,
                 a2^(1/2)=a2s,
                 a2=a2s^2,
                 (a2s^2)^(9/2)=a2s^9,
                 (a2s^2)^(11/2)=a2s^11,
                 (a2s^2)^(13/2)=a2s^13,
                 (a2s^2)^(-9/2)=a2s^(-9),
                 (a2s^2)^(-11/2)=a2s^(-11),
            phi[ll]);
 phi[ll] := normal(expand(phi[ll]));
 phi[ll] := subs(ln(QL)=logS-2*ln(2)-3*ln(a2s)+logaa,phi[ll]):
 phi[ll] := normal(expand(phi[ll]));
 if ll=0 then
  num := numer(phi[ll]): num := expand(num):
  den := denom(phi[ll]):
  num := collect(num,[logS,ndN,j1p,j1q,jp]);
  phi[ll] := num/den;
 fi;
od;
###
### output
###
file := `edgeP_A`:
fortran(phi,filename=file);
fortran(sng,filename=file);