(*           SPHERICAL HARMONICS IN SPHERICAL POLAR 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)];

Null

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=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[β]]]];