# Building a 2D function by combining 1D functions

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.953065429366962E^(-25.926092562282516 \ (-0.0017564449181118602 + xdash)^2) + 0.9338124308343011 E^(-2.8693935369154184 (0.195211521347951 + xdash)^2));
OneDimFunc2 = 0.03136086274687064 (30.953065429366962E^(-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?

• Since your original 1D slice is symmetric and you just want to revolve it around the z-axis in you original plot you just have an expression in polar coordinates where the slice from $(0, \infty)$ of your 1D is $r$ and the $\phi$ is as ususal. Mar 19, 2018 at 17:42
• The code that I have posted here uses just one 1D function (for the sake of simplicity) so maybe you understood that my 1D slice is symmetric. But there are six different 1D functions that I would like to use and those are not symmetric. Therefore, I cannot just revolve it around the $z$-axis. Mar 19, 2018 at 19:03
• Do you have any idea how they depend parametrically on $\varphi$? If so you can include that in an analytic expression and then convert from polar to cartesian. Otherwise I think you'll have to fake it via interpolation of some form. Mar 19, 2018 at 19:06
• You might want to look up transfinite interpolation. Mar 19, 2018 at 19:36
• @b3m2a1 I just have these equations on their own terms. I do not have any idea about their parametric dependence on the rotation angle $\phi$. Mar 20, 2018 at 0:13

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

• How do I interpolate twoDdata with InterpolationOrder-> 3?. It shows an error that " Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All" Mar 22, 2018 at 11:43
• @dykes, i don't know how to do interpolation on unstructured grids. Maybe some of these Q/As may be useful.
– kglr
Mar 22, 2018 at 20:34
• It's almost working : Show[ParametricPlot3D[Evaluate@oneDfuncs[x], {x, -1, 1}], ListSurfacePlot3D[ Flatten[Transpose[(oneDfuncs /@ Range[-1, 1, 2/99])], 1], MaxPlotPoints -> 80, BoxRatios -> {1, 1, 1/2}, ColorFunction -> Function[{x, y, z}, Hue[z]]]] Mar 22, 2018 at 23:25