# How to obtain the correct interpolation for the following data?

I have a dataset DistrTemp.dat. It represents a set of data in the form {{x,y,z,f[x,y,z]}}. Here, x defines a system, while f[x,y,z] is a probability to find some system in the vicinity of y,z and has a property $$\int \limits_{0}^{\pi} dy \int \limits_{x}^{3500} dz \ f(x,y,z) = 1$$ My goal is to obtain continuous function F[x,y,z] that I can use for the integration measure.

I try to interpolate the dataset for particular values of x:

SetDirectory[NotebookDirectory[]]
DistrData = Import["DistrTemp.dat","Table"];
DistrInterpolation005[\[Theta]S_, ES_] =
Interpolation[{Log10[#[]], Log10[#[]], #[]} & /@
Select[DistrData, #[] == 0.05 &], InterpolationOrder -> 1][
Log10[\[Theta]S], Log10[ES]]


Then I compare the obtainted interpolation with the initial dataset:

Plot3D[DistrInterpolation005[\[Theta]S, 10^(ES)], {\[Theta]S, 0,
Pi}, {ES, Log10[0.05], Log10}, ImageSize -> Large,
PlotRange -> All]
ListPlot3D[{#[], Log10[#[]], #[]} & /@
Select[DistrData, #[] == 0.05 &], PlotRange -> All,
PlotStyle -> Directive[Red], ImageSize -> Large]


I find that the interpolation has been done incorrectly, since some domains are lost, while some domains are overestimated: I tried to use WeightedData:

TableCoordinates = Select[DistrData, #[] == 0.05 &][[All, {2, 3}]]
TableWeights = Flatten[DataTemp[[All, 4]]]
F1 = WeightedData[TableCoordinates, TableWeights]
DistrWeighted =
Table[f[F1], {f, {HistogramDistribution, SmoothKernelDistribution,
EmpiricalDistribution}}]
Table[Plot3D[
PDF[i, {\[Theta]S, Log10[ES]}], {ES, Log10[0.05],
Log10}, {\[Theta]S, 0, Pi},
PlotLabel -> Row[{"Mean: ", Mean[i]}], ImageSize -> Large], {i,
DistrWeighted}]


However, neither of the distributions match the desired distribution: How to construct the interpolation correctly?

Interpolation does fit all the points, which suggests that it probably is a plotting resolution issue. Try a higher number of PlotPoints. For example,

lp = ListPlot3D[{#[], Log10[#[]], #[]} & /@
Select[DistrData, #[] == 0.05 &], PlotRange -> All,
PlotStyle -> {Directive[Green], Opacity[0.5]}, ImageSize -> Large];
plt = Plot3D[
DistrInterpolation005[\[Theta]S, 10^(ES)], {\[Theta]S, 0, Pi}, {ES,
Log10[0.05], Log10}, ImageSize -> Large, PlotRange -> All,
PlotPoints -> 400];
Show[lp, plt] 