Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question

1 Answer 1

up vote 1 down vote accepted

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}}]

Mathematica graphics

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

Mathematica graphics

share|improve this answer
    
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. –  Vaggelis_Z Mar 5 '13 at 17:31
    
I think we are approaching! However, when I create the tables dtax, daty and datz I get the following error message –  Vaggelis_Z Mar 5 '13 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.... –  Vaggelis_Z Mar 5 '13 at 17:44
    
that's because your potential is singular at the origin... –  chris Mar 5 '13 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? –  Vaggelis_Z Mar 5 '13 at 17:50

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.