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

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