transmittance[x_, σe_] :=
{
Exp[-σe * x] 
}

scattering[x_] :=
{
σs * transmittance[x, σe] 
}

(*Improved integration of (Section 5.6.3)*)
Integrate[scattering[xx],{xx,0,d}]

(*Schlick phase function with g remapping into k (Section 2.3)*)
phaseFuncSchlick[G_, A_] :=
{
k = 1.55*G - 0.55*G*G*G;
tmp = 1.0 + k*Cos[A];
-(1-k*k)/((4.0*π*tmp*tmp))
}[[1]]

phaseHG[G_, A_] :=
{
(1-G*G)/(4.0*Pi*(1+G*G-2*G*Cos[A])^1.5)
}[[1]]

phaseRayl[A_] :=
{
3*(1+Cos[A]*Cos[A])/(16.0*Pi)
}[[1]]

PolarPlot[phaseRayl[ a],{a,0,2Pi}, PlotRange->{{-0.15,0.15},{-0.1,0.1}}]

Show[
PolarPlot[phaseHG[0.0, a],{a,0,2Pi}, PlotStyle->Blue],
PolarPlot[phaseHG[0.3, a],{a,0,2Pi}, PlotStyle->Blue] ,
PolarPlot[phaseHG[0.6, a],{a,0,2Pi}, PlotStyle->Blue],
PolarPlot[phaseFuncSchlick[0.0, a],{a,0,2Pi}, PlotStyle->Red],
PolarPlot[phaseFuncSchlick[0.3, a],{a,0,2Pi}, PlotStyle->Red] ,
PolarPlot[phaseFuncSchlick[0.6, a],{a,0,2Pi}, PlotStyle->Red],
PlotRange->All
]

Show[
PolarPlot[phaseHG[0.9, a],{a,0,2Pi}, PlotStyle->Blue],
PolarPlot[phaseFuncSchlick[0.9, a],{a,0,2Pi}, PlotStyle->Red],
PlotRange->{{-0.2,27},{-1,1}}
]

Show[
Plot[transmittance[x, 0.5],{x,-0,5}, PlotStyle->Blue],
Plot[transmittance[x, 1.0],{x,-1,5}, PlotStyle->Red],
Plot[transmittance[x, 2.0],{x,-1,5}, PlotStyle->Orange],
PlotRange->{{-0.5,5},{-0.05,1.0}},
AxesOrigin->{0,0}
]

(*Integral of the phase function over a sphere*)
Integrate[
phaseHG[0.0, θ] *Sin[θ]
,{θ,0,π},{ϕ,0,2 Pi}]

(*Integral of cos(a) over hemisphere, show why luminance<->illuminance conversion only has π factor, because illuminance already include integral with Cos[θ] visibility/projected angle*)
Integrate[
Cos[θ] *Sin[θ]
,{θ,0,π/2},{ϕ,0,2π}]

(*Show how the exponential curve behave as compared to the linear interpolation of transmittance (to reach t=0.2 over 1meter). This seems simple but is the reason the trapozoidal rule is not enough to be energy conservative when integrating scattered luminance.*)
Show[
Plot[1-0.8*x,{x,0,1}, PlotStyle->Orange],
Plot[transmittance[x, -Log[0.2]],{x,0,1}, PlotStyle->Green],
PlotRange->{{-0.05,1.05},{-0.05,1.05}},
AxesOrigin->{0,0}
]