(*  AXISYMMETRIC HARMONIC POTENTIALS SINGULAR ON THE NEGATIVE z-AXIS   *)

(* 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 accelerate evaluation, only potentials for N less than 5 are evaluated automatically. To extend to higher orders, uncomment appropriate lines below. *)

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

R=Sqrt[(r^2+z^2)];
F0=1;
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=Log[R+z];
W1=z*Log[R+z]-R;
W2=F2*Log[R+z]+C21*R*z+C22*r^2;
W3=F3*Log[R+z]+C31*R*z^2+C32*R*r^2+C33*z*r^2;
W4=F4*Log[R+z]+C41*R*z^3+C42*R*z*r^2+C43*z^2*r^2+C44*r^4;
W5=F5*Log[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*Log[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*Log[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*Log[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*Log[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[D[W1,r,r]+(1/r)*D[W1,r]+D[W1,z,z]];
DW2[r_,z_]=Simplify[(D[W2,r,r]+(1/r)*D[W2,r]+D[W2,z,z])*(R+z)^2*R];
DW3[r_,z_]=Simplify[(D[W3,r,r]+(1/r)*D[W3,r]+D[W3,z,z])*(R+z)^2*R];
DW4[r_,z_]=Simplify[(D[W4,r,r]+(1/r)*D[W4,r]+D[W4,z,z])*(R+z)^2*R];
(*DW5[r_,z_]=Simplify[(D[W5,r,r]+(1/r)*D[W5,r]+D[W5,z,z])*(R+z)^2*R];
DW6[r_,z_]=Simplify[(D[W6,r,r]+(1/r)*D[W6,r]+D[W6,z,z])*(R+z)^2*R];
DW7[r_,z_]=Simplify[(D[W7,r,r]+(1/r)*D[W7,r]+D[W7,z,z])*(R+z)^2*R];
DW8[r_,z_]=Simplify[(D[W8,r,r]+(1/r)*D[W8,r]+D[W8,z,z])*(R+z)^2*R];DW9[r_,z_]=Simplify[(D[W9,r,r]+(1/r)*D[W9,r]+D[W9,z,z])*(R+z)^2*R];*)
Null

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

C21=Det[{C21/.solution2}];
C22=Det[{C22/.solution2}];
C31=Det[{C31/.solution3}];
C32=Det[{C32/.solution3}];
C33=Det[{C33/.solution3}];
C41=Det[{C41/.solution4}];
C42=Det[{C42/.solution4}];
C43=Det[{C43/.solution4}];
C44=Det[{C44/.solution4}];
(*C51=Det[{C51/.solution5}];
C52=Det[{C52/.solution5}];
C53=Det[{C53/.solution5}];
C54=Det[{C54/.solution5}];
C55=Det[{C55/.solution5}];
C61=Det[{C61/.solution6}];
C62=Det[{C62/.solution6}];
C63=Det[{C63/.solution6}];
C64=Det[{C64/.solution6}];
C65=Det[{C65/.solution6}];
C66=Det[{C66/.solution6}];
C71=Det[{C71/.solution7}];
C72=Det[{C72/.solution7}];
C73=Det[{C73/.solution7}];
C74=Det[{C74/.solution7}];
C75=Det[{C75/.solution7}];
C76=Det[{C76/.solution7}];
C77=Det[{C77/.solution7}];
C81=Det[{C81/.solution8}];
C82=Det[{C82/.solution8}];
C83=Det[{C83/.solution8}];
C84=Det[{C84/.solution8}];
C85=Det[{C85/.solution8}];
C86=Det[{C86/.solution8}];
C87=Det[{C87/.solution8}];
C88=Det[{C88/.solution8}];
C91=Det[{C91/.solution9}];
C92=Det[{C92/.solution9}];
C93=Det[{C93/.solution9}];
C94=Det[{C94/.solution9}];
C95=Det[{C95/.solution9}];
C96=Det[{C96/.solution9}];
C97=Det[{C97/.solution9}];
C98=Det[{C98/.solution9}];
C99=Det[{C99/.solution9}];*)

C21=Simplify[C21];
C22=Simplify[C22];
C31=Simplify[C31];
C32=Simplify[C32];
C33=Simplify[C33];
C41=Simplify[C41];
C42=Simplify[C42];
C43=Simplify[C43];
C44=Simplify[C44];
(*C51=Simplify[C51];
C52=Simplify[C52];
C53=Simplify[C53];
C54=Simplify[C54];
C55=Simplify[C55];
C61=Simplify[C61];
C62=Simplify[C62];
C63=Simplify[C63];
C64=Simplify[C64];
C65=Simplify[C65];
C66=Simplify[C66];
C71=Simplify[C71];
C72=Simplify[C72];
C73=Simplify[C73];
C74=Simplify[C74];
C75=Simplify[C75];
C76=Simplify[C76];
C77=Simplify[C77];
C81=Simplify[C81];
C82=Simplify[C82];
C83=Simplify[C83];
C84=Simplify[C84];
C85=Simplify[C85];
C86=Simplify[C86];
C87=Simplify[C87];
C88=Simplify[C88];
C91=Simplify[C91];
C92=Simplify[C92];
C93=Simplify[C93];
C94=Simplify[C94];
C95=Simplify[C95];
C96=Simplify[C96];
C97=Simplify[C97];
C98=Simplify[C98];
C99=Simplify[C99]*)

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