#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 multiplied by constants to
#remove factors of 2,8,16 etc.
#
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 F_n*ln(r)+G_n.
#
W0:=ln(r):
W1:=z*ln(r):
W2:=F2*ln(r)+C21*r^2:
W3:=F3*ln(r)+C31*z*r^2:
W4:=F4*ln(r)+C41*r^4+C42*z^2*r^2:
W5:=F5*ln(r)+C51*z*r^4+C52*z^3*r^2:
W6:=F6*ln(r)+C61*r^6+C62*z^2*r^4+C63*z^4*r^2:
W7:=F7*ln(r)+C71*z*r^6+C72*z^3*r^4+C73*z^5*r^2:
W8:=F8*ln(r)+C81*r^8+C82*z^2*r^6+C83*z^4*r^4+C84*z^6*r^2:
W9:=F9*ln(r)+C91*z*r^8+C92*z^3*r^6+C93*z^5*r^4+C94*z^7*r^2:
#
#Apply the Laplace operator
#
#
DW2:=simplify(diff(W2,r,r)+(1/r)*diff(W2,r)+diff(W2,z,z)):
DW3:=simplify(diff(W3,r,r)+(1/r)*diff(W3,r)+diff(W3,z,z)):
DW4:=simplify(diff(W4,r,r)+(1/r)*diff(W4,r)+diff(W4,z,z)):
DW5:=simplify(diff(W5,r,r)+(1/r)*diff(W5,r)+diff(W5,z,z)):
DW6:=simplify(diff(W6,r,r)+(1/r)*diff(W6,r)+diff(W6,z,z)):
DW7:=simplify(diff(W7,r,r)+(1/r)*diff(W7,r)+diff(W7,z,z)):
DW8:=simplify(diff(W8,r,r)+(1/r)*diff(W8,r)+diff(W8,z,z)):
DW9:=simplify(diff(W9,r,r)+(1/r)*diff(W9,r)+diff(W9,z,z)):
#
#Choose constants to make functions harmonic.
#
solution2:=solve({coeff(DW2,r,0)=0},{C21}):
solution3:=solve({coeff(DW3,r,0)=0},{C31}):
solution4:=solve({coeff(DW4,r,2)=0,coeff(DW4,r,0)=0},{C41,C42}):
solution5:=solve({coeff(DW5,r,2)=0,coeff(DW5,r,0)=0},{C51,C52}):
solution6:=solve({coeff(DW6,r,4)=0,coeff(DW6,r,2)=0,
coeff(DW6,r,0)=0},{C61,C62,C63}):
solution7:=solve({coeff(DW7,r,4)=0,coeff(DW7,r,2)=0,
coeff(DW7,r,0)=0},{C71,C72,C73}):
solution8:=solve({coeff(DW8,r,6)=0,coeff(DW8,r,4)=0,
coeff(DW8,r,2)=0,coeff(DW8,r,0)=0},{C81,C82,C83,C84}):
solution9:=solve({coeff(DW9,r,6)=0,coeff(DW9,r,4)=0,
coeff(DW9,r,2)=0,coeff(DW9,r,0)=0},{C91,C92,C93,C94}):
#
#Singular harmonic functions with ln(r) 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):
#