(*      SINGULAR AND BOUNDED SPHERICAL HARMONICS IN CYLINDRICAL COORDINATES            *)

(* 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)];

P1[x_]=P01;
P2[x_]=P02;
P3[x_]=P03;
P4[x_]=P04;
P5[x_]=P05;
P6[x_]=P06;
P7[x_]=P07;
P8[x_]=P08;
P9[x_]=P09;

(* HARMONIC AXISYMMETRIC BOUNDED POTENTIALS IN TERMS OF r,z *)

R1=Sqrt[(r^2+z^2)];
F0=1
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=Log[r];
W1=z*Log[r];
W2=F2*Log[r]+C21*r^2;
W3=F3*Log[r]+C31*z*r^2;
W4=F4*Log[r]+C41*r^4+C42*z^2*r^2;
W5=F5*Log[r]+C51*z*r^4+C52*z^3*r^2;
W6=F6*Log[r]+C61*r^6+C62*z^2*r^4+C63*z^4*r^2;
W7=F7*Log[r]+C71*z*r^6+C72*z^3*r^4+C73*z^5*r^2;
W8=F8*Log[r]+C81*r^8+C82*z^2*r^6+C83*z^4*r^4+C84*z^6*r^2;
W9=F9*Log[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=Expand[D[W2,r,r]+(1/r)*D[W2,r]+D[W2,z,z]];
DW3=Expand[D[W3,r,r]+(1/r)*D[W3,r]+D[W3,z,z]];
DW4=Expand[D[W4,r,r]+(1/r)*D[W4,r]+D[W4,z,z]];
DW5=Expand[D[W5,r,r]+(1/r)*D[W5,r]+D[W5,z,z]];
DW6=Expand[D[W6,r,r]+(1/r)*D[W6,r]+D[W6,z,z]];
DW7=Expand[D[W7,r,r]+(1/r)*D[W7,r]+D[W7,z,z]];
DW8=Expand[D[W8,r,r]+(1/r)*D[W8,r]+D[W8,z,z]];
DW9=Expand[D[W9,r,r]+(1/r)*D[W9,r]+D[W9,z,z]];

(* Choose constants to make functions harmonic *)

solution2=Solve[{Coefficient[DW2,r,0]==0},{C21}];
solution3=Solve[{Coefficient[DW3,r,0]==0},{C31}];
solution4=Solve[{Coefficient[DW4,r,2]==0,Coefficient[DW4,r,0]==0},{C41,C42}];
solution5=Solve[{Coefficient[DW5,r,2]==0,Coefficient[DW5,r,0]==0},{C51,C52}];
solution6=Solve[{Coefficient[DW6,r,4]==0,Coefficient[DW6,r,2]==0,
Coefficient[DW6,r,0]==0},{C61,C62,C63}];
solution7=Solve[{Coefficient[DW7,r,4]==0,Coefficient[DW7,r,2]==0,
Coefficient[DW7,r,0]==0},{C71,C72,C73}];
solution8=Solve[{Coefficient[DW8,r,6]==0,Coefficient[DW8,r,4]==0,
Coefficient[DW8,r,2]==0,Coefficient[DW8,r,0]==0},{C81,C82,C83,C84}];
solution9=Solve[{Coefficient[DW9,r,6]==0,Coefficient[DW9,r,4]==0,
Coefficient[DW9,r,2]==0,Coefficient[DW9,r,0]==0},{C91,C92,C93,C94}];

C21=Det[{C21/.solution2}];
C31=Det[{C31/.solution3}];
C41=Det[{C41/.solution4}];
C42=Det[{C42/.solution4}];
C51=Det[{C51/.solution5}];
C52=Det[{C52/.solution5}];
C61=Det[{C61/.solution6}];
C62=Det[{C62/.solution6}];
C63=Det[{C63/.solution6}];
C71=Det[{C71/.solution7}];
C72=Det[{C72/.solution7}];
C73=Det[{C73/.solution7}];
C81=Det[{C81/.solution8}];
C82=Det[{C82/.solution8}];
C83=Det[{C83/.solution8}];
C84=Det[{C84/.solution8}];
C91=Det[{C91/.solution9}];
C92=Det[{C92/.solution9}];
C93=Det[{C93/.solution9}];
C94=Det[{C94/.solution9}];

S0=W0
S1=W1
S2=W2
S3=W3
S4=W4
S5=W5
S6=W6
S7=W7
S8=W8
S9=W9