0
$\begingroup$

I would like to reproduce the solid blue line of Figure 2a in this article. For this, I have used the code found here. Notice, however, that I am not succeeding. Could you help me? Here is my code:

l[wx_, wz_] := wz/wx;

f[X_] := (1 + 2 X^2)/(1 - X^2) - 
   3 X^2 ArcTanh[Sqrt[1 - X^2]]/(1 - X^2)^(3/2);

yqf[h_, M_, as_, add_] := 
  32*4 Pi*h^2*as/(3 M) Sqrt[as^3/Pi] (1 + 3 add^2/(2 as^2));

az[h_, M_, wz_] := Sqrt[h/(M*wz)];

etot[x_, y_, as_] := 
 2 az[h, M, wz]^2/x^2 + az[h, M, wz]^2/y^2 + 
  x^2/(8 l[wx, wz]^2 az[h, M, wz]^2) + y^2/(16 az[h, M, wz]^2) + 
  8 n*az[h, M, wz]^3/(Sqrt[2 Pi] x^2 y) (as/az[h, M, wz] - 
     add *f[x/y]/az[h, M, wz]) + 
  128 n^(3/2) yqf[h, M, as, 
     add]/(25 Sqrt[5] Pi^(9/4) x^3 y^(3/2) h*wz)

a0 = 5.29*10^(-11);
wx = 2 Pi*45;
wz = 2 Pi*133;
h = 1.054*10^(-34);
M = 2.72*10^(-23);
add = 130*a0;
n = 15*10^3;

Clear[x, y, as]
ddata = Table[{x, y, as} /. 
   Last[NMinimize[{etot[x, y, as], x > 0 && y > 0}, {x, y}]], {as, 
   40*a0, 90*a0, 10*a0}]

mminima = Re[etot[#[[1]], #[[2]], #[[3]]]] & /@ ddata

uvalues = Table[as, {as, 40*a0, 90*a0, 10*a0}];
ListPlot[Transpose[{uvalues/a0, mminima}], Joined -> True, 
 Frame -> True]

Here, $M$ is Dy-164 mass and $h$ is the $\hbar$.

$\endgroup$
  • $\begingroup$ It seems to be working. What is the problem that you observed? $\endgroup$ – MelaGo Oct 8 at 0:49
  • $\begingroup$ @MelaGo My plot has not negative values for instance. $\endgroup$ – Dinesh Shankar Oct 8 at 2:15
  • $\begingroup$ Are you sure M shouldn't be 2.72*10^-25? $\endgroup$ – MelaGo Oct 8 at 3:50
1
$\begingroup$

I think the error was the value for M. (Also cleaned up the code a bit.)

ClearAll["Global`*"]

a0 = 5.29*10^(-11);
wx = 2 Pi*45;
wz = 2 Pi*133;
h = 1.054*10^(-34);
(*  M=2.72*10^(-23);  *)
M = 163.9 1.66 10^-27;
add = 130*a0;
n = 15*10^3;
az = Sqrt[h/(M*wz)];
l = wz/wx;

f[X_] := (1 + 2 X^2)/(1 - X^2) - 3 X^2 ArcTanh[Sqrt[1 - X^2]]/(1 - X^2)^(3/2);

yqf[as_] := 32*4 Pi*h^2*as/(3 M) Sqrt[as^3/Pi] (1 + 3 add^2/(2 as^2));

etot[x_, y_, as_] := 2 az^2/x^2 + az^2/y^2 + x^2/(8 l^2 az^2) + y^2/(16 az^2) + 
  8 n*az^3/(Sqrt[2 Pi] x^2 y) (as/az - add*f[x/y]/az) + 
  128 n^(3/2) yqf[as]/(25 Sqrt[5] Pi^(9/4) x^3 y^(3/2) h*wz)

ddata = Table[
  minsol = NMinimize[{etot[x, y, as ], x > 0 && y > 0}, {x, y}];
  Re[{x, y, as/a0 , minsol[[1]]} /. Last[minsol]], {as, 40 a0, 90 a0, a0}]

 (* {{2.54948*10^-7, 5.21789*10^-6, 40., -87.6828}, 
    {2.61239*10^-7, 5.25204*10^-6, 41., -81.3895},
    {2.67699*10^-7, 5.28499*10^-6, 42., -75.5096}, 
    {2.74338*10^-7, 5.31675*10^-6, 43., -70.0129}, etc *)

ListPlot[ddata[[All, {3, 4}]], Joined -> True, PlotRange -> {Automatic, 10}]

enter image description here

This curve looks funny, but it does resemble the Max of the blue and black solid curves in Fig 2 in the paper.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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