# Surface plot for matrices (Matlab-Mathematica)

Now I am try to convert matlab file to mathematica,

This is the original Matlab command, surf, for making surface plots:

surf (xei , yei , zei , EE,'FaceColor ','interp ','FaceAlpha ' ,1.0 , 'LineWidth ' ,... 301
0.01 , 'EdgeAlpha ' ,0.2 , 'EdgeColor ','w');


If I want to plot the same matrix in Mathematica, which correct Plot Function shall I use?

For a better understanding, I put some convert codes here.

n = 70;
m = n/2;
CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0,
0}, {63.9, 63.9, 165.7, 0, 0, 0}, {0, 0, 0, 2*79.6, 0, 0}, {0, 0,
0, 0, 2*79.6, 0}, {0, 0, 0, 0, 0, 2*79.6}};
phi = N[Range[0, Pi, Pi/(m - 1)]];
theta = N[Range[0, 2*Pi, 2*Pi/(n - 1)]];
EE = ConstantArray[0, {m, n}];

xei = ConstantArray[0, {m, n}];
yei = ConstantArray[0, {m, n}];
zei = ConstantArray[0, {m, n}];

SS = Inverse[CC];

For[i = 1, i < n + 1, i++, For[j = 1, j < m + 1, j++,
xhen = {{Sin[phi[[j]]]*Cos[theta[[i]]],
Sin[phi[[j]]]*Sin[theta[[i]]], Cos[phi[[j]]]}};
(*xheng=Transpose[xhen];*)
xiaoa = KroneckerProduct[Transpose[xhen], xhen];
dbV = {{xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]],
Sqrt[2]*xiaoa[[2, 3]], Sqrt[2]*xiaoa[[1, 3]],
Sqrt[2]*xiaoa[[1, 2]]}};
(*dbV=Transpose[dbV];*)
FFXX = SS . Transpose[dbV];
(*FXXX=Transpose[dbV];*)
EE[[j, i]] = 1/(dbV . FFXX);]]

For[i = 1, i < n + 1, i++, For[j = 1, j < m + 1, j++,
xei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Cos[theta[[i]]];
yei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Sin[theta[[i]]];
zei[[j, i]] = EE[[j, i]]*Cos[phi[[j]]];]]
ListPlot3D[Transpose[{xei, yei, zei}]]

• Try ListPlot3D[Transpose[{xei, yei, zei}]]. Jul 8, 2018 at 13:25
• @HenrikSchumacher, xei,yei,zei are all not vector/list, if I use your method, then ListPlot3D::arrayerr: {<<1>>} must be a valid array or a list of valid arrays. Jul 8, 2018 at 13:34
• So, what is xei, yei, zei? Example data is always welcome. Jul 8, 2018 at 13:37
• Hey @HenrikSchumacher, I put complete code here, thanks a lot for your support! Jul 8, 2018 at 13:47
• Try also ListSurfacePlot3D (along with Flatten as in Anton Antonov's answer). Jul 8, 2018 at 16:20

Edit 4

We can project it on xy plane by setting z=0.

n = 250;
m = n/2;
CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0,
0}, {63.9, 63.9, 165.7, 0, 0, 0}, {0, 0, 0, 2*79.6, 0, 0}, {0, 0,
0, 0, 2*79.6, 0}, {0, 0, 0, 0, 0, 2*79.6}};
phi = N[Range[0, Pi, Pi/(m - 1)]];
theta = N[Range[0, 2*Pi, 2*Pi/(n - 1)]];
EE = xei = yei = zei = ConstantArray[0, {m, n}];
SS = Inverse[CC];

Table[xhen = {Sin[phi[[j]]]*Cos[theta[[i]]],
Sin[phi[[j]]]*Sin[theta[[i]]], Cos[phi[[j]]]};
xiaoa = KroneckerProduct[xhen, xhen];
dbV = {xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]],
Sqrt[2]*xiaoa[[2, 3]], Sqrt[2]*xiaoa[[1, 3]],
Sqrt[2]*xiaoa[[1, 2]]};
FFXX = SS.dbV;
EE[[j, i]] = 1/(dbV.FFXX), {i, n}, {j, m}];

Table[xei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Cos[theta[[i]]];
yei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Sin[theta[[i]]];
zei[[j, i]] = EE[[j, i]]*Cos[phi[[j]]];, {i, n}, {j, m}];

data = Transpose[Flatten /@ {xei, yei, zei, EE}];
data[[All, 3]] = 0;
color = Blend["ThermometerColors", Rescale[#, MinMax@EE]] &;

Legended[Graphics3D[{color[Last[#]], AbsolutePointSize[20],
Point[Most[#]]} & /@ data, ImageSize -> Large,
ViewPoint -> Above, Boxed -> False],
BarLegend[{color, MinMax@EE}, LegendMarkerSize -> {22, 400}]]


Or

data = Transpose[Flatten /@ {xei, yei, zei, EE}];
data = data[[All, {1, 2, 4}]];

Legended[Graphics[{color[Last[#]], AbsolutePointSize[20],
Point[Most[#]]} & /@ data, ImageSize -> Large],
BarLegend[{color, MinMax@EE}, LegendMarkerSize -> {22, 400}]]


Or you achieve this by using @xzczd

CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0,
0}, {63.9, 63.9, 165.7, 0, 0, 0}, {0, 0, 0, 2 79.6, 0, 0}, {0, 0,
0, 0, 2 79.6, 0}, {0, 0, 0, 0, 0, 2 79.6}};
SS = Inverse[CC];
xhen = {Sin[phi] Cos[theta], Sin[phi] Sin[theta], Cos[phi]};
xiaoa = Outer[Times, xhen, xhen];
dbV = {xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]],
Sqrt[2] xiaoa[[2, 3]], Sqrt[2] xiaoa[[1, 3]],
Sqrt[2] xiaoa[[1, 2]]};
EE = 1/dbV.SS.dbV;

SphericalPlot3D[EE, {phi, 0, Pi}, {theta, 0, 2 Pi},
ColorFunction ->
Function[{x, y, z, phi, theta, r},
ColorData["ThermometerColors"][r]], ViewPoint -> Above,
Boxed -> False, Axes -> False, PlotPoints -> 100, Mesh -> None]


Edit 3

You can increase number of point by increasing n.

 {col, row} = ImageDimensions@parula;
ParulaMMA =
Module[{colorlist},
colorlist =
Catenate@
ImageData@ImageTake[parula, {Round[row/2], Round[row/2]}, All];
Evaluate[Blend[RGBColor @@@ colorlist, Rescale[#, MinMax@EE]] &]];

n = 250;
m = n/2;
CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0,
0}, {63.9, 63.9, 165.7, 0, 0, 0}, {0, 0, 0, 2*79.6, 0, 0}, {0, 0,
0, 0, 2*79.6, 0}, {0, 0, 0, 0, 0, 2*79.6}};
phi = N[Range[0, Pi, Pi/(m - 1)]];
theta = N[Range[0, 2*Pi, 2*Pi/(n - 1)]];
EE = xei = yei = zei = ConstantArray[0, {m, n}];
SS = Inverse[CC];

Table[
xhen = {Sin[phi[[j]]]*Cos[theta[[i]]],
Sin[phi[[j]]]*Sin[theta[[i]]], Cos[phi[[j]]]};
xiaoa = KroneckerProduct[xhen, xhen];
dbV = {xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]],
Sqrt[2]*xiaoa[[2, 3]], Sqrt[2]*xiaoa[[1, 3]],
Sqrt[2]*xiaoa[[1, 2]]};
FFXX = SS.dbV;
EE[[j, i]] = 1/(dbV.FFXX), {i, n}, {j, m}];

Table[xei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Cos[theta[[i]]];
yei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Sin[theta[[i]]];
zei[[j, i]] = EE[[j, i]]*Cos[phi[[j]]];, {i, n}, {j, m}];

data = Transpose[Flatten /@ {xei, yei, zei, EE}];

Legended[Graphics3D[{ParulaMMA[Last[#]], AbsolutePointSize[20],
Point[Most[#]]} & /@ data, Axes -> True, ImageSize -> Large,
PlotRange -> {{-200, 200}, {-200, 200}, {-200, 200}}],
BarLegend[{ParulaMMA, MinMax@EE}, LegendMarkerSize -> {22, 400}]]


Or

Mathematica's built in gradient color

color = Blend["ThermometerColors", Rescale[#, MinMax@EE]] &;

Legended[Graphics3D[{color[Last[#]], AbsolutePointSize[20],
Point[Most[#]]} & /@ data, Axes -> True, ImageSize -> Large,
PlotRange -> {{-200, 200}, {-200, 200}, {-200, 200}}],
BarLegend[{color, MinMax@EE}, LegendMarkerSize -> {22, 400}]]


Edit 2 Here is another way to visualize it. see this for parula color.

data = Transpose[Flatten /@ {xei, yei, zei, EE}];


{col, row} = ImageDimensions@parula;
ParulaMMA =
Module[{colorlist},
colorlist =
Catenate@
ImageData@ImageTake[parula, {Round[row/2], Round[row/2]}, All];
Evaluate[Blend[RGBColor @@@ colorlist, Rescale[#, MinMax@EE]] &]];

Legended[Graphics3D[{ParulaMMA[Last[#]], Sphere[Most[#], 20]} & /@
data, Axes -> True, ImageSize -> Large,
PlotRange -> {{-200, 200}, {-200, 200}, {-200, 200}}],
BarLegend[{ParulaMMA, MinMax@EE}, LegendMarkerSize -> {22, 400}]]


This is the closest I can get.

data = Transpose[Flatten /@ {xei, yei, zei, EE}];

ListPointPlot3D[List /@ Most /@ data,
PlotStyle -> ({AbsolutePointSize[22],
Blend["BlueGreenYellow", Rescale[#, MinMax@EE]]} & /@
Last /@ data), BoxRatios -> {1, 1, 1}]


• Thanks a lot for figuring this out! It looks much better, but surface is not smoothly in Mathematica. Jul 9, 2018 at 8:11
• If I run this plot, I get bugs, namely: ImageDimensions::imginv: Expecting an image or graphics instead of parula. Set::shape: Lists {col,row} and ImageDimensions[parula] are not the same shape. ImageTake::imginv: Expecting an image or graphics instead of parula. Jul 9, 2018 at 15:09
• You need to define parula. got to here mathworks.com/help/matlab/ref/colormap.html and copy parula and paste it in Mathematica notebook Jul 9, 2018 at 15:22
• Great, thanks a lot, this problem is solved! Jul 9, 2018 at 15:43
• @Okkes Dulgerci, is there also a simple command that we can get a 2D results(projection:3D-2D) Aug 10, 2018 at 11:05

I should have read your question and the surf help page more carefully -- your are using four matrices {xei, yei, zei, EE}.

Here is one way to use EE for the surface color:

colorRules =
Most[#] -> Last[#] & /@ Transpose[Flatten /@ {xei, yei, zei, EE}];

ListSurfacePlot3D[Transpose[Flatten /@ {xei, yei, zei}],
ColorFunction -> Function[{x, y, z}, {x, y, z} /. colorRules]]


I assume this can be improved.

## Second try

I understand that the following plots are "volume" density and contour plots, not a surface density plot which you want. But I think they are close...

Using ListDensityPlot3D:

ListDensityPlot3D[Transpose[Flatten /@ {xei, yei, zei, EE}]]


Another plot that is close is made with ListContourPlot3D:

ListContourPlot3D[Transpose[Flatten /@ {xei, yei, zei, EE}],
Contours -> 30, Mesh -> None, PlotLegends -> "Expressions"]


## First take

(More of a comment to find out what OP wants...)

Evaluating the code you posted we have these dimensions for {xei, yei, zei}:

Dimensions /@ {xei, yei, zei}
(* {{35, 70, 1, 1}, {35, 70, 1, 1}, {35, 70, 1, 1}} *)


I strongly suspect that what you want is this plot:

ListPlot3D[Transpose[Flatten /@ {xei, yei, zei}], PlotRange -> All]


• But It still contains the bugs as the cube is not closed form. looking for your reply! Thanks Jul 8, 2018 at 20:37
• the mesh gitter also not correct. Jul 8, 2018 at 20:52
• why 'ListSurfacePlot3D' is not correct in Mathematica? Jul 9, 2018 at 10:29
• @FEAPMAN "why 'ListSurfacePlot3D' is not correct in Mathematica?" -- I do not know is it correct or not. Meaning, ListSurfacePlot3D might have a bug, or the problem might be in your data and a particular "feature" of ListSurfacePlot3D... It is a good idea to investigate, but I am not doing that anytime soon. Jul 9, 2018 at 14:06

OK, though OP doesn't reply to my comment, I'd like to post a possible answer. If OP just wants to obtain that plot and "plot the same matrix" is not necessary, then the solution is quite simple. We don't need to discretize the analytic solution and transform the coordinate etc., because we have SphericalPlot3D:

CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0, 0}, {63.9, 63.9, 165.7, 0,
0, 0}, {0, 0, 0, 2 79.6, 0, 0}, {0, 0, 0, 0, 2 79.6, 0}, {0, 0, 0, 0, 0, 2 79.6}};
SS = Inverse[CC];
xhen = {Sin[phi] Cos[theta], Sin[phi] Sin[theta], Cos[phi]};
xiaoa = Outer[Times, xhen, xhen];
dbV = {xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]], Sqrt[2] xiaoa[[2, 3]],
Sqrt[2] xiaoa[[1, 3]], Sqrt[2] xiaoa[[1, 2]]};
EE = 1/dbV.SS.dbV;

SphericalPlot3D[EE, {phi, 0, Pi}, {theta, 0, 2 Pi},
ColorFunction -> Function[{x, y, z, phi, theta, r}, ColorData["AvocadoColors"][r]]]


I believe the definition of EE can be simplified further, but without background information it's somewhat beyond my reach.

• this plot results can show us the Anisotropic Elasticity (similar to the work: "Anisotropic Elasticity" demonstrations.wolfram.com/AnisotropicElasticity Wolfram Demonstrations Project Published: March 7 2011) Jul 10, 2018 at 12:10
• The original code is wrotten in Matlab which needs to discretize the analytic solution and transform the coordinate etc. Now Mathematica show its strongly numerical ability in this case, cheers! Jul 10, 2018 at 12:13
• @FEAPMAN Is EE defined by some kind of rotation? Can you elaborate a little on this part? Jul 10, 2018 at 15:05
• EE is elasticity tensor. Aug 10, 2018 at 11:03
• @ABCDEMMM I know, but what I'm interested in is, why is it defined in this way?xiaoa and dbV` looks like something related to rotation or coordinate transform, and I believe they can be simplified, but without more background information it's beyond my reach. Aug 11, 2018 at 14:21