# ContourPlot of a complicated function

I have the function $pot_1(x,y,z)$ which is given in analytical form. So, in order to plot the contours at the primary planes $(x,y)$, $(x,z)$ and $(y,z)$, I use the following Mathematica code:

Clear["Global"];
Off[General::spell];
pot1 = -Mn/Sqrt[x^2 + y^2 + z^2 + cn^2];
V = pot1;
Vxy = V /. {z -> 0};
Vxz = V /. {y -> 0};
Vyz = V /. {x -> 0};
Mn=0.08;cn=0.25;
S1 = ContourPlot[Vxy, {x, -10, 10}, {y, -10, 10},
FrameLabel -> {"x", "y"}, RotateLabel -> False, Contours -> 20,
ContourStyle -> Black, ContourShading -> False, PlotRange -> All]
S2 = ContourPlot[Vxz, {x, -10, 10}, {z, -10, 10},
FrameLabel -> {"x", "z"}, RotateLabel -> False, Contours -> 20,
ContourStyle -> Black, ContourShading -> False, PlotRange -> All]
S3 = ContourPlot[Vyz, {y, -10, 10}, {z, -10, 10},
FrameLabel -> {"y", "z"}, RotateLabel -> False, Contours -> 20,
ContourStyle -> Black, ContourShading -> False, PlotRange -> All]


OK, so far so good. Now, I have one more function $pot_2(x,y,z)$ which is much more complicated than $pot_1$ and we can obtain its value only on specific points $(x_0,y_0,z_0)$. This function is given by

pot2[G_, Mt_, a_, b_, c_, x_, y_, z_] := Module[{},
m[u_] := Sqrt[x^2/(a^2 + u) + y^2/(b^2 + u) + z^2/(c^2 + u)];
λ = If[m[0] > 1, u /. FindRoot[m[u]^2 == 1, {u, 1}], 0];
Δ[u_] := Sqrt[(a^2 + u) (b^2 + u) (c^2 + u)];
w000 = NIntegrate[1/Δ[u], {u, λ, \[Infinity]}];
w100 = NIntegrate[1/Δ[u] 1/(a^2 + u), {u, λ, \[Infinity]}];
w001 = NIntegrate[1/Δ[u] 1/(c^2 + u), {u, λ, \[Infinity]}];
w010 = 2/Δ[λ] - w100 - w001;
w110 = (w010 - w100)/(a^2 - b^2); w011 = (w001 - w010)/(b^2 - c^2);
w101 = (w100 - w001)/(c^2 - a^2);
w200 = 1/3 (2/(Δ[λ] (a^2 + λ)) - w110 - w101);
w020 = 1/3 (2/(Δ[λ] (b^2 + λ)) - w011 - w110);
w002 = 1/3 (2/(Δ[λ] (c^2 + λ)) - w101 - w011);
w111 = (w110 - w011)/(c^2 - a^2);
w120 = (w020 - w110)/(a^2 - b^2); w012 = (w002 - w011)/(b^2 - c^2);
w201 = (w200 - w101)/(c^2 - a^2);
w210 = (w110 - w200)/(a^2 - b^2); w021 = (w011 - w020)/(b^2 - c^2);
w102 = (w101 - w002)/(c^2 - a^2);
w300 = 1/5 (2/(Δ[λ] (a^2 + λ)^2) - w210 - w201);
w030 = 1/5 (2/(Δ[λ] (b^2 + λ)^2) - w021 - w120);
w003 = 1/5 (2/(Δ[λ] (c^2 + λ)^2) - w102 - w012);
cc = 15/16 G Mt;
pot = -(cc/6) (w000 - 6 x^2 y^2 z^2 w111 +
x^2 (x^2 (3 w200 - x^2 w300) +
3 (y^2 (2 w110 - y^2 w120 - x^2 w210) - w100)) +
y^2 (y^2 (3 w020 - y^2 w030) +
3 (z^2 (2 w011 - z^2 w012 - y^2 w021) - w010)) +
z^2 (z^2 (3 w002 - z^2 w003) +
3 (x^2 (2 w101 - x^2 w201 - z^2 w102) - w001)));
Print["pot2(x,y,z) = ", pot]
]


When $G=1$, $Mt=0.1$, $a=7$, $b=1.5$, $c=0.6$ and $-10 \leq x,y,z \leq 10$ I want to plot the contours on the three primary planes of the function $V(x,y,z) = pot1 + pot2$.

Any suggestions?

Having replaced the Print statement in the last line of your module by

pot]


and defined

pot1[x_, y_, z_] = -Mn/Sqrt[x^2 + y^2 + z^2 + cn^2] /. Mn -> 0.08 /. cn -> 0.25;


you just make a table out of it?

 datz = Table[pot2[1, 0.1, 9, 1.5, 0.6, x, y, 0]+pot1[ x, y, 0],
{x, -10, 10, 0.55}, {y, -10, 10, 0.55}];
ListContourPlot[datz, DataRange -> {{-10, 9.8}, {-10, 9.8}}]


If you want all 3 slices

  datz = Table[pot2[1, 0.1, 9, 1.5, 0.6, x, y, 0] +
pot1[x, y, 0], {x, -10, 10, 0.55}, {y, -10, 10, 0.55}, ColorFunction -> "Pastel"];
daty = Table[pot2[1, 0.1, 9, 1.5, 0.6, x, 0, y] +
pot1[x, 0, y], {x, -10, 10, 0.55}, {y, -10, 10, 0.55}, ColorFunction -> "Heat"];
datx = Table[pot2[1, 0.1, 9, 1.5, 0.6, 0, x, y] +
pot1[0, x, y], {x, -10, 10, 0.55}, {y, -10, 10, 0.55}, ColorFunction -> "Temperature"];


then

 {ListContourPlot[datx, DataRange -> {{-10, 9.8}, {-10, 9.8}}],
ListContourPlot[daty, DataRange -> {{-10, 9.8}, {-10, 9.8}}],
ListContourPlot[datz, DataRange -> {{-10, 9.8}, {-10, 9.8}}]} // GraphicsRow
`

does the trick

• Unfortunately, this is not what I want. Why did you set z=3? Also I want to plot the contours of the sum pot1 + pot2 at the three planes not just pot2 for a single value of z as you did. Mar 5, 2013 at 17:31
• I think we are approaching! However, when I create the tables dtax, daty and datz I get the following error message Mar 5, 2013 at 17:44
• FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.... Mar 5, 2013 at 17:44
• that's because your potential is singular at the origin... Mar 5, 2013 at 17:46
• I am not sure about pot2 but pot1 does not have a singularity at (0,0,0) since the constant cn is always present and non zero. Anyway, how can I overtake this issue? Mar 5, 2013 at 17:50