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

(*INTEGRAL CALCULATION OF THE PRODUCT OF THREE SPHERIC HARMONICS IN A REAL FORM*)

(* The integral to be calculated is:(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, p,q]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ*)

(*Spherical harmonics in real form can be denoted as:*)
(*(Subscript[Y, n,m])(θ,λ)=If[m≥0,c[n,m] (Subscript[P, n,m])(Sin[θ]) Cos[m λ],c[n,m] (Subscript[P, n,m])(Sin[θ]) Sin[|m| λ]], with c[n,m]=η[m] (Sqrt[(((2n+1)(n-m)!)/(ϵ[m] π(n+m)!))]), η[m]=If[m<0,(((-1))^(|m|)),1], ϵ[m]=If[m==0,4,2]*)

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

(*The possible integrals with respect to λ other than zero,which can be given, are:*)

(*(∫_0^(2π))Cos[m λ] Cos[s λ] Cos[q λ] \[DifferentialD]λ*)

(*m≥0 ∧ s≥0 ∧ q≥0*)

int1Bis[m_,s_,q_]:=Which[m==0&&s==0&&q==0,2π,(m==0&&s==q)||(s==0&&m==q)||(q==0&&m==s),π,(q==m+s)||(s==m+q)||(m==s+q),(π/2),(q≠m+s)&&(s≠m+q)&&(m≠s+q),0]

(*(∫_0^(2π))Cos[m λ] Sin[s λ] Sin[q λ] \[DifferentialD]λ ∨ (∫_0^(2π))Cos[s λ] Sin[m λ] Sin[q λ] \[DifferentialD]λ ∨ (∫_0^(2π))Cos[q λ] Sin[s λ] Sin[m λ] \[DifferentialD]λ*)

(*(m≥0 ∧ s<0 ∧ q<0) ∨ (s≥0 ∧ m<0 ∧ q<0) ∨ (q≥0 ∧ s<0 ∧ m<0)*)

int2[m_,s_,q_]:=Which[(s==0)||(q==0),0,m==0&&s==q,π,(q==m+s)||(s==m+q),(π/2),m==s+q,-(π/2),(q≠m+s)&&(s≠m+q)&&(m≠s+q),0]

int2Bis[m_,s_,q_]:=Which[m≤0&&s≤0&&q≥0,int2[q,Abs[m],Abs[s]],m≤0&&q≤0&&s≥0,int2[s,Abs[m],Abs[q]],q≤0&&s≤0&&m≥0,int2[m,Abs[q],Abs[s]],(m≥0&&s≥0&&q≤0)||(m≥0&&s≤0&&q≥0)||(m≤0&&s≥0&&q≥0),0,m<0&&s<0&&q<0,0]

(*The integral resulting from the above is defined below:*)

intR[m_,s_,q_]:=If[m≥0&&s≥0&&q≥0,int1Bis[m,s,q],int2Bis[m,s,q]]

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

(* (∫_(-1)^1)(x^α) (((1-(x^2)))^(z/2))\[DifferentialD]x =If[Mod[z,2]==0,((1+(((-1))^α))/2)(((α-1)!! ((z/2))! (2^((z/2)+1)))/((α+z+1)!!)),((1+(((-1))^α))/2)(((α-1)!! z!! π)/((α+z+1)!!))]*)

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)!!))]

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

(*Construction of the "n + m"-th derivative of ((((x^2)-1))^n)*)
(*((d^(n+m))/(d (x^(n+m))))((((x^2)-1))^n) = (n+m)! (((n!))^2) (∑_(α=0)^(n-m))((∑_(j=0)^(n-m))((1/(j!(n-m-j)!(m+j)!(n-j)!))((∑_(t=β)^j)(((-1))^t)(({{j}, {t}}))(({{n-m-j}, {n-m-α-t}})))))(x^α) , where β=If[j-α<0,0,j-α], "n" and "m" belong to the union of natural numbers,the zero and "m≤n"*)

(*suma1[α_,n_,m_,j_]:=(∑_(t=β)^j)(((-1))^t)(({{j}, {t}}))(({{n-m-j}, {n-m-α-t}}))*)

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_]:=(∑_(j=0)^(n-m))((1/(j!(n-m-j)!(m+j)!(n-j)!))((∑_(t=β)^j)(((-1))^t)(({{j}, {t}}))(({{n-m-j}, {n-m-α-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_]:=(∑_(u=δ)^g)((Subscript[Λ, u]) (Subscript[E, ξ-u])),where δ=If[r-s-ξ<0,Abs[r-s-ξ],0],g=If[ξ>n-m,n-m,ξ],(Subscript[Λ, u])=suma2[u,n,m],(Subscript[E, ξ-u])=suma2[ξ-u,r,s]*)

χ[ξ_,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_]:=(∑_(uu=δδ)^gg)((Subscript[X, uu]) (Subscript[H, ζ-uu])), where δδ=If[i-z-ζ<0,Abs[i-z-ζ],0],gg=If[ζ>n+r-m-s,n+r-m-s,ζ],(Subscript[X, uu])=χ[uu,n,m,r,s],(Subscript[H, ζ-uu])= suma2[ζ-uu,i,z]*)

χχ[ζ_,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}]]

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

ϵ[m_]:=If[m==0,4,2]

cR[n_,m_]:=(Sqrt[(((2n+1)(n-m)!)/(ϵ[m]π(n+m)!))])

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

(*intFornerR[j_,zz_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,zz]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ*)

intFornerR[j_Integer,zz_Integer,n_Integer,m_Integer,r_Integer,s_Integer]:=intR[m,s,zz](((n+Abs[m])!n!(r+Abs[s])!r!(j+Abs[zz])!j!)/(2^(n+r+j)))cR[n,Abs[m]] cR[r,Abs[s]] cR[j,Abs[zz]] Sum[χχ[ζ,n,Abs[m],r,Abs[s],j,Abs[zz]]solInt[ζ,Abs[m]+Abs[s]+Abs[zz]] ,{ζ,0,n+r+j-(Abs[m]+Abs[s]+Abs[zz])}]

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

(*ASCERTAINMENT OF INTEGRAL CALCULATION*)

(*Definition of associated Legendre functions*)
(*(Subscript[P, n,m])(x)=(((1-(x^2)))^(m/2)) ((d^m)/(d (x^m)))((Subscript[P, n])(x))=(1/((2^n) n!)) (((1-(x^2)))^(m/2)) ((d^(n+m))/(d (x^(n+m))))((((x^2)-1))^n)*)

funAsLe[n_,m_,θ_]:=(((-1))^m)LegendreP[n,m,Sin[θ]]

(*Definition of spherical harmonics in real form*)

YF[n_,m_,θ_,λ_]:=If[m≥0,If[m>0,(Sqrt[(((2n+1)(n-m)!)/(2π(n+m)!))])(funAsLe[n,m,θ]/.{Sqrt[1-Cos[x_]^2]->Sin[x],Sqrt[1-Sin[x_]^2]->Cos[x]})Cos[m λ],(Sqrt[((2n+1)/(4π))])(LegendreP[n,Sin[θ]]/.{Sqrt[1-Cos[x_]^2]->Sin[x],Sqrt[1-Sin[x_]^2]->Cos[x]})],(Sqrt[(((2n+1)(n-m)!)/(2π(n+m)!))])((((-1))^(Abs[m]))funAsLe[n,m,θ]/.{Sqrt[1-Cos[x_]^2]->Sin[x],Sqrt[1-Sin[x_]^2]->Cos[x]})Sin[Abs[m] λ]]

(*integralR[j_,z_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,zz]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ*)

integralR[j_,zz_,n_,m_,r_,s_]:=Integrate[Integrate[YF[j,zz,θ,λ] YF[n,m,θ,λ] YF[r,s,θ,λ] Cos[θ],{θ,-(π/2),(π/2)}],{λ,0,2π}]

(*Examples*)

integralR[5,4,6,-2,7,-2]

intFornerR[5,4,6,-2,7,-2]

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

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

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

(*GENERALIZATION OF THE DEVELOPMENT,IN SERIES OF HARMONICS IN REAL FORM,OF THE PRODUCT OF SPHERICAL HARMONICS IN REAL FORM.*)

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

(*If m>0 ∧ s>0*)

(*a11R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,m+s]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,m+s,n,m,r,s]*)

a11R[j_,n_,m_,r_,s_]:=Module[{z},z=m+s;intFornerR[j,z,n,m,r,s]]

(*b11R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,|m-s|]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,|m-s|,n,m,r,s]*)

b11R[j_,n_,m_,r_,s_]:=Module[{z},z=Abs[m-s];intFornerR[j,z,n,m,r,s]]

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

(*If m<0 ∧ s<0*)

(*a22R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,|m+s|]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,|m+s|,n,m,r,s]*)

a22R[j_,n_,m_,r_,s_]:=Module[{z},z=Abs[m+s];intFornerR[j,z,n,m,r,s]]

(*b22R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,|m-s|]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,|m-s|,n,m,r,s]*)

b22R[j_,n_,m_,r_,s_]:=Module[{z},z=Abs[m-s];intFornerR[j,z,n,m,r,s]]

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

(*If m>0 ∧ s<0*)

(*a12R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,-m+s]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,-m+s,n,m,r,s]*)

a12R[j_,n_,m_,r_,s_]:=Module[{z},z=-m+s;intFornerR[j,z,n,m,r,s]]

(*b12R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,-|m+s|]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,-|m+s|,n,m,r,s]*)

b12R[j_,n_,m_,r_,s_]:=Module[{z},z=-Abs[m+s];intFornerR[j,z,n,m,r,s]]

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

(*If m<0 ∧ s>0*)

(*a21R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,m-s]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,m-s,n,m,r,s]*)

a21R[j_,n_,m_,r_,s_]:=Module[{z},z=m-s;intFornerR[j,z,n,m,r,s]]

(*b21R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,-|m+s|]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,-|m+s|,n,m,r,s]*)

b21R[j_,n_,m_,r_,s_]:=Module[{z},z=-Abs[m+s];intFornerR[j,z,n,m,r,s]]

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

(*If m=0 ∧ s>0*)

(*a01R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,0])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,s]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,s,n,0,r,s]*)

a01R[j_,n_,m_,r_,s_]:=Module[{z},z=s;intFornerR[j,z,n,m,r,s]]

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

(*If m>0 ∧ s=0*)

(*a10R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,0])(θ,λ) (Subscript[Y, j,m]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,m,n,m,r,0]*)

a10R[j_,n_,m_,r_,s_]:=Module[{z},z=m;intFornerR[j,z,n,m,r,s]]

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

(*If m=0 ∧ s<0*)

(*a02R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,0])(θ,λ) (Subscript[Y, r,s])(θ,λ) (Subscript[Y, j,s]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,s,n,0,r,s]*)

a02R[j_,n_,m_,r_,s_]:=Module[{z},z=s;intFornerR[j,z,n,m,r,s]]

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

(*If m<0 ∧ s=0*)

(*a20R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,0])(θ,λ) (Subscript[Y, j,m]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,m,n,m,r,0]*)

a20R[j_,n_,m_,r_,s_]:=Module[{z},z=m;intFornerR[j,z,n,m,r,s]]

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

(*If m=0 ∧ s=0*)

(*a00R[j_,n_,m_,r_,s_]:=(∫_(λ=0)^(2π))(∫_(θ=-(π/2))^(π/2)) (Subscript[Y, n,0])(θ,λ) (Subscript[Y, r,0])(θ,λ) (Subscript[Y, j,0]) Cos[θ] \[DifferentialD]θ \[DifferentialD]λ=intFornerR[j,0,n,0,r,0]*)

a00R[j_,n_,m_,r_,s_]:=Module[{z},z=m;intFornerR[j,z,n,m,r,s]]

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

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

(*If m>0 ∧ s>0*)

(*prod11[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=m+s)^(n+r))a11R[i,n,m,r,s] (Subscript[Y, i,m+s])(θ,λ)+(∑_(p=|m-s|)^(n+r))b11R[p,n,m,r,s] (Subscript[Y, p,|m-s|])(θ,λ)*)

prod11[n_,m_,r_,s_]:=Sum[a11R[i,n,m,r,s]*FF[i,m+s],{i,m+s,n+r}]+Sum[b11R[p,n,m,r,s]*FF[p,Abs[m-s]],{p,Abs[m-s],n+r}]

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

(*If m<0 ∧ s<0*)

(*prod22[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=-m-s)^(n+r))a22R[i,n,m,r,s] (Subscript[Y, i,-m-s])(θ,λ)+(∑_(p=|m-s|)^(n+r))b22R[p,n,m,r,s] (Subscript[Y, p,|m-s|])(θ,λ)*)

prod22[n_,m_,r_,s_]:=Sum[a22R[i,n,m,r,s]*FF[i,-m-s],{i,-m-s,n+r}]+Sum[b22R[p,n,m,r,s]*FF[p,Abs[m-s]],{p,Abs[m-s],n+r}]

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

(*If m>0 ∧ s<0*)

(*If m+s=0*)

(*prod12[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=m-s)^(n+r))a12R[i,n,m,r,s] (Subscript[Y, i,-m+s])(θ,λ)*)

(*If m+s>0*)

(*prod12[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=m-s)^(n+r))a12R[i,n,m,r,s] (Subscript[Y, i,-m+s])(θ,λ)+(∑_(p=|m+s|)^(n+r))b12R[p,n,m,r,s] (Subscript[Y, p,-|m+s|])(θ,λ)*)

(*If m+s<0*)

(*prod12[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=m-s)^(n+r))a12R[i,n,m,r,s] (Subscript[Y, i,-m+s])(θ,λ)+(∑_(p=|m+s|)^(n+r))b12R[p,n,m,r,s] (Subscript[Y, p,-|m+s|])(θ,λ)*)

prod12[n_,m_,r_,s_]:=If[m+s==0,Sum[a12R[i,n,m,r,s]*FF[i,-m+s],{i,m-s,n+r}],If[m+s>0,Sum[a12R[i,n,m,r,s]*FF[i,-m+s],{i,m-s,n+r}]+Sum[b12R[p,n,m,r,s]*FF[p,-Abs[m+s]],{p,Abs[m+s],n+r}],Sum[a12R[i,n,m,r,s]*FF[i,-m+s],{i,m-s,n+r}]+Sum[b12R[p,n,m,r,s]*FF[p,-Abs[m+s]],{p,Abs[m+s],n+r}]]]

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

(*If m<0 ∧ s>0*)

(*If m+s=0*)

(*prod21[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=-m+s)^(n+r))a21R[i,n,m,r,s] (Subscript[Y, i,-s+m])(θ,λ)*)

(*If m+s>0*)

(*prod21[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=-m+s)^(n+r))a21R[i,n,m,r,s] (Subscript[Y, i,-s+m])(θ,λ)+(∑_(p=|m+s|)^(n+r))b21R[p,n,m,r,s] (Subscript[Y, p,-|m+s|])(θ,λ)*)

(*If m+s<0*)

(*prod21[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=-m+s)^(n+r))a21R[i,n,m,r,s] (Subscript[Y, i,-s+m])(θ,λ)+(∑_(p=|m+s|)^(n+r))b21R[p,n,m,r,s] (Subscript[Y, p,-|m+s|])(θ,λ)*)

prod21[n_,m_,r_,s_]:=If[m+s==0,Sum[a21R[i,n,m,r,s]*FF[i,-s+m],{i,-m+s,n+r}],If[m+s>0,Sum[a21R[i,n,m,r,s]*FF[i,-s+m],{i,-m+s,n+r}]+Sum[b21R[p,n,m,r,s]*FF[p,-Abs[m+s]],{p,Abs[m+s],n+r}],Sum[a21R[i,n,m,r,s]*FF[i,-s+m],{i,-m+s,n+r}]+Sum[b21R[p,n,m,r,s]*FF[p,-Abs[m+s]],{p,Abs[m+s],n+r}]]]

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

(*If m=0 ∧ s>0*)

(*prod01[n_,m_,r_,s_]:= (Subscript[Y, n,0])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=s)^(n+r))a01R[i,n,0,r,s] (Subscript[Y, i,s])(θ,λ)*)

prod01[n_,m_,r_,s_]:=Sum[ a01R[i,n,m,r,s]*FF[i,s],{i,s,n+r}]

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

(*If m>0 ∧ s=0*)

(*prod10[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,0])(θ,λ)= (∑_(i=m)^(n+r))a10R[i,n,m,r,0] (Subscript[Y, i,m])(θ,λ)*)

prod10[n_,m_,r_,s_]:=Sum[ a10R[i,n,m,r,s]*FF[i,m],{i,m,n+r}]

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

(*If m=0 ∧ s<0*)

(*prod02[n_,m_,r_,s_]:= (Subscript[Y, n,0])(θ,λ) (Subscript[Y, r,s])(θ,λ)= (∑_(i=-s)^(n+r))a02R[i,n,0,r,s] (Subscript[Y, i,s])(θ,λ)*)

prod02[n_,m_,r_,s_]:=Sum[ a02R[i,n,m,r,s]*FF[i,s],{i,-s,n+r}]

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

(*If m<0 ∧ s=0*)

(*prod20[n_,m_,r_,s_]:= (Subscript[Y, n,m])(θ,λ) (Subscript[Y, r,0])(θ,λ)= (∑_(i=-m)^(n+r))a20R[i,n,m,r,0] (Subscript[Y, i,m])(θ,λ)*)

prod20[n_,m_,r_,s_]:=Sum[ a20R[i,n,m,r,s]*FF[i,m],{i,-m,n+r}]

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

(*If m=0 ∧ s=0*)

(*prod00[n_,m_,r_,s_]:= (Subscript[Y, n,0])(θ,λ) (Subscript[Y, r,0])(θ,λ)= (∑_(i=|n-r|)^(n+r))a00R[i,n,0,r,0] (Subscript[Y, i,0])(θ,λ)*)

prod00[n_,m_,r_,s_]:=Sum[ a00R[i,n,m,r,s]*FF[i,m],{i,Abs[n-r],n+r}]

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

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

prodArmFornerRc[n_Integer,m_Integer,r_Integer,s_Integer]:=prodArmFornerRc[n,m,r,s]=Which[n==0&&m==0&&r==0&&s==0,(1/(4π)),n==0&&m==0,(1/(2(Sqrt[π])))FF[r,s],r==0&&s==0,(1/(2(Sqrt[π])))FF[n,m],m==0&&s==0,prod00[n,m,r,s],m<0&&s==0,prod20[n,m,r,s],m==0&&s<0,prod02[n,m,r,s],m>0&&s==0,prod10[n,m,r,s],m==0&&s>0,prod01[n,m,r,s],m>0&&s<0,prod12[n,m,r,s],m<0&&s>0,prod21[n,m,r,s],m<0&&s<0,prod22[n,m,r,s],m>0&&s>0,prod11[n,m,r,s],n<0||r<0,Print["Neither n nor r can be negative numbers"],n<Abs[m]||r<Abs[s],Print["Neither n nor r can be smaller than |m| and |s|, respectively"]]

Clear[prodArmFornerRc]

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

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

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

(*ANOTHER WAY TO PROGRAM THE CODE:*)

(*GENERALIZATION OF THE DEVELOPMENT,IN SERIES OF SPHERICAL HARMONICS IN REAL FORM,OF THE PRODUCT OF SPHERICAL HARMONICS IN REAL FORM*)

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

aR[j_,n_,m_,r_,s_]:=Module[{z},z=Abs[m+s];intFornerR[j,z,n,m,r,s]]

bR[j_,n_,m_,r_,s_]:=Module[{z},z=Abs[m-s];intFornerR[j,z,n,m,r,s]]

aR1[j_,n_,m_,r_,s_]:=Module[{z},z=-m+s;intFornerR[j,z,n,m,r,s]]

bR1[j_,n_,m_,r_,s_]:=Module[{z},z=-Abs[m+s];intFornerR[j,z,n,m,r,s]]

aR2[j_,n_,m_,r_,s_]:=Module[{z},z=m-s;intFornerR[j,z,n,m,r,s]]

prodR[n_,m_,r_,s_]:=
Which[(m>0&&s>0)||(m<0&&s<0),Sum[aR[j,n,m,r,s]*FF[j,Abs[m+s]],{j,Abs[m+s],n+r}]+Sum[bR[j,n,m,r,s]*FF[j,Abs[m-s]],{j,Abs[m-s],n+r}],m>0&&s<0,Sum[aR1[j,n,m,r,s]*FF[j,-m+s],{j,m-s,n+r}]+Sum[bR1[j,n,m,r,s]*FF[j,-Abs[m+s]],{j,Abs[m+s],n+r}],m<0&&s>0,Sum[aR2[j,n,m,r,s]*FF[j,m-s],{j,-m+s,n+r}]+Sum[bR1[j,n,m,r,s]*FF[j,-Abs[m+s]],{j,Abs[m+s],n+r}],m==0&&s>0,Sum[aR[j,n,0,r,s]*FF[j,s],{j,s,n+r}],m==0&&s<0,Sum[aR1[j,n,0,r,s]*FF[j,s],{j,-s,n+r}],m>0&&s==0,Sum[aR[j,n,m,r,0]*FF[j,m],{j,m,n+r}],m<0&&s==0,Sum[aR2[j,n,m,r,0]*FF[j,m],{j,-m,n+r}],m==0&&s==0,Sum[aR[j,n,0,r,0]*FF[j,0],{j,Abs[n-r],n+r}],n==0&&m==0,((FF[r,s])/(2(Sqrt[π]))),r==0&&s==0,((FF[n,m])/(2(Sqrt[π]))),n==0&&m==0&&r==0&&s==0,(1/(4π))]

prodArmFornerRb[n_Integer,m_Integer,r_Integer,s_Integer]:=prodArmFornerRb[n,m,r,s]=Which[(n≥Abs[m])&&(r≥Abs[s]),prodR[n,m,r,s],n<0||r<0,Print["Neither n nor r can be negative numbers"],n<Abs[m]||r<Abs[s],Print["Neither n nor r can be smaller than |m| and |s|, respectively"]]

Clear[prodArmFornerRb]

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

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

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

(*DEFINITIVE FORM OF THE CODE*)

prodDefR[n_,m_,r_,s_]:=Which[(m>0&&s>0)||(m<0&&s<0),Sum[intFornerR[i,Abs[m+s],n,m,r,s]FF[i,Abs[m+s]],{i,Abs[m+s],n+r}]+Sum[intFornerR[j,Abs[m-s],n,m,r,s]FF[j,Abs[m-s]],{j,Abs[m-s],n+r}],(m>0&&s<0)||(m<0&&s>0),Sum[intFornerR[i,-(Abs[m]+Abs[s]),n,m,r,s]FF[i,-(Abs[m]+Abs[s])],{i,Abs[m]+Abs[s],n+r}]+Sum[intFornerR[j,-Abs[m+s],n,m,r,s]FF[j,-Abs[m+s]],{j,Abs[m+s],n+r}],(m==0&&s>0)||(m==0&&s<0),Sum[intFornerR[i,s,n,0,r,s]FF[i,s],{i,Abs[s],n+r}],(m>0&&s==0)||(m<0&&s==0),Sum[intFornerR[i,m,n,m,r,0]FF[i,m],{i,Abs[m],n+r}],m==0&&s==0,Sum[intFornerR[i,0,n,0,r,0]FF[i,0],{i,Abs[n-r],n+r}],n==0&&m==0,((FF[r,s])/(2(Sqrt[π]))),r==0&&s==0,((FF[n,m])/(2(Sqrt[π]))),n==0&&m==0&&r==0&&s==0,(1/(4π))]

prodArmFDefR[n_Integer,m_Integer,r_Integer,s_Integer]:=prodArmFDefR[n,m,r,s]=Which[(n≥Abs[m])&&(r≥Abs[s]),prodDefR[n,m,r,s],n<0||r<0,Print["Neither n nor r can be negative numbers"],n<Abs[m]||r<Abs[s],Print["Neither n nor r can be smaller than |m| and |s|, respectively"]]

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