(* Functions to compute the SH vectors *)
K[l_,m_]:=Sqrt[(2*l+1)((l-Abs[m])!)/(4*Pi * ((l+Abs[m])!))];
SH[l_, m_, t_, p_] := Piecewise[{{Sqrt[2]*K[l,m]*Cos[m*p]*LegendreP[l,m,Cos[t]],m>0},{Sqrt[2]*K[l,m]*Sin[Abs[m]*p]*LegendreP[l,Abs[m],Cos[t]],m<0},{K[l,0]*LegendreP[l,0,Cos[t]],m==0}}];

(* Phong BRDF *)
Phong[r_, t_] = (r+1)/ (2*Pi) * Cos[t]^r;

(* Those functions compute the (log) SH vectors of the Phong BRDF :
   brdfSH : classic projection of the BRDF in the SH basis;
     brdfLogSHEps : projection in log-SH, depending on the parameter epsilon;
     brdfLogSHEpsHemi : projection in log-SH, ignoring the lower hemisphere;
     brdfLogSHEpsNotNormalized : projection of cos^r in log-SH (that is, unnormalized Phong), considering the whole sphere.
 *)
brdfSH[r_] := Flatten[N[Table[Integrate[Phong[r,t]Sin[t]SH[l,m,t,p],{t, 0, Pi/2},{p,0, 2*Pi}], {l,0,3}, {m,-l,l}]]]
brdfLogSHEps[r_,epsilon_] :=Flatten[N[Table[
Integrate[Log[Max[Phong[r,t], epsilon]]Sin[t]SH[l,m,t,p],{t,0,Pi/2},{p,0,2*Pi}] + Integrate[Log[epsilon]Sin[t]SH[l,m,t,p],{t,Pi/2,Pi},{p,0,2*Pi}],
{l,0,3},{m,-l,l}]]];
brdfLogSHEpsHemi[r_,epsilon_] :=Flatten[N[Table[
Integrate[Log[Max[Phong[r,t], epsilon]]Sin[t]SH[l,m,t,p],{t,0,Pi/2},{p,0,2*Pi}] ,
{l,0,3},{m,-l,l}]]];
brdfLogSHEpsNotNormalized[r_,epsilon_] :=Flatten[N[Table[
Integrate[Log[Max[Cos[t]^r, epsilon]]Sin[t]SH[l,m,t,p],{t,0,Pi/2},{p,0,2*Pi}] + Integrate[Log[epsilon]Sin[t]SH[l,m,t,p],{t,Pi/2,Pi},{p,0,2*Pi}],
{l,0,3},{m,-l,l}]]];

(* Here I precompute the LogSH of Phong BRDF at discrete values of r and epsilon. I write the results in different files for each value of epsilon. Then I read the file in a C program and compute Exp(LogSH(BRDF)). I write the results of the C program in the corresponding files. *)
Do[
LogBRDFTable = Table[brdfLogSHEps[r,epsilon],{r,1,5}];
fileLog = OpenWrite[NotebookDirectory[]<>"fileLogBRDFcomp_"<>ToString[epsilon]<>".txt"];
Do[Do[WriteString[fileLog,TextString[LogBRDFTable[[i,j]]] <> " "], {j,1,16}] WriteString[fileLog,"\n"],{i,1,5} ];
Close[NotebookDirectory[]<>"fileLogBRDFcomp_"<>ToString[epsilon]<>".txt"];
Print[ToString[epsilon]<>" done."];
,{epsilon, { 0.01, 0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}}];

(* Same thing for brdfLogSHEpsHemi *)
Do[
LogBRDFTable = Table[brdfLogSHEpsHemi[r,epsilon],{r,1,5}];
fileLog = OpenWrite[NotebookDirectory[]<>"fileLogBRDFhemi_"<>ToString[epsilon]<>".txt"];
Do[Do[WriteString[fileLog,TextString[LogBRDFTable[[i,j]]] <> " "], {j,1,16}] WriteString[fileLog,"\n"],{i,1, 5} ];
Close[NotebookDirectory[]<>"fileLogBRDFhemi_"<>ToString[epsilon]<>".txt"];
Print[ToString[epsilon]<>" done."];
,{epsilon, { 0.01, 0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}}];

(* Same thing for brdfLogSHEpsNotNormalized *)
Do[
LogBRDFTableNotNormalized = Table[brdfLogSHEpsNotNormalized[r,epsilon],{r,{1,2,3,4,5}}];
fileLog = OpenWrite[NotebookDirectory[]<>"fileLogBRDFnotNormalized_"<>ToString[epsilon]<>".txt"];
Do[Do[WriteString[fileLog,TextString[LogBRDFTableNotNormalized[[i,j]]] <> " "], {j,1,16}] WriteString[fileLog,"\n"],{i,1, 5} ];
Close[NotebookDirectory[]<>"fileLogBRDFnotNormalized_"<>ToString[epsilon]<>".txt"];
Print[ToString[epsilon]<>" done."];
,{epsilon, { 0.01, 0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}}];

(* Here I read back the output of the C program that has applied the exponentiation on the previously generated files. This way I'll be able to compare normal SH with Exp(Log-SH).*)ExpLogBRDFTable = Table[ReadList[NotebookDirectory[]<>"fileExpLogBRDFcomp_"<>ToString[epsilon]<>".txt", Table[Number,{16}]]
,{epsilon, { 0.01, 0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}}];

ExpLogBRDFhemiTable = Table[ReadList[NotebookDirectory[]<>"fileExpLogBRDFhemi_"<>ToString[epsilon]<>".txt", Table[Number,{16}]]
,{epsilon, { 0.01, 0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}}];

ExpLogBRDFnnTable = Table[ReadList[NotebookDirectory[]<>"fileExpLogBRDFnotNormalized_"<>ToString[epsilon]<>".txt", Table[Number,{16}]]
,{epsilon, { 0.01, 0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}}];

(* Here is an example of the values I get for the basic Exp(Log-SH) *)
ExpLogBRDFTable

(* In order to compare my results, I compute the normal SH vectors for the Phong exponents 1 to 5 *)
SHBRDFTable = Table[brdfSH[r],{r,1,5}]

(* I compare the values using Euclidean distance. Here are the numerical values. *)
Table[EuclideanDistance[ExpLogBRDFTable[[e,#]], SHBRDFTable[[#]]], {e,1,19}]&/@{ 1,2,3,4,5}//Evaluate

(* And here the plot. *)
valuesX = {{1,0.01},{2,0.02}, {3,0.03}, {4,0.04}, {5,0.05}, {6,0.06}, {7,0.07}, {8,0.08}, {9,0.09}, {10,0.1}, {11,0.2}, {12,0.3}, {13,0.4},  {14,0.5}, {15,0.6}, {16,0.7},{17,0.8}, {18,0.9}, {19,1}};
rdata = {1,2,3,4,5};
p1 = ListLinePlot[Table[{valuesX[[e,2]],EuclideanDistance[ExpLogBRDFTable[[e,#]], SHBRDFTable[[#]]]}, {e,1,19}]&/@{ 1,2,3,4,5}//Evaluate, 
AxesLabel->{Style["ϵ",14],Style["(Subscript[L, 2]) error",14]},ImageSize->Large, Ticks->Automatic, 
PlotLegends->Placed[{Style["  r = "<>ToString[rdata[[#]]],14]&/@{ 1,2,3,4,5}},Below], PlotStyle->{ {RGBColor[0.996078431372549,0.3607843137254902,0.027450980392156862],AbsoluteThickness[3]}, {RGBColor[0.5411764705882353,0.7137254901960784,0.027450980392156862],AbsoluteThickness[3]},{RGBColor[0.15294117647058825,0.11372549019607843,0.49019607843137253],AbsoluteThickness[3]}, {RGBColor[0.00784313725490196,0.5098039215686274,0.9294117647058824],AbsoluteThickness[3]},{RGBColor[0.9058823529411765,0.027450980392156862,0.12941176470588237],AbsoluteThickness[3]}}, BaseStyle->{FontSize->14}, Joined-> True]

Export[NotebookDirectory[]<>"errorStudyGraph12345.pdf",p1];