NEWTONIAN FLUX COMPUTATION

$Assumptions=ν∈Reals&&0<ν<1/4&&ϕ∈Reals

PRELIMINARIES

ϵ[el_,em_]=Mod[el+em,2];
n0[el_,em_]=(((I em))^el)((8π)/((2el+1)!!))(Sqrt[(((el+1)(el+2))/(el(el-1)))]);
n1[el_,em_]=-(((I em))^el)((16 I π)/((2el+1)!!))(Sqrt[(((2el+1)(el+2)((el^2)-(em^2)))/(el(el-1)(el+1)(2el-1)))]);
n[el_,em_]=If[Mod[el+em,2]==0,n0[el,em],n1[el,em]];
Y[el_,em_]=SphericalHarmonicY[el-ϵ[el,em],-em,(π/2),ϕ];

(Subscript[c, ℓ+ϵ])(ν) in terms of (Subscript[m, 1]), (Subscript[m, 2])

M=m1+m2;
X1=m1/M;
X2=m2/M;

c[el_,em_]=(X2^(el+ϵ[el,em]-1))+(((-1))^em)(X1^(el+ϵ[el,em]-1));
Simplify[Table[c[l,m],{l,2,8},{m,1,l}]]

(Subscript[c, ℓ+ϵ])(ν) in terms of ν

Clear[x1,x2]

sys=Flatten[Solve[x1+x2==1&&x1 x2==ν,{x1,x2}]]

x1=sys[[1,2]];
x2=sys[[2,2]];
x1Alt=sys[[3,2]];
x2Alt=sys[[4,2]];
X1=x1
X2=x2

c[el_,em_]=(X2^(el+ϵ[el,em]-1))+(((-1))^em)(X1^(el+ϵ[el,em]-1));
Simplify[Table[c[l,m],{l,2,8},{m,1,l}]]

CHECK FlmNewt i.e. Newtonian Flux

DEFINE Flm using Eqs.(2) and (4) of 0811.2069
NOTE: the paper seems to be missing a factor of c^3, see F22

FNewt[el_,em_]=With[{Ω=(((c^3)(x^(3/2)))/(G M))},Simplify[(2/(16 π G ))(em^2)(Ω^2)(((G (M/(c^2))ν n[el,em]c[el,em](x^((el+ϵ[el,em])/2))))^2)Y[el,em]Y[el,-em]]];

FNewt[2,2]

FlmNewtMath=Flatten[(G/((ν^2)(c^2)))Simplify[Table[FNewt[l,m],{l,2,8},{m,1,l}]]];

C’s FlmNewt from TEOBResumSFlux.c

CNlm={0.17777777777777778,6.4,0.0007936507936507937,0.5079365079365079,8.678571428571429,2.2675736961451248*10^-05,0.008062484252960444,1.0414285714285714,14.447971781305114,5.010421677088344*10^-08,0.0006384836014465644,0.03106534090909091,1.9614216236438458,25.688197074915823,8.898517797026696*10^-10,4.464920289836115*10^-06,0.003830520129221428,0.08778390483440988,3.584607999290539,46.98226573426573,1.0695333890657086*10^-12,2.3656367312710969*10^-07,4.972309783123969*10^-05,0.013643024456211269,0.21542115380351798,6.472046810332524,87.1329124642076,1.2233225038333268*10^-14,9.27700678929842*10^-10,4.468579053461083*10^-06,0.0002675102834551229,0.03813282951314997,0.4894825318738884,11.627213264559293,162.79300083906728};
x5=x*x*x*x*x;
x6=x*x5;
x7=x*x6;
x8=x*x7;
x9=x*x8;
x10=x*x9;
x11=x*x10;
x12=x*x11;
sp2=(1-4ν);
sp4=sp2*(((1-2ν))^2);
sp3=(((1-3ν))^2);
sp5=(((1-5ν+5(ν^2)))^2);
sp6=sp2*(((3(ν^2)-4ν+1))^2);
sp7=(((1-7ν+14(ν^2)-7(ν^3)))^2);
sp8=sp2*(((1-6ν+10(ν^2)-4(ν^3)))^2);
spx={sp2*x6,x5,sp2*x6,sp3*x7,sp2*x6,sp4*x8,sp3*x7,sp4*x8,sp3*x7,sp4*x8,sp5*x9,sp4*x8,sp5*x9,sp4*x8,sp6*x10,sp5*x9,sp6*x10,sp5*x9,sp6*x10,sp5*x9,sp6*x10,sp7*x11,sp6*x10,sp7*x11,sp6*x10,sp7*x11,sp6*x10,sp8*x12,sp7*x11,sp8*x12,sp7*x11,sp8*x12,sp7*x11,sp8*x12,(((7*(ν^3)-14*(ν^2)+7*ν-1))^2)*x11};
FlmNewtC=CNlm*spx;

C++’s FlmNewt from TEOBResumSFlux.cpp

Nlm0=8./45*sp2*x6;
Nlm1=32./5.*x5;
Nlm2=1./1260*sp2*x6;
Nlm3=32./63*sp3*x7;
Nlm4=243./28*sp2*x6;
Nlm5=1./44100*sp4*x8;
Nlm6=32./3969*sp3*x7;
Nlm7=729./700*sp4*x8;
Nlm8=8192./567*sp3*x7;
Nlm9=1./19958400*sp4*x8;
Nlm10=256./400950*sp5*x9;
Nlm11=2187./70400*sp4*x8;
Nlm12=131072./66825*sp5*x9;
Nlm13=1953125./76032*sp4*x8;
Nlm14=1./1123782660*sp6*x10;
Nlm15=128./28667925*sp5*x9;
Nlm16=59049./15415400*sp6*x10;
Nlm17=4194304./47779875*sp5*x9;
Nlm18=48828125./13621608*sp6*x10;
Nlm19=839808./17875*sp5*x9;
Nlm20=1./(9.3498717312*10^11)*sp6*x10;
Nlm21=32./(1.35270135*10^8)*sp7*x11;
Nlm22=1594323./(3.2064032*10^10)*sp6*x10;
Nlm23=4194304./(3.07432125*10^8)*sp7*x11;
Nlm24=1220703125./(5.666588928*10^9)*sp6*x10;
Nlm25=5668704./875875*sp7*x11;
Nlm26=96889010407./1111968000*sp6*x10;
Nlm27=1./(8.174459284992*10^13)*sp8*x12;
Nlm28=32./(3.4493884425*10^10)*sp7*x11;
Nlm29=177147./(3.96428032*10^10)*sp8*x12;
Nlm30=4194304./(1.5679038375*10^10)*sp7*x11;
Nlm31=30517578125./(8.00296713216*10^11)*sp8*x12;
Nlm32=51018336./(1.04229125*10^8)*sp7*x11;
Nlm33=4747561509943./(4.083146496*10^11)*sp8*x12;
Nlm34=274877906944./1688511825.*sp7*x11;

FlmNewtCpp={Nlm0,Nlm1,Nlm2,Nlm3,Nlm4,Nlm5,Nlm6,Nlm7,Nlm8,Nlm9,
Nlm10,Nlm11,Nlm12,Nlm13,Nlm14,Nlm15,Nlm16,Nlm17,Nlm18,Nlm19,
Nlm20,Nlm21,Nlm22,Nlm23,Nlm24,Nlm25,Nlm26,Nlm27,Nlm28,Nlm29,
Nlm30,Nlm31,Nlm32,Nlm33,Nlm34};

COMPARE Mathematica and C/C++

(* FIRST: check CNlm's, if they agree to Machine Epsilon, replace them with Mathematica's exact expressions *)
CNlmMath=FullSimplify[(FlmNewtMath/spx)];
Norm[CNlm-CNlmMath]
CNlmCpp=FullSimplify[(FlmNewtCpp/spx)];
Norm[CNlmCpp-CNlmMath]

(* COMPARE EXACTLY *)

DiffCMath=FullSimplify[CNlmMath*spx-FlmNewtMath]
DiffCppMath=FullSimplify[CNlmMath*spx-FlmNewtMath]

(* COMPARE Numerically *)

DiffCMath=Norm[FullSimplify[FlmNewtCpp-FlmNewtMath]/.x->1/.ν->1/10]
DiffCppMath=Norm[FullSimplify[FlmNewtC-FlmNewtMath]/.x->1/.ν->1/10]
DiffCppC=Norm[FullSimplify[FlmNewtC-FlmNewtCpp]/.x->1/.ν->1/10]

EXCEPTIONAL MODES: k=0,2,4,5,79,11,13 
COMPARE C and C++ for ALL 4 CASES

cases={"BBH_NOspin", "BNS_NOspin", "BBH_spin", "BNS_spin"};
Clear[i];
Do[
case=cases[[i]];
(* FROM Flux.cpp:265-320 *)
sp2=If[case=="BBH_NOspin",1-4ν,1];
sp4=If[(case=="BBH_spin")||(case=="BNS_spin"),(((1-2ν))^2),(1-4ν)(((1-2ν))^2)];
Nlm0=8./45*sp2*x6;
Nlm2=1./1260*sp2*x6;
Nlm4=243./28*sp2*x6;
Nlm5=1./44100*sp4*x8;
Nlm7=729./700*sp4*x8;
Nlm9=If[(case=="BBH_spin")||(case=="BNS_spin"),1./19958400*(1-4ν)(((1-2ν))^2)*x8,1./19958400*sp4*x8];
Nlm11=If[(case=="BBH_spin")||(case=="BNS_spin"),2187./70400*x8*(1-4ν)(((1-2ν))^2),2187./70400*sp4*x8];
Nlm13=If[(case=="BBH_spin")||(case=="BNS_spin"),1953125./76032*x8*(1-4ν)(((1-2ν))^2),1953125./76032*sp4*x8];
FlmNewtCpp={Nlm0,Nlm1,Nlm2,Nlm3,Nlm4,Nlm5,Nlm6,Nlm7,Nlm8,Nlm9,
Nlm10,Nlm11,Nlm12,Nlm13,Nlm14,Nlm15,Nlm16,Nlm17,Nlm18,Nlm19,
Nlm20,Nlm21,Nlm22,Nlm23,Nlm24,Nlm25,Nlm26,Nlm27,Nlm28,Nlm29,
Nlm30,Nlm31,Nlm32,Nlm33,Nlm34};
(*From Flux.c:245-265 *)
FlmNewtC[[1]]=If[(case=="BBH_spin")||(case=="BNS_spin")||(case=="BNS_NOspin"),CNlm[[1]]*x6,FlmNewtMath[[1]]];
FlmNewtC[[3]]=If[(case=="BBH_spin")||(case=="BNS_spin")||(case=="BNS_NOspin"),CNlm[[3]]*x6,FlmNewtMath[[3]]];
FlmNewtC[[5]]=If[(case=="BBH_spin")||(case=="BNS_spin")||(case=="BNS_NOspin"),CNlm[[5]]*x6,FlmNewtMath[[5]]];
FlmNewtC[[6]]=If[(case=="BBH_spin")||(case=="BNS_spin"),CNlm[[6]]*(((1-2ν))^2)*x8,FlmNewtMath[[6]]];
FlmNewtC[[8]]=If[(case=="BBH_spin")||(case=="BNS_spin"),CNlm[[8]]*(((1-2ν))^2)*x8,FlmNewtMath[[8]]];
NormAll=Norm[FullSimplify[2((FlmNewtC-FlmNewtCpp)/FlmNewtMath)]/.x->1/2/.ν->1/10];
Print[cases[[i]],"\t\tNorm[Rel Diff C-C++] = ",NormAll];

,{i,1,4}]

(OLDER version pre-loop)

Change i from 1 to 4 and Execute the new 3 sub-sections together

i=4;
case=cases[[i]]

FROM Flux.cpp:265-320

{spx[[10]],spx[[12]],spx[[14]]}/(x^8)
sp2=If[case=="BBH_NOspin",1-4ν,1];
sp4=If[(case=="BBH_spin")||(case=="BNS_spin"),(((1-2ν))^2),(1-4ν)(((1-2ν))^2)];
{sp4,sp4,sp4}

Nlm0=8./45*sp2*x6;
Nlm2=1./1260*sp2*x6;
Nlm4=243./28*sp2*x6;
Nlm5=1./44100*sp4*x8;
Nlm7=729./700*sp4*x8;
Nlm9=1./19958400*sp4*x8
Nlm11=2187./70400*sp4*x8;
Nlm13=1953125./76032*sp4*x8;
Nlm9=If[(case=="BBH_spin")||(case=="BNS_spin"),1./19958400*(1-4ν)sp4*x8,1./19958400*sp4*x8]
Nlm11=If[(case=="BBH_spin")||(case=="BNS_spin"),Nlm11,2187./70400*sp4*x8]
Nlm13=If[(case=="BBH_spin")||(case=="BNS_spin"),Nlm13,1953125./76032*sp4*x8]

FlmNewtCpp={Nlm0,Nlm1,Nlm2,Nlm3,Nlm4,Nlm5,Nlm6,Nlm7,Nlm8,Nlm9,
Nlm10,Nlm11,Nlm12,Nlm13,Nlm14,Nlm15,Nlm16,Nlm17,Nlm18,Nlm19,
Nlm20,Nlm21,Nlm22,Nlm23,Nlm24,Nlm25,Nlm26,Nlm27,Nlm28,Nlm29,
Nlm30,Nlm31,Nlm32,Nlm33,Nlm34};

From Flux.c:245-265

FlmNewtC[[1]]=If[(case=="BBH_spin")||(case=="BNS_spin")||(case=="BNS_NOspin"),CNlm[[1]]*x6,FlmNewtMath[[1]]];
FlmNewtC[[3]]=If[(case=="BBH_spin")||(case=="BNS_spin")||(case=="BNS_NOspin"),CNlm[[3]]*x6,FlmNewtMath[[3]]];
FlmNewtC[[5]]=If[(case=="BBH_spin")||(case=="BNS_spin")||(case=="BNS_NOspin"),CNlm[[5]]*x6,FlmNewtMath[[5]]];
FlmNewtC[[6]]=If[(case=="BBH_spin")||(case=="BNS_spin"),CNlm[[6]]*(((1-2ν))^2)*x8,FlmNewtMath[[6]]];
FlmNewtC[[8]]=If[(case=="BBH_spin")||(case=="BNS_spin"),CNlm[[8]]*(((1-2ν))^2)*x8,FlmNewtMath[[8]]];

NormAll=Norm[FullSimplify[2((FlmNewtC-FlmNewtCpp)/FlmNewtMath)]/.x->1/.ν->1/10];
ProbModesNorm=Take[2((FlmNewtC-FlmNewtCpp)/FlmNewtMath),{10,14,2}]/.x->1/.ν->1/10;
NonPromModesNorm=Norm[Drop[2((FlmNewtC-FlmNewtCpp)/FlmNewtMath),{10,14,2}]/.x->1/.ν->1/10];
If[NormAll≥(10^(-14)),Print["Problematic modes: (5,1), (5,3), (5,5)\nRel. Diff C and C++ = ", ProbModesNorm,
"\nThe other modes are FINE: ",NonPromModesNorm],Print["Norm of Rel Diff between C and C++\n",NormAll]]

CHECK hlmNewt i.e. Newtonian waveform

AMPLITUDE COMPUTATION
Easier since we take Norm of everything so minus signs matter not
Results here FULLY AGREE with Waveform.c:24

Amps=Table[Abs[n[l,m]*Y[l,m]],{l,2,8},{m,1,l}]~N~18/.ϕ->0

MyAmps=Flatten[Amps]

CForm[MyAmps]

PHASES
Most of my phases DISAGREE with C’s by ±π
I divide phase difference by 2π and Mod by 1, if the result is ≠0 then phases disagree
NOTE: C has no phases for Modes (8,1) to (8,7). At the BOTTOM I construct these.

Table[((n[l,m]*Y[l,m])/(Abs[n[l,m]*Y[l,m]Exp[I m ϕ]])),{l,2,4},{m,1,l}]

Phases=Arg[Simplify[Table[((c[l,m]n[l,m]*Y[l,m]Exp[I m ϕ])/(Abs[c[l,m]n[l,m]*Y[l,m]])),{l,2,8},{m,1,l}],ϕ∈Reals&&ν<1/4]]

MyPhases=Flatten[Phases]

MyPhasesforC =MyPhases~N~18

Compare the Phase Difference

Cphases=Rationalize[{3*Pi/2,Pi,3*Pi/2,0.,Pi/2,Pi/2,0.,3*Pi/2,Pi,Pi/2,Pi,3*Pi/2,0.,Pi/2,3*Pi/2,Pi,Pi/2,0,3*Pi/2,Pi,3*Pi/2,0.,Pi/2,Pi,3*Pi/2,0.,Pi/2,0,0,0,0,0,0,0,Pi}];
Abs[MyPhases-Cphases]
Print[Style["If it's not zero, PHASES DISAGREE:\n",20],Mod[((MyPhases-Cphases)/(2π)),1]]

ATTEMPT TO RECONCILE THE PHASE DIFFERENCES
MY THEORY: The Alm[k] that goes into hlmNewt->ampli[k] in C differs from(Subscript[c, ℓ+ϵ])(ν) of Table V of 0811.206 by some (-) signs.
 Therefore there is ±π Phase difference for most modes.
Specifically, there are some differences in the signs of the(Subscript[c, ℓ+ϵ])(ν) and p2, p3, p4,p6,p7, p8 of C. Here I call (Subscript[c, ℓ+ϵ])(ν) by clms.

(* Define the clms *)
clms=Flatten[Simplify[Table[c[l,m],{l,2,8},{m,1,l}]]];

ONLY p5 of C and (Subscript[c, 6])(ν) of Table V. FULLY AGREE therefore there is no phase disagreement for modes containing p5, which are k=11,13,16,18,20
First, I show this phase agreement

p5=1-5*ν+5*(ν^2); (* exactly as above so there should be no phase difference *)
kvals={11,13,16,18,20};
Print["p5 agrees with (Subscript[c, 6])",clms[[kvals]]/p5]
(* Showing that the phases are equal for these modes *)
Print[Style["These 5 phases AGREE:\n",20],MyPhases[[kvals]]-Cphases[[kvals]]]

I think the Phase DISAGREEMENT could be due to sign(m1-m2) in X12 introduced above Table V
If we use m2>m1 then (Subscript[X, 2])-(Subscript[X, 1])=(Sqrt[1-4ν])
Let me introduce a factor (((-1))^r) into (Subscript[p, 2]) which I can set by r=±1

p1=1;
p2=(((-1))^r)(Sqrt[1-4ν]); (* Because X2≥ X1, it is X2-X1 = (Sqrt[1-4ν])  *)
p3=(3*ν-1); (* p3 is off by -1 = Pi phase difference *)
p4=(2*ν-1)*p2; 
p6=(1-4*ν+3*(ν^2))*p2;
p7=7*(ν^3)-14*(ν^2)+7*ν-1; (* This is off by -1 = Pi phase difference*)
p8=(4*(ν^3)-10*(ν^2)+6*ν-1)*p2;
(* This is Alm[35] of Waveform.c:111 *)
Alm={p2,p1,p2,p3,p2,p4,p3,p4,p3,p4,p5,p4,p5,p4,p6,p5,p6,p5,p6,p5,p6,p7,p6,p7,p6,p7,p6,p8,p7,p8,p7,p8,p7,p8,p7};

$Assumptions=ν∈Reals&&0<ν<1/4

ℓ=2 CHECK: (2,2) mode has p1=1 so no problem there, BUT (2,1) has p2 =±(Sqrt[1-4ν]), so PhaseDiff = π

kvals={1,2};
PhaseDiff=MyPhases[[kvals]]-Cphases[[kvals]];
Print["Phase Diff for (2,1),(2,2) = ",PhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->1]]-Simplify[Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["Extra minus sign into p2 FIXES (2,1)\nPhase Diff for (2,1), (2,2) = ",FixedPhaseDiff]

ℓ=3 CHECK: (3,1) and (3,3) have p2’s. (3,2) uses p3=(3ν-1) BUT it should be 1-3ν so we expect all 3 phases to be off by ±π

kvals={3,4,5};
PhaseDiff=MyPhases[[kvals]]-Cphases[[kvals]];
Print["Phase Diff for ℓ=3 modes = ",PhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->1]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["Extra minus sign into p2 FIXES (3,1), (3,3)\nPhase Diff for ℓ=3 modes = ",FixedPhaseDiff]

ℓ=4 CHECK: (4,2), (4,4) have p3’s which are OFF by -1 so that’s ±π .
(4,1), (4,3) have p4’s which are also OFF by -1, BUT the phase of (4,1), (4,3) somehow AGREE mod 2π 
I think this is due to the ± ambiguity of X12 in the C code as I show below.

kvals=Range[6,9];
PhaseDiff=MyPhases[[kvals]]-Cphases[[kvals]];
Print["Phase Diff for ℓ=4 modes = ",PhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->0]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["No minus sign in p2 which goes into p4, (4,1), (4,3) BROKEN\nPhase Diff for ℓ=4 modes = ",FixedPhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->1]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["Extra minus sign into p2->p4 FIXES (4,1), (4,3)\nPhase Diff for ℓ=4 modes = ",FixedPhaseDiff]

ℓ=5 CHECK:  Phases ALL AGREE. (5,2), (5,4) have p5’s so phases agree as shown above.
(5,1), (5,3),(5,5) have p4’s which are OFF by -1, so these phases should be off by π 
BUT they AGREE (mod 2π) due to the ± ambiguity of X12 mentioned above.

kvals=Range[10,14];
PhaseDiff=MyPhases[[kvals]]-Cphases[[kvals]];
Print["Phase Diff for ℓ=5 modes mod 2Pi = ",Mod[(PhaseDiff/(2π)),1]]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->0]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["No minus sign in p2 which goes into p4, (5,1), (5,3), (5,5) BROKEN\nPhase Diff for ℓ=5 modes = ",FixedPhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->1]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["Extra minus sign into p2->p4 FIXES (5,1),(5,3),(5,5)\nPhase Diff for ℓ=5 modes = ",FixedPhaseDiff]

ℓ=6 CHECK:  (6,2), (6,4),(6,6) have p5’s so phases agree.
(6,1), (6,3),(6,5) have p6’s which contain p2 thus the ± ambiguity

kvals=Range[15,20];
PhaseDiff=MyPhases[[kvals]]-Cphases[[kvals]];
Print["Phase Diff for ℓ=6 modes = ",PhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->0]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["No minus sign in p2 which goes into p6, (6,1), (6,3), (6,5) BROKEN\nPhase Diff for ℓ=6 modes = ",FixedPhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->1]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["Extra minus sign into p2->p6 FIXES (6,1),(6,3),(6,5)\nPhase Diff for ℓ=6 modes = ",FixedPhaseDiff]

ℓ=7 CHECK:  (7,2),(7,4),(7,6) contain p7 which is OFF by -1 so π phase difference
(7,1), (7,3),(7,5),(7,7) have p6’s which contain p2 thus the ± ambiguity

kvals=Range[21,27];
PhaseDiff=MyPhases[[kvals]]-Cphases[[kvals]];
Print["Phase Diff for ℓ=7 modes = ",PhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->0]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["No minus sign in p2 which goes into p6, (7,2), (7,4), (7,6) BROKEN\nPhase Diff for ℓ=7 modes = ",FixedPhaseDiff]
FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*Cphases[[kvals]]]/.r->1]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["Extra minus sign into p2->p6 FIXES (7,2), (7,4), (7,6)\nPhase Diff for ℓ=7 modes = ",FixedPhaseDiff]

l=8 PHASES
For this there is no data in any of the codes, but my Mathematica output above combined with sign change of p7 and ±ambiguity  in p8 can be used to obtain the final result for Waveform.c:40
It turns out adding π to the Mathematica output gives the phase compatible with Waveform.c

kvals=Range[28,35];
MyPhases[[kvals]]
Testl8Phases=π*Mod[((MyPhases[[kvals]]+π)/π),2]

CphasesNew=Flatten[Join[{Cphases[[1;;27]]},{Testl8Phases}]];

FixedPhaseDiff=Simplify[Arg[Alm[[kvals]]*Exp[I*CphasesNew[[kvals]]]/.r->0]-Arg[clms[[kvals]]*Exp[I*MyPhases[[kvals]]]]];
Print["Re-constructed phases for C, (8,1) ..(8,8)\nPhase Diff for ℓ=8 modes = ",FixedPhaseDiff]

REMARK: Matlab did have Phase=π for (8,8) and my above reconstructed (8,8) also is π  hence providing ONE consistency check

Cphases[[-1]]-CphasesNew[[-1]]

EXTRA

Mod[((MyPhases[[3;;5]]+Arg[Simplify[((clms[[3;;5]])/(Alm[[3;;5]]))]]-Cphases[[3;;5]])/(2π))/.r->1,1]

Block[{k1=6,k2=9},Mod[((MyPhases[[k1;;k2]]+Arg[Simplify[((clms[[k1;;k2]])/(Alm[[k1;;k2]]))]]-Cphases[[k1;;k2]])/(2π)),1]]/.r->1

Block[{k1=10,k2=14},Mod[((MyPhases[[k1;;k2]]+Arg[Simplify[((clms[[k1;;k2]])/(Alm[[k1;;k2]]))]]-Cphases[[k1;;k2]])/(2π)),1]]

Block[{k1=15,k2=20},Mod[((MyPhases[[k1;;k2]]+Arg[Simplify[((clms[[k1;;k2]])/(Alm[[k1;;k2]]))]]-Cphases[[k1;;k2]])/(2π)),1]]

Block[{k1=21,k2=26},Mod[((MyPhases[[k1;;k2]]+Arg[Simplify[((clms[[k1;;k2]])/(Alm[[k1;;k2]]))]]-Cphases[[k1;;k2]])/(2π)),1]]

Block[{k1=34,k2=34},Mod[((MyPhases[[k1;;k2]]+Arg[Simplify[((clms[[k1;;k2]])/(Alm[[k1;;k2]]))]]-Cphases[[k1;;k2]])/(2π)),1]]