ClearAll["Global`*"]

(************************)
Spin = 1;
(**********************
Intermediate definitions
**********************)

Bhel[χ_] := (1/2)*fz2Val /;( Abs[χ] == 2);
Bhel[χ_] := (1/2)*fz1Val /; (Abs[χ] == 1);
Bhel[χ_] := fz0Val /; (Abs[χ] == 0);

chimin = 0;
chimax = 0;
chistep = 0;

If [Spin == 1,chimin = -1];
If [Spin == 1,chimax = 1];
If [Spin == 1,chistep = 1];

If [Spin == 2,chimin = -2];
If [Spin == 2,chimax = 2];
If [Spin == 2,chistep = 1];

FJij[J_,i_,j_] := Sum[Bhel[χ]*FullSimplify[TrigExpand[WignerD[{J,i,χ},-ArcCos[z]]]]*FullSimplify[TrigExpand[WignerD[{J,j,χ},-ArcCos[z]]]],{χ,chimin,chimax,chistep}] /;((J ≥ i)&&(J ≥ j));
FJij[J_,i_,j_] :=0 /;((J < Abs[i])||(J < Abs[j]));
(* ij correspond to i = m-n for A_mn and j = m'-n' for A_m'n' *)

(************************)

If [Spin == 0,fpm = 0];
If [Spin == 0,fmp = 0];
If [Spin == 0,ϕpm = 0];If [Spin == 0,phipm = 0];
If [Spin == 0,ϕmp = 0];If [Spin == 0,phimp = 0];

If [Spin == 0,fp0 = 0];
If [Spin == 0,f0p = 0];
If [Spin == 0,f0m = 0];
If [Spin == 0,fm0 = 0];
If [Spin == 0,ϕp0 = 0];If [Spin == 0,phip0 = 0];
If [Spin == 0,ϕ0p = 0];If [Spin == 0,phi0p = 0];
If [Spin == 0,ϕ0m = 0];If [Spin == 0,phi0m = 0];
If [Spin == 0,ϕm0 = 0];If [Spin == 0,phim0 = 0];

If [Spin == 0,fz0Val = 1];

(*****)

If [Spin == 1,fpm = 0];
If [Spin == 1,fmp = 0];
If [Spin == 1,ϕpm = 0];If [Spin == 1,phipm = 0];
If [Spin == 1,ϕmp = 0];If [Spin == 1,phimp = 0];

If [Spin == 1,fz2Val = 0];

(** spin-independent terms! **)

Array[terms,45];

terms[1] =FullSimplify[ FJij[Spin,0,0]*((1/4) f00 (-1+(x^2)) (-1+(y^2)))];
terms[2] = FullSimplify[ FJij[Spin,0,0]*((1/16) fmm (-1+2 R1Val x-(x^2)) (-1+2 R2Val y-(y^2)))];
terms[3] =  FullSimplify[FJij[Spin,0,0]*((1/16) fpp (1+2 R1Val x+(x^2)) (1+2 R2Val y+(y^2)))];
terms[4] = FullSimplify[FJij[Spin,1,1]*(-(1/8) fp0 (1+2 R1Val x+(x^2)) (-1+(y^2)))];
terms[5] = FullSimplify[FJij[Spin,1,1]*((1/8) f0m (-1+(x^2)) (-1+2 R2Val y-(y^2)))];
terms[6] =FullSimplify[ FJij[Spin,1,1]*(-(1/8) f0p (-1+(x^2)) (1+2 R2Val y+(y^2)))];
terms[7] = FullSimplify[FJij[Spin,1,1]*((1/8) fm0 (-1+2 R1Val x-(x^2)) (-1+(y^2)))];
terms[8] = FullSimplify[FJij[Spin,2,2]*(-(1/16) fpm (1+2 R1Val x+(x^2)) (-1+2 R2Val y-(y^2)))];
terms[9] = FullSimplify[FJij[Spin,2,2]*(-(1/16) fmp (-1+2 R1Val x-(x^2)) (1+2 R2Val y+(y^2)))];
Print["through 9..."]
terms[10] = FullSimplify[FJij[Spin,0,0]*((1/4) (Sqrt[f00]) (Sqrt[fmm]) (-R1Val+x) (Sqrt[1-(x^2)]) (-R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ-ϕmm])];
terms[11] = FullSimplify[ FJij[Spin,0,0]*((1/4) (Sqrt[f00]) (Sqrt[fpp]) (R1Val+x) (Sqrt[1-(x^2)]) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+ϕpp])];
terms[12] = FullSimplify[FJij[Spin,0,1]*(-(((Sqrt[f00]) (Sqrt[fp0]) (R1Val+x) (Sqrt[1-(x^2)]) (-1+(y^2)) Cos[Φ1-ϕp0])/(2 (Sqrt[2]))))];
terms[13] = FullSimplify[FJij[Spin,0,1]*((((Sqrt[f00]) (Sqrt[f0m]) (-1+(x^2)) (-R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ-ϕ0m+Φ1])/(2 (Sqrt[2]))))];
terms[14] =  FullSimplify[FJij[Spin,0,-1]*(-(((Sqrt[f00]) (Sqrt[f0p]) (-1+(x^2)) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+ϕ0p+Φ1])/(2 (Sqrt[2]))))];
terms[15] =  FullSimplify[FJij[Spin,0,-1]*((((Sqrt[f00]) (Sqrt[fm0]) (-R1Val+x) (Sqrt[1-(x^2)]) (-1+(y^2)) Cos[Φ1+ϕm0])/(2 (Sqrt[2]))))];
terms[16] = FullSimplify[FJij[Spin,0,2]*((1/4) (Sqrt[f00]) (Sqrt[fpm]) (R1Val+x) (Sqrt[1-(x^2)]) (R2Val-y) (Sqrt[1-(y^2)]) Cos[Φ+2 Φ1-ϕpm])];
terms[17] = FullSimplify[ FJij[Spin,0,-2]*((1/4) (Sqrt[f00]) (Sqrt[fmp]) (R1Val-x) (Sqrt[1-(x^2)]) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+2 Φ1+ϕmp])]; (* should be equal to FJij[2,0,-2] *)
Print["through 17..."]
terms[18] =  FullSimplify[FJij[Spin,0,0]*((1/8) (Sqrt[fmm]) (Sqrt[fpp]) (-1+(x^2)) (-1+(y^2)) Cos[2 Φ-ϕmm+ϕpp])];
terms[19] = FullSimplify[FJij[Spin,0,1]*((((Sqrt[fmm]) (Sqrt[fp0]) (-1+(x^2)) (-R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ-Φ1-ϕmm+ϕp0])/(4 (Sqrt[2]))))];
terms[20] =  FullSimplify[FJij[Spin,0,1]*((((Sqrt[f0m]) (Sqrt[fmm]) (-R1Val+x) (Sqrt[1-(x^2)]) (-1+2 R2Val y-(y^2)) Cos[ϕ0m-Φ1-ϕmm])/(4 (Sqrt[2]))))];
terms[21] =  FullSimplify[FJij[Spin,0,-1]*((((Sqrt[f0p]) (Sqrt[fmm]) (-R1Val+x) (Sqrt[1-(x^2)]) (-1+(y^2)) Cos[2 Φ+ϕ0p+Φ1-ϕmm])/(4 (Sqrt[2]))))];
terms[22] = FullSimplify[ FJij[Spin,0,-1]*((((Sqrt[fm0]) (Sqrt[fmm]) (-1+2 R1Val x-(x^2)) (-R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+Φ1+ϕm0-ϕmm])/(4 (Sqrt[2]))))];
terms[23] =  FullSimplify[FJij[Spin,0,2]*((1/8) (Sqrt[fmm]) (Sqrt[fpm]) (-1+(x^2)) (-1+2 R2Val y-(y^2)) Cos[2 Φ1+ϕmm-ϕpm])];
terms[24] =  FullSimplify[FJij[Spin,0,-2]*((1/8) (Sqrt[fmm]) (Sqrt[fmp]) (-1+2 R1Val x-(x^2)) (-1+(y^2)) Cos[2 Φ+2 Φ1-ϕmm+ϕmp])];
Print["through 24..."]
terms[25] =  FullSimplify[FJij[Spin,0,1]*((((Sqrt[fp0]) (Sqrt[fpp]) (1+2 R1Val x+(x^2)) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+Φ1-ϕp0+ϕpp])/(4 (Sqrt[2]))))];
terms[26] =  FullSimplify[FJij[Spin,0,1]*(-(((Sqrt[f0m]) (Sqrt[fpp]) (R1Val+x) (Sqrt[1-(x^2)]) (-1+(y^2)) Cos[2 Φ-ϕ0m+Φ1+ϕpp])/(4 (Sqrt[2]))))];
terms[27] = FullSimplify[ FJij[Spin,0,-1]*((((Sqrt[f0p]) (Sqrt[fpp]) (R1Val+x) (Sqrt[1-(x^2)]) (1+2 R2Val y+(y^2)) Cos[ϕ0p+Φ1-ϕpp])/(4 (Sqrt[2]))))];
terms[28] =  FullSimplify[FJij[Spin,0,-1]*(-(((Sqrt[fm0]) (Sqrt[fpp]) (-1+(x^2)) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ-Φ1-ϕm0+ϕpp])/(4 (Sqrt[2]))))];
terms[29] =  FullSimplify[FJij[Spin,0,2]*(-(1/8) (Sqrt[fpm]) (Sqrt[fpp]) (1+2 R1Val x+(x^2)) (-1+(y^2)) Cos[2 Φ+2 Φ1-ϕpm+ϕpp])];
terms[30] =  FullSimplify[FJij[Spin,0,-2]*(-(1/8) (Sqrt[fmp]) (Sqrt[fpp]) (-1+(x^2)) (1+2 R2Val y+(y^2)) Cos[2 Φ1+ϕmp-ϕpp])];  (* should be equal to FJij[2,0,-2] *)
Print["through 30..."]
terms[31] =  FullSimplify[FJij[Spin,1,1]*((1/4) (Sqrt[f0m]) (Sqrt[fp0]) (R1Val+x) (Sqrt[1-(x^2)]) (R2Val-y) (Sqrt[1-(y^2)]) Cos[Φ-ϕ0m+ϕp0])];
terms[32] =  FullSimplify[FJij[Spin,1,-1]*((1/4) (Sqrt[f0p]) (Sqrt[fp0]) (R1Val+x) (Sqrt[1-(x^2)]) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+ϕ0p+2 Φ1-ϕp0])];
terms[33] = FullSimplify[FJij[Spin,1,-1]*((1/4) (Sqrt[fm0]) (Sqrt[fp0]) (-1+(x^2)) (-1+(y^2)) Cos[2 Φ1+ϕm0-ϕp0])];
terms[34] =  FullSimplify[FJij[Spin,1,2]*((((Sqrt[fp0]) (Sqrt[fpm]) (1+2 R1Val x+(x^2)) (R2Val-y) (Sqrt[1-(y^2)]) Cos[Φ+Φ1+ϕp0-ϕpm])/(4 (Sqrt[2]))))];
terms[35] = FullSimplify[FJij[Spin,1,-2]*(-(((Sqrt[fmp]) (Sqrt[fp0]) (-1+(x^2)) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+3 Φ1+ϕmp-ϕp0])/(4 (Sqrt[2]))))];
Print["through 35..."]
terms[36] =  FullSimplify[FJij[Spin,1,-1]*((1/4) (Sqrt[f0m]) (Sqrt[f0p]) (-1+(x^2)) (-1+(y^2)) Cos[2 Φ-ϕ0m+ϕ0p+2 Φ1])];
terms[37] =  FullSimplify[FJij[Spin,1,-1]*((1/4) (Sqrt[f0m]) (Sqrt[fm0]) (-R1Val+x) (Sqrt[1-(x^2)]) (-R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ-ϕ0m+2 Φ1+ϕm0])];
terms[38] =  FullSimplify[FJij[Spin,1,2]*(-(((Sqrt[f0m]) (Sqrt[fpm]) (R1Val+x) (Sqrt[1-(x^2)]) (-1+2 R2Val y-(y^2)) Cos[ϕ0m+Φ1-ϕpm])/(4 (Sqrt[2]))))];
terms[39] = FullSimplify[ FJij[Spin,1,-2]*((((Sqrt[f0m]) (Sqrt[fmp]) (-R1Val+x) (Sqrt[1-(x^2)]) (-1+(y^2)) Cos[2 Φ-ϕ0m+3 Φ1+ϕmp])/(4 (Sqrt[2]))))];
Print["through 39..."]
terms[40] =  FullSimplify[FJij[Spin,-1,-1]*((1/4) (Sqrt[f0p]) (Sqrt[fm0]) (R1Val-x) (Sqrt[1-(x^2)]) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+ϕ0p-ϕm0])];
terms[41] = FullSimplify[ FJij[Spin,-1,2]*(-(((Sqrt[f0p]) (Sqrt[fpm]) (R1Val+x) (Sqrt[1-(x^2)]) (-1+(y^2)) Cos[2 Φ+ϕ0p+3 Φ1-ϕpm])/(4 (Sqrt[2]))))];
terms[42] = FullSimplify[ FJij[Spin,-1,-2]*((((Sqrt[f0p]) (Sqrt[fmp]) (R1Val-x) (Sqrt[1-(x^2)]) (1+2 R2Val y+(y^2)) Cos[ϕ0p-Φ1-ϕmp])/(4 (Sqrt[2]))))];
Print["through 42..."]
terms[43] = FullSimplify[FJij[Spin,-1,2]*((((Sqrt[fm0]) (Sqrt[fpm]) (-1+(x^2)) (-R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+3 Φ1+ϕm0-ϕpm])/(4 (Sqrt[2]))))];
terms[44] = FullSimplify[ FJij[Spin,-1,-2]*(-(((Sqrt[fm0]) (Sqrt[fmp]) (-1+2 R1Val x-(x^2)) (R2Val+y) (Sqrt[1-(y^2)]) Cos[Φ+Φ1-ϕm0+ϕmp])/(4 (Sqrt[2]))))];
Print["through 44..."]
terms[45] = FullSimplify[ FJij[Spin,2,-2]*((1/8) (Sqrt[fmp]) (Sqrt[fpm]) (-1+(x^2)) (-1+(y^2)) Cos[2 Φ+4 Φ1+ϕmp-ϕpm])];

fout = "test.txt";
foutIntAll = "test_IntAll.txt";
foutIntNoH1 = "test_IntNoH1.txt";
foutIntNoH2 = "test_IntNoH2.txt";
foutIntNoHs = "test_IntNoHs.txt";
foutIntNoPhi = "test_IntNoPhi.txt";
foutIntNoPhi1 = "test_IntNoPhi1.txt";

For[i=1,i≤45,i=i+1,
	  If[Mod[i,10] == 0, Print[i,"..."];]

WriteString[fout,"value += "];
PutAppend[CForm[FullSimplify[
	terms[i]]]/.{ Φ->Phi,Φ1->Phi1,x->h1,y->h2,z->hs,ϕmm->phimm,ϕpp->phipp,ϕp0->phip0,ϕ0m->phi0m,ϕ0p->phi0p,ϕm0->phim0,ϕpm->phipm,ϕmp->phimp},fout];
WriteString[fout,"; \n"];

]

(*********)
(* Integration Limits *)
(*********)
xlo = -1;
xhi = 1;
ylo = -1;
yhi = 1;
zlo = -1;
zhi = 1;

PhiLo = -Pi;
PhiHi = Pi;
Phi1Lo = -Pi;
Phi1Hi = Pi;

(*********)

For[i=1,i≤45,i=i+1,
	  If[Mod[i,10] == 0, Print[i,"..."];]
	If[terms[i]==0,PutAppend[0,foutIntAll],Print["False!"],
WriteString[foutIntAll,"value += "];
PutAppend[CForm[FullSimplify[
	Integrate[terms[i],{x,xlo,xhi},{y,ylo,yhi},{z,zlo,zhi},{Φ,PhiLo,PhiHi},{Φ1,Phi1Lo,Phi1Hi}]]]/.{ Φ->Phi,Φ1->Phi1,x->h1,y->h2,z->hs,ϕmm->phimm,ϕpp->phipp,ϕp0->phip0,ϕ0m->phi0m,ϕ0p->phi0p,ϕm0->phim0,ϕpm->phipm,ϕmp->phimp},foutIntAll];
   WriteString[foutIntAll,"; \n"];
]
]

For[i=1,i≤45,i=i+1,
	  If[Mod[i,10] == 0, Print[i,"..."];]
	If[terms[i]==0,PutAppend[0,foutIntNoH1],Print["False!"],
WriteString[foutIntNoH1,"value += "];
PutAppend[CForm[FullSimplify[
	Integrate[terms[i],{y,ylo,yhi},{z,zlo,zhi},{Φ,PhiLo,PhiHi},{Φ1,Phi1Lo,Phi1Hi}]]]/.{ Φ->Phi,Φ1->Phi1,x->h1,y->h2,z->hs,ϕmm->phimm,ϕpp->phipp,ϕp0->phip0,ϕ0m->phi0m,ϕ0p->phi0p,ϕm0->phim0,ϕpm->phipm,ϕmp->phimp},foutIntNoH1];
WriteString[foutIntNoH1,"; \n"];
]
]

For[i=1,i≤45,i=i+1,
	  If[Mod[i,10] == 0, Print[i,"..."];]
	If[terms[i]==0,PutAppend[0,foutIntNoH2],Print["False!"],
WriteString[foutIntNoH2,"value += "];
PutAppend[CForm[FullSimplify[
	Integrate[terms[i],{x,xlo,xhi},{z,zlo,zhi},{Φ,PhiLo,PhiHi},{Φ1,Phi1Lo,Phi1Hi}]]]/.{ Φ->Phi,Φ1->Phi1,x->h1,y->h2,y->hs,ϕmm->phimm,ϕpp->phipp,ϕp0->phip0,ϕ0m->phi0m,ϕ0p->phi0p,ϕm0->phim0,ϕpm->phipm,ϕmp->phimp},foutIntNoH2];
WriteString[foutIntNoH2,"; \n"];
]
]

For[i=1,i≤45,i=i+1,
	  If[Mod[i,10] == 0, Print[i,"..."];]
	If[terms[i]==0,PutAppend[0,foutIntNoHs],Print["False!"],
WriteString[foutIntNoHs,"value += "];
PutAppend[CForm[FullSimplify[
	Integrate[terms[i],{x,xlo,xhi},{y,ylo,yhi},{Φ,PhiLo,PhiHi},{Φ1,Phi1Lo,Phi1Hi}]]]/.{ Φ->Phi,Φ1->Phi1,x->h1,y->h2,z->hs,ϕmm->phimm,ϕpp->phipp,ϕp0->phip0,ϕ0m->phi0m,ϕ0p->phi0p,ϕm0->phim0,ϕmm->phipm,ϕmp->phimp},foutIntNoHs];
WriteString[foutIntNoHs,"; \n"];
]
]

For[i=1,i≤45,i=i+1,
	  If[Mod[i,10] == 0, Print[i,"..."];]
	If[terms[i]==0,PutAppend[0,foutIntNoPhi],Print["False!"],
WriteString[foutIntNoPhi,"value += "];
PutAppend[CForm[FullSimplify[
	Integrate[terms[i],{x,-1,1},{y,-1,1},{z,-1,1},{Φ1,Phi1Lo,Phi1Hi} ]]]/.{ Φ->Phi,Φ1->Phi1,x->h1,y->h2,z->hs,ϕmm->phimm,ϕpp->phipp,ϕp0->phip0,ϕ0m->phi0m,ϕ0p->phi0p,ϕm0->phim0,ϕpm->phipm,ϕmp->phimp},foutIntNoPhi];
WriteString[foutIntNoPhi,"; \n"];
]
]

For[i=1,i≤45,i=i+1,
	  If[Mod[i,10] == 0, Print[i,"..."];]
	If[terms[i]==0,PutAppend[0,foutIntNoPhi1],Print["False!"],
WriteString[foutIntNoPhi1,"value += "];
PutAppend[CForm[FullSimplify[
	Integrate[terms[i],{x,xlo,xhi},{y,ylo,yhi},{z,zlo,zhi},{Φ,PhiLo,PhiHi}]]]/.{ Φ->Phi,Φ1->Phi1,x->h1,y->h2,z->hs,ϕmm->phimm,ϕpp->phipp,ϕp0->phip0,ϕ0m->phi0m,ϕ0p->phi0p,ϕm0->phim0,ϕpm->phipm,ϕmp->phimp},foutIntNoPhi1];
WriteString[foutIntNoPhi1,"; \n"];
]
]