Building a 2D function by combining 1D functions

Question

I have six 1D continuous functions that I want to use in order to find a 2D function. Each of these 1D functions has to be placed at a unique angle $\phi$ w.r.t the $x$-axis as given in the image below. In the image, for example, a 1D function in black color is at an angle $\phi = -1.54337$ rad and the red colored 1D function is at an angle $\phi = 0.5$ rad.

This way $n$ 1D functions can be used. Here, I have used six different 1D functions and have placed them at different angles. I'm quite new to Mathematica, therefore, my code is a little messy. The Edited code to construct the data is here:

OneDimFunc1 = 0.03136086274687064` (30.953065429366962`E^(-25.926092562282516` \ (-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-2.8693935369154184` (0.195211521347951` + xdash)^2));
OneDimFunc2 = 0.03136086274687064` (30.953065429366962`E^(-22 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-2 (0.195211521347951` + xdash)^2));
OneDimFunc3 = 0.03136086274687064` (30.953065429366962` E^(-15 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-2.5 (0.195211521347951` + xdash)^2));
OneDimFunc4 = 0.03136086274687064` (30.953065429366962` E^(-27 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-2.9 (0.195211521347951` + xdash)^2));
OneDimFunc5 = 0.03136086274687064` (30.953065429366962` E^(-17 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-1.8 (0.195211521347951` + xdash)^2));
OneDimFunc6 = 0.03136086274687064` (30.953065429366962` E^(-23 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-3.2 (0.195211521347951` + xdash)^2));

NumData = 100;         (* Number of data points *)
xD = Array[# &, NumData, {-1, 1}];

zDt1 = Flatten[Abs[OneDimFunc1] /. {xdash -> xD}]; 
x2Ddata1 = Array[# &, NumData, {-1, 1}] Cos[-1.54337];  (* \[Phi] = -1.54337*)
y2Ddata1 = Array[# &, NumData, {-1, 1}] Sin[-1.54337];
Zdt1 = Transpose[{x2Ddata1, y2Ddata1, zDt1}]; 

zDt2 = Flatten[Abs[OneDimFunc2] /. {xdash -> xD}];  
x2Ddata2 = Array[# &, NumData, {-1, 1}] Cos[-1];   (* \[Phi] = 1*)
y2Ddata2 = Array[# &, NumData, {-1, 1}] Sin[-1];
Zdt2 = Transpose[{x2Ddata2, y2Ddata2, zDt2}];

zDt3 = Flatten[Abs[OneDimFunc3] /. {xdash -> xD}]; 
x2Ddata3 = Array[# &, NumData, {-1, 1}] Cos[-0.5];  (* \[Phi] = -0.5*)
y2Ddata3 = Array[# &, NumData, {-1, 1}] Sin[-0.5];
Zdt3 = Transpose[{x2Ddata3, y2Ddata3, zDt3}];

zDt4 = Flatten[Abs[OneDimFunc4] /. {xdash -> xD}];          
x2Ddata4 = Array[# &, NumData, {-1, 1}] Cos[-0];   (* \[Phi] = 0*)
y2Ddata4 = Array[# &, NumData, {-1, 1}] Sin[-0]; 
Zdt4 = Transpose[{x2Ddata4, y2Ddata4, zDt4}];

zDt5 = Flatten[Abs[OneDimFunc5] /. {xdash -> xD}];     
x2Ddata5 = Array[# &, NumData, {-1, 1}] Cos[0.5];  (* \[Phi] = 0.5*)
y2Ddata5 = Array[# &, NumData, {-1, 1}] Sin[0.5];
Zdt5 = Transpose[{x2Ddata5, y2Ddata5, zDt5}];

zDt6 = Flatten[Abs[OneDimFunc6] /. {xdash -> xD}];    
x2Ddata6 = Array[# &, NumData, {-1, 1}] Cos[1];  (* \[Phi] = 1*)
y2Ddata6 = Array[# &, NumData, {-1, 1}] Sin[1]; 
Zdt6 = Transpose[{x2Ddata6, y2Ddata6, zDt6}];

TwoDimdata = Join[Zdt1, Zdt2, Zdt3, Zdt4, Zdt5, Zdt6];  (* Combined data points *)
ListPointPlot3D[TwoDimdata]

After combining the data I fit 2D Gaussian functions to the data using nonlinear fitting as follows:

NBasisFunc2D = 3;
TwoDimModel = Sum[a[i] Exp[-( ((xout - center1[i])/b[i])^2 + ((yout - center2[i])/b[i])^2)], {i, 1, NBasisFunc2D, 1}];  (* Bivariate Gaussians *)
TwoDimParameters = Flatten[Transpose[{Join[Array[a, NBasisFunc2D][[;; ;; 1]], Array[b, NBasisFunc2D][[;; ;; 1]], Array[center1, NBasisFunc2D][[;; ;; 1]], Array[center2, NBasisFunc2D][[;; ;; 1]]]}]];
nlm2D = NonlinearModelFit[TwoDimdata, TwoDimModel, TwoDimParameters, {xout, yout}, MaxIterations -> 1000];

FittedPrams = nlm2D["BestFitParameters"]      

Show[Plot3D[TwoDimModel /. FittedPrams, {xout, -0.45, 0.45}, {yout, -0.45, 0.45}], Graphics3D[{Red, Point /@ TwoDimdata}]]

(Edit: The 2D Gaussian fit is poor after using six different 1D funtions)

I am not sure whether this is a correct way to find out the 2D function. Can I directly use the 1D continuous functions to find out the 2D function?. If yes then how do I do it?. Or should I use cubic spline interpolation instead of least squares fit?.

UPDATE:

I tried the Matlab griddata to interpolate my data. From Matlab documentation: vq = griddata(x,y,v,xq,yq) fits a surface of the form $v = f(x,y)$ to the scattered data in the vectors $(x,y,v)$. The griddata function interpolates the surface at the query points specified by $(xq,yq)$ and returns the interpolated values, $vq$. The Matlab code is:

[xq,yq] = meshgrid(-1:0.01:1);
vq = griddata(x,y,Fxy,xq,yq,'cubic');
surface(xq,yq,vq)
hold on;
plot3(x,y,Fxy,'or','MarkerFaceColor', 'r')
hold off;

where x,y and Fxy is taken from TwoDimdata in Mathematica:

x = TwoDimdata[[All, 1]] ;
y = TwoDimdata[[All, 2]];
Fxy = TwoDimdata[[All, 3]];  

When I try the same in Mathematica I recieve an error:

Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All

How do I do this Matlab equivalent interpolation in Mathematica?

Answer 1

This is an alternative way to construct your data:

ϕ = {-1.54337, -1, -0.5, 0, .5, 1.};
aa = {-25.926092562282516`, -22, -15, -27, -17, -23};
bb = {-2.8693935369154184`, -2, -2.5, -2.9, -1.8, -3.2};
{α, β, γ, δ, ρ} =  {0.03136086274687064, 30.953065429366962`, -0.0017564449181118602`, 
   0.9338124308343011`, 0.195211521347951`};
cc = Transpose[{ϕ, aa, bb}];
NumData = 100;
ClearAll[oneDfuncs]
oneDfuncs[x_] := {x Cos@#, x Sin@#, Abs[α(β E^(#2 (γ + x)^2) + δ E^(#3 (ρ + x)^2))]}&@@@cc;
twoDdata = Transpose[oneDfuncs /@ Range[-1, 1, 2/99]];

Show[ParametricPlot3D[Evaluate@oneDfuncs[x], {x, -1, 1}, 
  PlotLegends -> ("OneDimFunc" <> ToString[#] & /@ Range[6])], 
 ListPointPlot3D[twoDdata]]

The data generated is the same as your TwoDimdata

Join @@ twoDdata == TwoDimdata

True