

       (* NON - 
    AXISYMMETRIC SPHERICAL HARMONICS IN SPHERICAL COORDINATES (UP TO m = 
          3 *)

(* 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, \[Beta]   *)

FO = 1
F1 = Simplify[TrigReduce[Simplify[R*P1[Cos[\[Beta]]]]]];
F2 = Simplify[TrigReduce[Simplify[R^2*P2[Cos[\[Beta]]]]]];
F3 = Simplify[TrigReduce[Simplify[R^3*P3[Cos[\[Beta]]]]]];
F4 = Simplify[TrigReduce[Simplify[R^4*P4[Cos[\[Beta]]]]]];
F5 = Simplify[TrigReduce[Simplify[R^5*P5[Cos[\[Beta]]]]]];
F6 = Simplify[TrigReduce[Simplify[R^6*P6[Cos[\[Beta]]]]]];
F7 = Simplify[TrigReduce[Simplify[R^7*P7[Cos[\[Beta]]]]]];
F8 = Simplify[TrigReduce[Simplify[R^8*P8[Cos[\[Beta]]]]]];
F9 = Simplify[TrigReduce[Simplify[R^9*P9[Cos[\[Beta]]]]]];
 

(* Harmonic axisymmetric singular potentials in terms of R, \[Beta]  *)

S0 = 1/R;
S1 = TrigReduce[Simplify[R^(-2)*P1[Cos[\[Beta]]]]];
S2 = TrigReduce[Simplify[R^(-3)*P2[Cos[\[Beta]]]]];
S3 = TrigReduce[Simplify[R^(-4)*P3[Cos[\[Beta]]]]];
S4 = TrigReduce[Simplify[R^(-5)*P4[Cos[\[Beta]]]]];
S5 = TrigReduce[Simplify[R^(-6)*P5[Cos[\[Beta]]]]];
S6 = TrigReduce[Simplify[R^(-7)*P6[Cos[\[Beta]]]]];
S7 = TrigReduce[Simplify[R^(-8)*P7[Cos[\[Beta]]]]];
S8 = TrigReduce[Simplify[R^(-9)*P8[Cos[\[Beta]]]]];
S9 = TrigReduce[Simplify[R^(-10)*P9[Cos[\[Beta]]]]];
 
(* Legendre functions P^1_i. We first cancluate dP_n/dx *)

dP1 = D[P1[x], x];
dP2 = D[P2[x], x];
dP3 = D[P3[x], x];
dP4 = D[P4[x], x];
dP5 = D[P5[x], x];
dP6 = D[P6[x], x];
dP7 = D[P7[x], x];
dP8 = D[P8[x], x];
dP9 = D[P9[x], x];
P11[x_] = -Sqrt[1 - x^2]*dP1;
P12[x_] = -Sqrt[1 - x^2]*dP2;
P13[x_] = -Sqrt[1 - x^2]*dP3;
P14[x_] = -Sqrt[1 - x^2]*dP4;
P15[x_] = -Sqrt[1 - x^2]*dP5;
P16[x_] = -Sqrt[1 - x^2]*dP6;
P17[x_] = -Sqrt[1 - x^2]*dP7;
P18[x_] = -Sqrt[1 - x^2]*dP8;
P19[x_] = -Sqrt[1 - x^2]*dP9;
(* The following functions will yield harmonic potentials if multiplied by \
either Sin[\[Theta]] or Cos[\[Theta]] *)

(* Bounded potentials *)

F11 = Simplify[R^1*P11[Cos[\[Beta]]]];
F12 = Simplify[R^2*P12[Cos[\[Beta]]]];
F13 = Simplify[R^3*P13[Cos[\[Beta]]]];
F14 = Simplify[R^4*P14[Cos[\[Beta]]]];
F15 = Simplify[R^5*P15[Cos[\[Beta]]]];
F16 = Simplify[R^6*P16[Cos[\[Beta]]]];
F17 = Simplify[R^7*P17[Cos[\[Beta]]]];
F18 = Simplify[R^8*P18[Cos[\[Beta]]]];
F19 = Simplify[R^9*P19[Cos[\[Beta]]]];

(* Singular potentials *)

S11 = Simplify[R^(-2)*P11[Cos[\[Beta]]]];
S12 = Simplify[R^(-3)*P12[Cos[\[Beta]]]];
S13 = Simplify[R^(-4)*P13[Cos[\[Beta]]]];
S14 = Simplify[R^(-5)*P14[Cos[\[Beta]]]];
S15 = Simplify[R^(-6)*P15[Cos[\[Beta]]]];
S16 = Simplify[R^(-7)*P16[Cos[\[Beta]]]];
S17 = Simplify[R^(-8)*P17[Cos[\[Beta]]]];
S18 = Simplify[R^(-9)*P18[Cos[\[Beta]]]];
S19 = Simplify[R^(-10)*P19[Cos[\[Beta]]]];

(* Repeat for P^2_i *)

d2P2 = D[dP2, x];
d2P3 = D[dP3, x];
d2P4 = D[dP4, x];
d2P5 = D[dP5, x];
d2P6 = D[dP6, x];
d2P7 = D[dP7, x];
d2P8 = D[dP8, x];
d2P9 = D[dP9, x];

P22[x_] = (1 - x^2)*d2P2;
P23[x_] = (1 - x^2)*d2P3;
P24[x_] = (1 - x^2)*d2P4;
P25[x_] = (1 - x^2)*d2P5;
P26[x_] = (1 - x^2)*d2P6;
P27[x_] = (1 - x^2)*d2P7;
P28[x_] = (1 - x^2)*d2P8;
P29[x_] = (1 - x^2)*d2P9;


(* The following functions will yield harmonic potentials if multiplied by \
either Sin[2*\[Theta]] or Cos[2*\[Theta]] *)

(* Bounded potentials *)

F22 = Simplify[R^2*P22[Cos[\[Beta]]]];
F23 = Simplify[R^3*P23[Cos[\[Beta]]]];
F24 = Simplify[R^4*P24[Cos[\[Beta]]]];
F25 = Simplify[R^5*P25[Cos[\[Beta]]]];
F26 = Simplify[R^6*P26[Cos[\[Beta]]]];
F27 = Simplify[R^7*P27[Cos[\[Beta]]]];
F28 = Simplify[R^8*P28[Cos[\[Beta]]]];
F29 = Simplify[R^9*P29[Cos[\[Beta]]]];

(* Singular potentials *)

S22 = Simplify[R^(-3)*P22[Cos[\[Beta]]]];
S23 = Simplify[R^(-4)*P23[Cos[\[Beta]]]];
S24 = Simplify[R^(-5)*P24[Cos[\[Beta]]]];
S25 = Simplify[R^(-6)*P25[Cos[\[Beta]]]];
S26 = Simplify[R^(-7)*P26[Cos[\[Beta]]]];
S27 = Simplify[R^(-8)*P27[Cos[\[Beta]]]];
S28 = Simplify[R^(-9)*P28[Cos[\[Beta]]]];
S29 = Simplify[R^(-10)*P29[Cos[\[Beta]]]];

(* Repeat for P^3_i *)

d3P3 = D[d2P3, x];
d3P4 = D[d2P4, x];
d3P5 = D[d2P5, x];
d3P6 = D[d2P6, x];
d3P7 = D[d2P7, x];
d3P8 = D[d2P8, x];
d3P9 = D[d2P9, x];


P33[x_] = -(1 - x^2)^(3/2)*d3P3;
P34[x_] = -(1 - x^2)^(3/2)*d3P4;
P35[x_] = -(1 - x^2)^(3/2)*d3P5;
P36[x_] = -(1 - x^2)^(3/2)*d3P6;
P37[x_] = -(1 - x^2)^(3/2)*d3P7;
P38[x_] = -(1 - x^2)^(3/2)*d3P8;
P39[x_] = -(1 - x^2)^(3/2)*d3P9;

(* The following functions will yield harmonic potentials if multiplied by \
either Sin[3*\[Theta]] or Cos[3*\[Theta]] *)

(* Bounded potentials *)

F33 = Simplify[R^3*P33[Cos[\[Beta]]]];
F34 = Simplify[R^4*P34[Cos[\[Beta]]]];
F35 = Simplify[R^5*P35[Cos[\[Beta]]]];
F36 = Simplify[R^6*P36[Cos[\[Beta]]]];
F37 = Simplify[R^7*P37[Cos[\[Beta]]]];
F38 = Simplify[R^8*P38[Cos[\[Beta]]]];
F39 = Simplify[R^9*P39[Cos[\[Beta]]]];
(* Singular potentials *)

S33 = Simplify[R^(-4)*P33[Cos[\[Beta]]]];
S34 = Simplify[R^(-5)*P34[Cos[\[Beta]]]];
S35 = Simplify[R^(-6)*P35[Cos[\[Beta]]]];
S36 = Simplify[R^(-7)*P36[Cos[\[Beta]]]];
S37 = Simplify[R^(-8)*P37[Cos[\[Beta]]]];
S38 = Simplify[R^(-9)*P38[Cos[\[Beta]]]];
S39 = Simplify[R^(-10)*P39[Cos[\[Beta]]]];

