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

(* CALCULATION AND DEVELOPMENT IN SPHERICAL HARMONICS SERIES IN REAL FORM OF THE FOLLOWING INTEGRALS: *)

(* 1) (∫_(λ1=0)^(2π))(∫_(θ1=-(π/2))^(π/2))(Subscript[P, k])(Cos[θ1]Cos[λ1]) (Subscript[P, r])(Sin[θ1]) (Subscript[P, n])(Cos[γ])Cos[θ1]\[DifferentialD]θ1\[DifferentialD]λ1*)

(* 2) (∫_(λ1=0)^(2π))(∫_(θ1=-(π/2))^(π/2))(Subscript[P, k])(Cos[θ1]Cos[λ1]) (Subscript[Y, r,s])(θ1,λ1) (Subscript[P, n])(Cos[γ])Cos[θ1]\[DifferentialD]θ1\[DifferentialD]λ1*)

(* 3) (∫_(λ1=0)^(2π))(∫_(θ1=-(π/2))^(π/2))(Subscript[P, k])(Cos[θ1]Cos[λ1]) (Subscript[P, r,s])(Sin[θ1]) (Subscript[P, n])(Cos[γ])Cos[θ1]\[DifferentialD]θ1\[DifferentialD]λ1*)

(* (Overscript[OM, ->])=(r Cos[θ] Cos[λ], r Cos[θ] Sin[λ], r Sin[θ]); (Overscript[OM, ->])1=(r1 Cos[θ1] Cos[λ1], r1 Cos[θ1] Sin[λ1], r1 Sin[θ1]) *)

(* Cos[γ]=(((Overscript[OM, ->]).(Overscript[OM, ->])1)/(||(Overscript[OM, ->])||.||(Overscript[OM, ->])1||))=Sin[θ]Sin[θ1]+Cos[θ]Cos[θ1]Cos[λ-λ1]*)

(* If λ=0, Cos[θ]=1 -> Sin[θ]=0, θ=0 -> Cos[γ]=Cos[θ] Cos[λ]*)

(* If λ1=0, Cos[θ1]=1 -> Sin[θ1]=0, θ1=0 -> Cos[γ]=Cos[θ1] Cos[λ1]*)

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

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

γ=Sin[θ]Sin[θ1]+Cos[θ]Cos[θ1]Cos[λ-λ1];

γ0=Cos[θ] Cos[λ];

γ10=Cos[θ1] Cos[λ1];

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

(*Definitions of associated Legendre functions and spherical harmonics in real form.*)

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

YF[n_,m_,θ_,λ_]:=Which[Abs[m]≤n,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] λ]],Abs[m]>0,0]

(* The addition theorem for spherical harmonics in real form *)

pnγ[n_]:=((4π)/(2n+1))Sum[YF[n,m,θ,λ]YF[n,m,θ1,λ1],{m,-n,n}]

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

(* (Subscript[P, n])(0)=(1/(2^n))(∑_(j=0)^n)(((-1))^j)(((((n!)/(j!(n-j)!))))^2), if "n" is even; (Subscript[P, n])(0)=0, if "n" is odd *)

(*  (Subscript[P, n,m])(0)=(((n+m)!n!)/(2^n))(∑_(j=0)^(n-m))((((-1))^j)/(j!(n-m-j)!(m+j)!(n-j)!)), if "n+m" is even and "n≠m"; (Subscript[P, n,n])(0)=(2n-1)!!; (Subscript[P, n,m])(0)=0, if "n+m" is odd *)

polLeF0[n_]:=If[Mod[n,2]==0,(1/(2^n))Sum[(((-1))^j)(((((n!)/(j!(n-j)!))))^2),{j,0,n}],0]

funAsLeF0[n_,m_]:=If[Mod[n+m,2]==0,If[n==m,(2n-1)!!,(((n+m)! n!)/(2^n))Sum[((((-1))^j)/(j!(n-m-j)!(m+j)!(n-j)!)),{j,0,n-m}]],0]

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

(* (Subscript[P, k])(Cos[θ1]Cos[λ1])=(Subscript[P, k])(0) (Subscript[P, k])(Sin[θ1]) + 2(Underoverscript[ ∑, (m=1), k])(((k-m)!)/((k+m)!)) (Subscript[P, k,m])(0) Cos[m λ1] (Subscript[P, k,m])(Sin[θ1])=polFuAsF[k,x]*)

polFuAsF[n_,x_]:=polLeF0[n] LegendreP[n,Sin[θ]]+2Sum[(((n-m)!)/((n+m)!)) funAsLeF0[n,m] funAsLe[n,m,θ] Cos[m λ],{m,1,n}]

(* (Subscript[P, k])(Cos[θ1]Cos[λ1])=(Sqrt[((4π)/(2k+1))])((Subscript[P, k])(0) (Subscript[Y, k,0])(θ1,λ1) + (Sqrt[2]) (Underoverscript[ ∑, (m=1), k])(Sqrt[(((k-m)!)/((k+m)!))]) (Subscript[P, k,m])(0) (Subscript[Y, k,m])(θ1,λ1))=polHaEsF[k,x,θ,λ]*)

polHaEsF[n_,x_,θ_,λ_]:=2(Sqrt[(π/(2n+1))])(polLeF0[n] YF[n,0,θ,λ]+(Sqrt[2])Sum[(Sqrt[(((n-m)!)/((n+m)!))]) funAsLeF0[n,m] YF[n,m,θ,λ],{m,1,n}])

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

(* intFornerR[p,q,n,m,r,s]=(∫_(λ1=0)^(2π))(∫_(θ1=-(π/2))^(π/2)) (Subscript[Y, n,m])(θ1,λ1) (Subscript[Y, r,s])(θ1,λ1) (Subscript[Y, p,q])(θ1,λ1) Cos[θ1] \[DifferentialD]θ1\[DifferentialD]λ1*)

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]

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]

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

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}]]

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

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

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

intFornerR[j_Integer,zz_Integer,n_Integer,m_Integer,r_Integer,s_Integer]:=If[Abs[zz]>j ||Abs[m]>n ||Abs[s ]>r,0,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])}]]

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

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

(* 1) (∫_(λ1=0)^(2π))(∫_(θ1=-(π/2))^(π/2))(Subscript[P, k])(Cos[θ1]Cos[λ1]) (Subscript[P, r])(Sin[θ1]) (Subscript[P, n])(Cos[γ])Cos[θ1]\[DifferentialD]θ1\[DifferentialD]λ1*)

intPoLeKopalR[n_,k_,r_]:=Integrate[Integrate[LegendreP[k,γ10] LegendreP[r,Sin[θ1]] LegendreP[n,γ] Cos[θ1],{θ1,-(π/2),(π/2)}],{λ1,0,2π}]

integralGenPoLeFR[n_,k_,r_]:=((16(π^2))/((2n+1)(Sqrt[(2k+1)(2r+1)])))(polLeF0[k]intFornerR[n,0,k,0,r,0] YF[n,0,θ,λ]+(Sqrt[2])Sum[(Sqrt[(((k-m)!)/((k+m)!))])funAsLeF0[k,m]intFornerR[n,m,k,m,r,0] YF[n,m,θ,λ],{m,1,k}])

intDesGenPoleFR[n_,k_,r_]:=((16(π^2))/((2n+1)(Sqrt[(2k+1)(2r+1)])))(polLeF0[k]intFornerR[n,0,k,0,r,0] yFR[n,0]+(Sqrt[2])Sum[(Sqrt[(((k-m)!)/((k+m)!))])funAsLeF0[k,m]intFornerR[n,m,k,m,r,0] yFR[n,m],{m,1,k}])

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

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

(* 2) (∫_(λ1=0)^(2π))(∫_(θ1=-(π/2))^(π/2))(Subscript[P, k])(Cos[θ1]Cos[λ1]) (Subscript[Y, r,s])(θ1,λ1) (Subscript[P, n])(Cos[γ])Cos[θ1]\[DifferentialD]θ1\[DifferentialD]λ1*)

intHaEsKopalR[n_,k_,r_,s_]:=Integrate[Integrate[LegendreP[k,γ10] YF[r,s,θ1,λ1] LegendreP[n,γ] Cos[θ1],{θ1,-(π/2),(π/2)}],{λ1,0,2π}]

integralGenHarEsFR[n_,k_,r_,s_]:=((8(π^(3/2)))/(2n+1))(Sqrt[(1/(2k+1))])Which[s<0,polLeF0[k] intFornerR[n,s,k,0,r,s] YF[n,s,θ,λ]+(Sqrt[2])Sum[(Sqrt[(((k-m)!)/((k+m)!))])funAsLeF0[k,m](intFornerR[n,-m+s,k,m,r,s] YF[n,-m+s,θ,λ]+intFornerR[n,-Abs[m+s],k,m,r,s] YF[n,-Abs[m+s],θ,λ]),{m,1,k}],s==0,polLeF0[k] intFornerR[n,0,k,0,r,0] YF[n,0,θ,λ]+(Sqrt[2])Sum[(Sqrt[(((k-m)!)/((k+m)!))])funAsLeF0[k,m]intFornerR[n,m,k,m,r,0] YF[n,m,θ,λ],{m,1,k}],s>0,polLeF0[k] intFornerR[n,s,k,0,r,s] YF[n,s,θ,λ]+(Sqrt[2])Sum[(Sqrt[(((k-m)!)/((k+m)!))])funAsLeF0[k,m](intFornerR[n,m+s,k,m,r,s] YF[n,m+s,θ,λ]+intFornerR[n,Abs[m-s],k,m,r,s] YF[n,Abs[m-s],θ,λ]),{m,1,k}]]

intDesGenHarEsFR[n_,k_,r_,s_]:=((8(π^(3/2)))/(2n+1))(Sqrt[(1/(2k+1))])Which[s<0,polLeF0[k] intFornerR[n,s,k,0,r,s] yFR[n,s]+(Sqrt[2])Sum[(Sqrt[(((k-m)!)/((k+m)!))])funAsLeF0[k,m](intFornerR[n,-m+s,k,m,r,s] yFR[n,-m+s]+intFornerR[n,-Abs[m+s],k,m,r,s] yFR[n,-Abs[m+s]]),{m,1,k}],s==0,polLeF0[k] intFornerR[n,0,k,0,r,0] yFR[n,0]+(Sqrt[2])Sum[(Sqrt[(((k-m)!)/((k+m)!))])funAsLeF0[k,m]intFornerR[n,m,k,m,r,0] yFR[n,m],{m,1,k}],s>0,polLeF0[k] intFornerR[n,s,k,0,r,s] yFR[n,s]+(Sqrt[2])Sum[(Sqrt[(((k-m)!)/((k+m)!))])funAsLeF0[k,m](intFornerR[n,m+s,k,m,r,s] yFR[n,m+s]+intFornerR[n,Abs[m-s],k,m,r,s] yFR[n,Abs[m-s]]),{m,1,k}]]

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

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

(* intprodFunAsLegF[n,m,r,s,p,q]=(∫_(-1)^1)(Subscript[P, n,m])(x) (Subscript[P, r,s])(x) (Subscript[P, p,q])(x) \[DifferentialD]x*)

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

intprodFunAsLegF[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],0]

(* invnormafunAsLe[n,m]=(1/(||(Subscript[P, n,m])(x)((||)^2)))*)

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

(* aF[n,m,r,s,p,q]=(((∫_(-1)^1)(Subscript[P, n,m])(x) (Subscript[P, r,s])(x) (Subscript[P, p,q])(x) \[DifferentialD]x)/(||(Subscript[P, n,m])(x)((||)^2)))*)

aF[n_,m_,r_,s_,p_,q_]:=intprodFunAsLegF[n,m,r,s,p,q]invnormafunAsLe[p,q]

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

(* 3) (∫_(λ1=0)^(2π))(∫_(θ1=-(π/2))^(π/2))(Subscript[P, k])(Cos[θ1]Cos[λ1]) (Subscript[P, r,s])(Sin[θ1]) (Subscript[P, n])(Cos[γ])Cos[θ1]\[DifferentialD]θ1\[DifferentialD]λ1*)

intFuAsKopalR[n_,k_,r_,s_]:=Integrate[Integrate[LegendreP[k,γ10] funAsLe[r,s,θ1] LegendreP[n,γ] Cos[θ1],{θ1,-(π/2),(π/2)}],{λ1,0,2π}]

sposF[n_,m_,s_,i_]:=If[m≤n && m+s≤i,Which[m≠0,2π(Sqrt[((2π(n-m)!)/((2n+1)(n+m)!))])intprodFunAsLegF[i,m+s,n,m,0,0] YF[n,m,θ,λ],m==0 ,4π(Sqrt[(π/(2n+1))])intprodFunAsLegF[i,s,n,0,0,0]YF[n,0,θ,λ]],0]

snegF[n_,m_,s_,i_]:=If[m≤n && Abs[s-m]≤i,Which[m≠0,2π(Sqrt[((2π(n-m)!)/((2n+1)(n+m)!))])intprodFunAsLegF[i,-m+s,n,m,0,0] YF[n,m,θ,λ],m==0,4π(Sqrt[(π/(2n+1))])intprodFunAsLegF[i,s,n,0,0,0]YF[n,0,θ,λ]],0]

mzeroF[j_,s_,n_]:=If[Abs[s]≤j,2π(Sqrt[((4π)/(2n+1))])intprodFunAsLegF[j,s,n,0,0,0] YF[n,0,θ,λ],0]

integralGenFuAsFR[n_,k_,r_,s_]:=Which[s<0,polLeF0[k] Sum[aF[k,0,r,s,j,s] mzeroF[j,s,n],{j,Abs[s],k+r}]+2Sum[(((k-m)!)/((k+m)!)) funAsLeF0[k,m] Sum[aF[k,m,r,s,i,s-m] snegF[n,m,s,i],{i,m-s,k+r}],{m,1,k}],s≥0,polLeF0[k] Sum[aF[k,0,r,s,j,s] mzeroF[j,s,n],{j,Abs[s],k+r}]+2Sum[(((k-m)!)/((k+m)!)) funAsLeF0[k,m] Sum[aF[k,m,r,s,i,m+s] sposF[n,m,s,i],{i,m+s,k+r}],{m,1,k}]]

sposFR[n_,m_,s_,i_]:=If[m≤n && m+s≤i,Which[m≠0,2π(Sqrt[((2π(n-m)!)/((2n+1)(n+m)!))])intprodFunAsLegF[i,m+s,n,m,0,0] yFR[n,m],m==0 ,4π(Sqrt[(π/(2n+1))])intprodFunAsLegF[i,s,n,0,0,0]yFR[n,0]],0]

snegFR[n_,m_,s_,i_]:=If[m≤n && Abs[s-m]≤i,Which[m≠0,2π(Sqrt[((2π(n-m)!)/((2n+1)(n+m)!))])intprodFunAsLegF[i,-m+s,n,m,0,0] yFR[n,m],m==0,4π(Sqrt[(π/(2n+1))])intprodFunAsLegF[i,s,n,0,0,0]yFR[n,0]],0]

mzeroFR[j_,s_,n_]:=If[Abs[s]≤j,2π(Sqrt[((4π)/(2n+1))])intprodFunAsLegF[j,s,n,0,0,0] yFR[n,0],0]

intDesGenFuAsFR[n_,k_,r_,s_]:=Which[s<0,polLeF0[k] Sum[aF[k,0,r,s,j,s] mzeroFR[j,s,n],{j,Abs[s],k+r}]+2Sum[(((k-m)!)/((k+m)!)) funAsLeF0[k,m] Sum[aF[k,m,r,s,i,s-m] snegFR[n,m,s,i],{i,m-s,k+r}],{m,1,k}],s≥0,polLeF0[k] Sum[aF[k,0,r,s,j,s] mzeroFR[j,s,n],{j,Abs[s],k+r}]+2Sum[(((k-m)!)/((k+m)!)) funAsLeF0[k,m] Sum[aF[k,m,r,s,i,m+s] sposFR[n,m,s,i],{i,m+s,k+r}],{m,1,k}]]

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

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