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

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

(* Harmonic axisymmetric singular potentials in terms of R, β  *)

S0=1/R;
S1=TrigReduce[Simplify[R^(-2)*P1[Cos[β]]]];
S2=TrigReduce[Simplify[R^(-3)*P2[Cos[β]]]];
S3=TrigReduce[Simplify[R^(-4)*P3[Cos[β]]]];
S4=TrigReduce[Simplify[R^(-5)*P4[Cos[β]]]];
S5=TrigReduce[Simplify[R^(-6)*P5[Cos[β]]]];
S6=TrigReduce[Simplify[R^(-7)*P6[Cos[β]]]];
S7=TrigReduce[Simplify[R^(-8)*P7[Cos[β]]]];
S8=TrigReduce[Simplify[R^(-9)*P8[Cos[β]]]];
S9=TrigReduce[Simplify[R^(-10)*P9[Cos[β]]]];

(* 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[θ] or Cos[θ] *)

(* Bounded potentials *)

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

(* Singular potentials *)

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

(* 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*θ] or Cos[2*θ] *)

(* Bounded potentials *)

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

(* Singular potentials *)

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

(* 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*θ] or Cos[3*θ] *)

(* Bounded potentials *)

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

(* Singular potentials *)

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