(*The foundation and development of this code can be found in Manuel Forner's book:<<Desarrollos en serie de los productos de algunas funciones especiales>>paper support:http://www.tenda.uji.es/pls/iglu/!GCPPA00.GCPPR0002?lg=CA&id_art=1110 e-Book-PDF:http://www.tenda.uji.es/pls/iglu/!GCPPA00.GCPPR0002?lg=CA&id_art=1171*)

(*DEFINITION AND NORM ASSOCIATED LEGENDRE FUNCTIONS*)

funAsLe[n_,m_,x_]:=(((-1))^m)LegendreP[n,m,x]

normafunAsLe[n_,m_]:=((2(n+m)!)/((2n+1)(n-m)!))

invnormafunAsLe[n_,m_]:=(((2n+1)(n-m)!)/(2(n+m)!))

(*INTEGRAL OF THE PRODUCT OF THREE ASSOCIATED LEGENDRE FUNCTIONS*)

solInt[α_,z_]:=If[Mod[z,2]==0,((1+(((-1))^α))/2)(((α-1)!! ((z/2))! (2^((z/2)+1)))/((α+z+1)!!)),((1+(((-1))^α))/2)(((α-1)!! z!! π)/((α+z+1)!!))]

suma1[α_,n_,m_,j_]:=Module[{t,β,γ},γ=j;β=If[j-α<0,0,j-α];Sum[(((-1))^t)((j!)/(t!(j-t)!))(((n-m-j)!)/((n-m-α-t)!(α+t-j)!)),{t,β,γ}]]

suma2[α_,n_,m_]:=Expand[Sum[((suma1[α,n,m,j])/(j!(n-m-j)!(m+j)!(n-j)!)),{j,0,n-m}]]

χ[ξ_,n_,m_,r_,s_]:=Module[{u,δ,l},δ=If[r-s-ξ<0,Abs[r-s-ξ],0];l=If[ξ>n-m,n-m,ξ];Sum[suma2[u,n,m] suma2[ξ-u,r,s],{u,δ,l}]]

χχ[ζ_,n_,m_,r_,s_,i_,z_]:=Module[{uu,δδ,ll},δδ=If[i-z-ζ<0,Abs[i-z-ζ],0];ll=If[ζ>n+r-m-s,n+r-m-s,ζ];Sum[χ[uu,n,m,r,s] suma2[ζ-uu,i,z],{uu,δδ,ll}]]

η[j_Integer]:=If[j<0,(((-1))^(Abs[j])),1]

intprodFunAsLeg[n_Integer,m_Integer,r_Integer,s_Integer,p_Integer,q_Integer]:=Which[(n≥Abs[m])&&(r≥Abs[s])&&(p≥Abs[q]),(((n+m)!n!(r+s)!r!(p+q)!p!)/(2^(n+r+p)))η[m]η[s]η[q]Sum[χχ[α,n,Abs[m],r,Abs[s],p,Abs[q]]solInt[α,Abs[m]+Abs[s]+Abs[q]],{α,0,n+r+p-(Abs[m]+Abs[s]+Abs[q])}],n<0||r<0||p<0,Print["Neither n nor r nor p can be negative numbers"],n<Abs[m]||r<Abs[s]||p<Abs[q],Print["Neither n nor r nor p can be smaller |m| , |s| and |q|, respectively"]]

(*DEVELOPMENT,IN ASSOCIATED LEGENDRE FUNCTIONS SERIES,OF THE PRODUCT OF ASSOCIATED LEGENDRE FUNCTIONS*)

producteFA[n_,m_,r_,s_]:=Simplify[funAsLe[n,m,x]funAsLe[r,s,x]]

(*Possibility 1*)(*Used for the development of the complex form*)

(*(Subscript[P, n,m])(x) (Subscript[P, r,s])(x)=(∑_(j=|m+s|)^(n+r))(Subscript[A, j])(Subscript[P, j,m+s])(x), on (Subscript[A, j])=(((∫_(-1)^1)(Subscript[P, n,m])(x) (Subscript[P, r,s])(x) (Subscript[P, j,m+s])(x)\[DifferentialD]x)/((∫_(-1)^1)((((Subscript[P, j,m+s])(x)))^2)\[DifferentialD]x)), n, r ∈ (Natural numbers ⋃ {0}), m, s ∈ Integers numbers, |m|≤n, |s|≤r*)

prodFunAsLeg[n_,m_,r_,s_]:=Sum[intprodFunAsLeg[n,m,r,s,j,m+s]invnormafunAsLe[j,m+s]funAsLe[j,m+s,x],{j,Abs[m+s],n+r}]

desprodFunAsLeg[n_,m_,r_,s_]:=Sum[intprodFunAsLeg[n,m,r,s,j,m+s]invnormafunAsLe[j,m+s]PF[j,m+s],{j,Abs[m+s],n+r}]

(*This form covers all the possible products of associated Legendre functions.But for the development of the spherical harmonic products in real form,it is necessary to look for other possibilities of development of the product of associated Legendre functions,based on the sign of the second subscript of associated Legendre function.Below some of these possibilities are implemented*)

(*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*)

(*Possibility 2*)

(*(Subscript[P, n,m])(x) (Subscript[P, r,s])(x)=(∑_(j=|m-s|)^(n+r))(Subscript[B, j])(Subscript[P, j,|m-s|])(x), where (Subscript[B, j])=(((∫_(-1)^1)(Subscript[P, n,m])(x) (Subscript[P, r,s])(x) (Subscript[P, j,|m-s|])(x)\[DifferentialD]x)/((∫_(-1)^1)((((Subscript[P, j,|m-s|])(x)))^2)\[DifferentialD]x)), n, r ∈ (Natural numbers ⋃ {0}), m, s ∈ (Natural numbers ⋃ {0}), m≤n, s≤r*)

prodFunAsLeg2[n_,m_,r_,s_]:=Sum[intprodFunAsLeg[n,m,r,s,j,Abs[m-s]]invnormafunAsLe[j,Abs[m-s]]funAsLe[j,Abs[m-s],x],{j,Abs[m-s],n+r}]

(*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*)

(*Possibility 3*)

(*(Subscript[P, n,m])(x) (Subscript[P, r,s])(x)=(∑_(j=m-s)^(n+r))(Subscript[A, j])(Subscript[P, j,s-m])(x), where (Subscript[A, j])=(((∫_(-1)^1)(Subscript[P, n,m])(x) (Subscript[P, r,s])(x) (Subscript[P, j,s-m])(x)\[DifferentialD]x)/((∫_(-1)^1)((((Subscript[P, j,s-m])(x)))^2)\[DifferentialD]x)),(Subscript[P, n,m])(x) (Subscript[P, r,s])(x)=(∑_(j=|m+s|)^(n+r))(Subscript[B, j])(Subscript[P, j,-|m+s|])(x), where (Subscript[B, j])=(((∫_(-1)^1)(Subscript[P, n,m])(x) (Subscript[P, r,s])(x) (Subscript[P, j,-|m+s|])(x)\[DifferentialD]x)/((∫_(-1)^1)((((Subscript[P, j,-|m+s|])(x)))^2)\[DifferentialD]x)), n, r, m ∈ (Natural numbers ⋃ {0}), s ∈ Negative integers, m≤n, |s|≤r*)

prodFunAsLeg3A[n_,m_,r_,s_]:=Sum[intprodFunAsLeg[n,m,r,s,j,s-m]invnormafunAsLe[j,s-m]funAsLe[j,s-m,x],{j,m-s,n+r}]

prodFunAsLeg3B[n_,m_,r_,s_]:=Sum[intprodFunAsLeg[n,m,r,s,j,-Abs[m+s]]invnormafunAsLe[j,-Abs[m+s]]funAsLe[j,-Abs[m+s],x],{j,Abs[m+s],n+r}]

(*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*)