##
##   Symmetric Galerkin 3D Laplace
##   Coincident LINEAR 
##
##   Log term
##
##
with(linalg):
readlib(fortran):
##
## define arrays
##
xp := array(1..3): yp := array(1..3): zp := array(1..3):
N  := array(1..3): 
LL := array(1..1):
##
## output from Hyp_cc_A
##
sq3 := sqrt(3):
##  (2/3), (1/3) from shape functions omitted
Lterm := 
   jp^2*sin(p)/(a_cc*cos(p)^2+a_ss*sin(p)^2-a_cs*sin(p)*cos(p))^(3/2);
##
## q = cotan(p)   dq = -1/sp^2 * dp   sp^2 factor already in there
##  switching [-inf,inf] limits gets rid of minus sign
##
Lterm := 
   jp^2/(a_cc*cotan(p)^2+a_ss-a_cs*cotan(p))^(3/2)/sin(p)^2;
Lterm := 
   jp^2/(acc*q^2+ass-acs*q)^(3/2);
#
##
## WLOG
##
xp[1] := 0: yp[1] := 0: zp[1] := 0:
            yp[2] := 0: zp[2] := 0:
                        zp[3] := 0:
N[1] := 0:  N[2] := 0: N[3] := 1:
jp   := xp[2]*yp[3]/(2*sq3):
##
acc := yp[2]**2/4+yp[1]**2/4+xp[2]**2/4-xp[1]*xp[2]/2+xp[1]**2/4
     +zp[2]**2/4-yp[1]*yp[2]/2-zp[1]*zp[2]/2+zp[1]**2/4:
acs := (-xp[1]*xp[3]-xp[2]**2/2+zp[1]**2/2+xp[1]**2/2-yp[2]**2/2
     +zp[3]*zp[2]+xp[3]*xp[2]+yp[1]**2/2-yp[1]*yp[3]+yp[3]*yp[2]
     -zp[2]**2/2-zp[1]*zp[3])/sq3:
ass := (-yp[3]*yp[2]-xp[3]*xp[2]-xp[1]*xp[3]+xp[1]*xp[2]/2+zp[1]**2/4
        +xp[2]**2/4+zp[2]**2/4+xp[1]**2/4+xp[3]**2-zp[1]*zp[3]-yp[1]*yp[3]
        +yp[1]**2/4+yp[3]**2-zp[3]*zp[2]+zp[3]**2+yp[2]**2/4
        +yp[1]*yp[2]/2+zp[1]*zp[2]/2)/3:
acc := expand(acc);
acs := expand(acs);
ass := expand(ass);
##
######
LL[1]   := int(Lterm,q=-infinity..infinity):
LL[1]   := subs((xp[2]^2)^(-1/2)=1/xp[2],LL[1]);
##########
# output
file := `C_Log_A`:
fortran(LL,filename=file);