#				spsing

#This file generates harmonic singular functions of the form R^nPn(cos(beta)*ln(R+R*cos(beta))+Wn(R,beta). 
#These functions are typically used for problems of the solid or hollow cone. 
#To use them, uncomment lines in each block up to and including the required function. 
#Higher order functions can be generated using a similar procedure.
#
#
#Legendre polynomials as functions of x
#
P0:=1:
P01:=x:
P02:=simplify((3*x*P01-1*P0)/2):
P03:=simplify((5*x*P02-2*P01)/3):
#P04:=simplify((7*x*P03-3*P02)/4):
#P05:=simplify((9*x*P04-4*P03)/5):
#P06:=simplify((11*x*P05-5*P04)/6):
#P07:=simplify((13*x*P06-6*P05)/7):
#P08:=simplify((15*x*P07-7*P06)/8):
#P09:=simplify((17*x*P08-8*P07)/9):
#
#Simplified forms of Legendre polynomials.
#
P1:=unapply(P01,x):
P2:=unapply(simplify(P02),x):
P3:=unapply(simplify(P03),x):
#P4:=unapply(simplify(P04),x):
#P5:=unapply(simplify(P05),x):
#P6:=unapply(simplify(P06),x):
#P7:=unapply(simplify(P07),x):
#P8:=unapply(simplify(P08),x):
#P9:=unapply(simplify(P09),x):
#
R1:=sqrt(r^2+z^2):
#
#Harmonic axisymmetric bounded potentials in terms of r,z
#
F1:=expand(simplify(R1*P1(z/R1))):
F2:=expand(simplify(R1^2*P2(z/R1))):
F3:=expand(simplify(R1^3*P3(z/R1))):
#F4:=expand(simplify(R1^4*P4(z/R1))):
#F5:=expand(simplify(R1^5*P5(z/R1))):
#F6:=expand(simplify(R1^6*P6(z/R1))):
#F7:=expand(simplify(R1^7*P7(z/R1))):
#F8:=expand(simplify(R1^8*P8(z/R1))):
#F9:=expand(simplify(R1^9*P9(z/R1))):
#
#Functions of the form R1^n*ln(R1+z)+W_n.
#
#Only even powers of r can occur.
#Terms with (for example) R1^3*z^2 are implicitly included with R1*z^4 and R1*z^2*r^2.
#There exists a bounded polynomial with any power of z (z^n) as a leading term, so no term with just z^n is needed. 
#More generally, we can add in an arbitrary multiplier of Fn as required.
#
W0:=ln(R1+z):
W1:=z*ln(R1+z)-R1:
W2:=F2*ln(R1+z)+C21*R1*z+C22*r^2:
W3:=F3*ln(R1+z)+C31*R1*z^2+C32*R1*r^2+C33*z*r^2:
#W4:=F4*ln(R1+z)+C41*R1*z^3+C42*R1*z*r^2+C43*z^2*r^2+C44*r^4:
#W5:=F5*ln(R1+z)+C51*R1*z^4+C52*R1*z^2*r^2+C53*R1*r^4+C54*z^3*r^2+C55*z*r^4:
#W6:=F6*ln(R1+z)+C61*R1*z^5+C62*R1*z^3*r^2+C63*R1*z*r^4+C64*z^4*r^2
#+C65*z^2*r^4+C66*r^6:
#W7:=F7*ln(R1+z)+C71*R1*z^6+C72*R1*z^4*r^2+C73*R1*z^2*r^4+C74*R1*r^6
#+C75*z^5*r^2+C76*z^3*r^4+C77*z*r^6:
#W8:=F8*ln(R1+z)+C81*R1*z^7+C82*R1*z^5*r^2+C83*R1*z^3*r^4+C84*R1*z*r^6
#+C85*z^6*r^2+C86*z^4*r^4+C87*z^2*r^6+C88*r^8:
#W9:=F9*ln(R1+z)+C91*R1*z^8+C92*R1*z^6*r^2+C93*R1*z^4*r^4+C94*R1*z^2*r^6
#+C95*R1*r^8+C96*z^7*r^2+C97*z^5*r^4+C98*z^3*r^6+C99*z*r^8:
#
#Apply the Laplace operator
#
#
DW1:=simplify(diff(W1,r,r)+(1/r)*diff(W1,r)+diff(W1,z,z)):
DW2:=unapply(simplify(diff(W2,r,r)+(1/r)*diff(W2,r)+diff(W2,z,z))*(R1+z)^2*R1,r,z):
DW3:=unapply(simplify(diff(W3,r,r)+(1/r)*diff(W3,r)+diff(W3,z,z))*(R1+z)^2*R1*r,r,z):
#DW4:=unapply(simplify(diff(W4,r,r)+(1/r)*diff(W4,r)+diff(W4,z,z))*(R1+z)^2*R1,r,z):
#DW5:=unapply(simplify(diff(W5,r,r)+(1/r)*diff(W5,r)+diff(W5,z,z))*(R1+z)^2*R1,r,z):
#DW6:=unapply(simplify(diff(W6,r,r)+(1/r)*diff(W6,r)+diff(W6,z,z))*(R1+z)^2*R1,r,z):
#DW7:=unapply(simplify(diff(W7,r,r)+(1/r)*diff(W7,r)+diff(W7,z,z))*(R1+z)^2*R1,r,z):
#DW8:=unapply(simplify(diff(W8,r,r)+(1/r)*diff(W8,r)+diff(W8,z,z))*(R1+z)^2*R1,r,z):
#DW9:=unapply(simplify(diff(W9,r,r)+(1/r)*diff(W9,r)+diff(W9,z,z))*(R1+z)^2*R1,r,z):
#
#Choose constants to make functions harmonic.
#
#solution1:=solve({coeff(DW1,r,0)=0},{C11}):
solution2:=solve({DW2(1,0)=0,DW2(0,1)=0},{C21,C22}):
solution3:=solve({DW3(1,0)=0,DW3(2,1)=0,DW3(1,2)=0},
{C31,C32,C33}):
#solution4:=solve({DW4(1,0)=0,DW4(2,1)=0,DW4(1,2)=0,DW4(3,1)=0},
#{C41,C42,C43,C44}):
#solution5:=solve({DW5(1,0)=0,DW5(2,1)=0,DW5(1,2)=0,DW5(3,1)=0,DW5(1,3)=0},
#{C51,C52,C53,C54,C55}):
#solution6:=solve({DW6(1,0)=0,DW6(2,1)=0,DW6(1,2)=0,DW6(3,1)=0,
#DW6(1,3)=0,DW6(1,4)=0},{C61,C62,C63,C64,C65,C66}):
#solution7:=solve({DW7(1,0)=0,DW7(2,1)=0,DW7(1,2)=0,DW7(3,1)=0,
#DW7(1,3)=0,DW7(0,1)=0,DW7(1,1)=0},{C71,C72,C73,C74,C75,C76,C77}):
#solution8:=solve({DW8(1,0)=0,DW8(2,1)=0,DW8(1,2)=0,DW8(3,1)=0,
#DW8(1,3)=0,DW8(0,1)=0,DW8(1,1)=0,DW8(1,4)=0},{C81,C82,C83,C84,C85,C86,C87,C88}):
#solution9:=solve({DW9(1,0)=0,DW9(2,1)=0,DW9(1,2)=0,DW9(3,1)=0,
#DW9(1,3)=0,DW9(0,1)=0,DW9(1,1)=0,DW9(1,4)=0,DW9(4,1)=0},
#{C91,C92,C93,C94,C95,C96,C97,C98,C99}):
#
#Singular axisymmetric harmonic functions with ln(R1+z) multipliers
#
S0:=W0;
S1:=W1;
S2:=subs(solution2, W2);
S3:=subs(solution3, W3);
#S4:=subs(solution4, W4):
#S5:=subs(solution5, W5):
#S6:=subs(solution6, W6):
#S7:=subs(solution7, W7):
#S8:=subs(solution8, W8):
#S9:=subs(solution9, W9):
#
#For verification of harmonicity, see `sing'

F00:=subs(sqrt(r^2+z^2)=R,S0);
F01:=subs(sqrt(r^2+z^2)=R,S1);
F02:=subs(sqrt(r^2+z^2)=R,S2);
F03:=subs(sqrt(r^2+z^2)=R,S3);
#F04:=subs(sqrt(r^2+z^2)=R,S4);
#F05:=subs(sqrt(r^2+z^2)=R,S5);
#F06:=subs(sqrt(r^2+z^2)=R,S6);
#F07:=subs(sqrt(r^2+z^2)=R,S7);
#F08:=subs(sqrt(r^2+z^2)=R,S8);
#F09:=subs(sqrt(r^2+z^2)=R,S9);


F0:=subs({r=R*sin(beta),z=R*cos(beta)},F00);
F1:=subs({r=R*sin(beta),z=R*cos(beta)},F01);
F2:=subs({r=R*sin(beta),z=R*cos(beta)},F02);
F3:=subs({r=R*sin(beta),z=R*cos(beta)},F03);
#F4:=subs({r=R*sin(beta),z=R*cos(beta)},F04);
#F5:=subs({r=R*sin(beta),z=R*cos(beta)},F05);
#F6:=subs({r=R*sin(beta),z=R*cos(beta)},F06);
#F7:=subs({r=R*sin(beta),z=R*cos(beta)},F07);
#F8:=subs({r=R*sin(beta),z=R*cos(beta)},F08);
#F9:=subs({r=R*sin(beta),z=R*cos(beta)},F09);