#
#This file generates harmonic singular functions of the form
#R^nPn(z/R)*ln(R+z)+Wn(r,z). These functions are typically
#used for problems of the half-space z>0.
#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);
#
R:=sqrt(r^2+z^2);
#
#Harmonic axisymmetric bounded potentials in terms of r,z
#
F1:=expand(simplify(R*P1(z/R)));
F2:=expand(simplify(R^2*P2(z/R)));
F3:=expand(simplify(R^3*P3(z/R)));
#F4:=expand(simplify(R^4*P4(z/R)));
#F5:=expand(simplify(R^5*P5(z/R)));
#F6:=expand(simplify(R^6*P6(z/R)));
#F7:=expand(simplify(R^7*P7(z/R)));
#F8:=expand(simplify(R^8*P8(z/R)));
#F9:=expand(simplify(R^9*P9(z/R)));
#
#Functions of the form R^n*ln(R+z)+W_n.
#
#Only even powers of r can occur.
#Terms with (for example) R^3*z^2 are implicitly included with R*z^4 
#and R*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(R+z);
W1:=z*ln(R+z)-R;
W2:=F2*ln(R+z)+C21*R*z+C22*r^2;
W3:=F3*ln(R+z)+C31*R*z^2+C32*R*r^2+C33*z*r^2;
#W4:=F4*ln(R+z)+C41*R*z^3+C42*R*z*r^2+C43*z^2*r^2+C44*r^4;
#W5:=F5*ln(R+z)+C51*R*z^4+C52*R*z^2*r^2+C53*R*r^4+C54*z^3*r^2+C55*z*r^4;
#W6:=F6*ln(R+z)+C61*R*z^5+C62*R*z^3*r^2+C63*R*z*r^4+C64*z^4*r^2
#+C65*z^2*r^4+C66*r^6;
#W7:=F7*ln(R+z)+C71*R*z^6+C72*R*z^4*r^2+C73*R*z^2*r^4+C74*R*r^6
#+C75*z^5*r^2+C76*z^3*r^4+C77*z*r^6;
#W8:=F8*ln(R+z)+C81*R*z^7+C82*R*z^5*r^2+C83*R*z^3*r^4+C84*R*z*r^6
#+C85*z^6*r^2+C86*z^4*r^4+C87*z^2*r^6+C88*r^8;
#W9:=F9*ln(R+z)+C91*R*z^8+C92*R*z^6*r^2+C93*R*z^4*r^4+C94*R*z^2*r^6
#+C95*R*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))*(R+z)^2*R,r,z);
DW3:=unapply(simplify(diff(W3,r,r)+(1/r)*diff(W3,r)+diff(W3,z,z))*(R+z)^2*R*r,r,z);
#DW4:=unapply(simplify(diff(W4,r,r)+(1/r)*diff(W4,r)+diff(W4,z,z))*(R+z)^2*R,r,z);
#DW5:=unapply(simplify(diff(W5,r,r)+(1/r)*diff(W5,r)+diff(W5,z,z))*(R+z)^2*R,r,z);
#DW6:=unapply(simplify(diff(W6,r,r)+(1/r)*diff(W6,r)+diff(W6,z,z))*(R+z)^2*R,r,z);
#DW7:=unapply(simplify(diff(W7,r,r)+(1/r)*diff(W7,r)+diff(W7,z,z))*(R+z)^2*R,r,z);
#DW8:=unapply(simplify(diff(W8,r,r)+(1/r)*diff(W8,r)+diff(W8,z,z))*(R+z)^2*R,r,z);
#DW9:=unapply(simplify(diff(W9,r,r)+(1/r)*diff(W9,r)+diff(W9,z,z))*(R+z)^2*R,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(R+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);
#
#Verify that they are harmonic.
#
#dW11:=simplify(diff(S1,r,r)+(1/r)*diff(S1,r)+diff(S1,z,z));
#dW21:=simplify(diff(S2,r,r)+(1/r)*diff(S2,r)+diff(S2,z,z));
#dW31:=simplify(diff(S3,r,r)+(1/r)*diff(S3,r)+diff(S3,z,z));
#dW41:=simplify(diff(S4,r,r)+(1/r)*diff(S4,r)+diff(S4,z,z));
#dW51:=simplify(diff(S5,r,r)+(1/r)*diff(S5,r)+diff(S5,z,z));
#dW61:=simplify(diff(S6,r,r)+(1/r)*diff(S6,r)+diff(S6,z,z));
#dW71:=simplify(diff(S7,r,r)+(1/r)*diff(S7,r)+diff(S7,z,z));
#dW81:=simplify(diff(S8,r,r)+(1/r)*diff(S8,r)+diff(S8,z,z));
#dW91:=simplify(diff(S9,r,r)+(1/r)*diff(S9,r)+diff(S9,z,z));