Problem Setting

I have defined the following functions.

atomNumber[R_, d_, n_] := 2*n*volume[R, d]

volume[R_, d_ ] := Pi R^2 d

\[Gamma]Fn[x_] := 2/x^2 (1 - Exp[-x^2] (BesselI[0, x^2] + BesselI[1, x^2]))

\[Eta]Eval[\[Lambda]_, d_,  R_, rc_, n_] := 
\[Lambda] *((atomNumber[R, d, n])^2/d^2) \[Gamma]Fn[R/(Sqrt[2] rc)] (1 - Exp[-d^2/(4*rc^2)])

\[Del]Fn[\[Omega]_, m_] := Sqrt[Quantity["ReducedPlanckConstant"]/(\[Omega]*m)]

The inputs to all these functions should eventually be physical units. I then combine these functions towards the following:

coherence[\[Omega]_, d_, R_, rc_, \[Lambda]_, n_, T_, m_] := 
 4*\[Eta]Eval[\[Lambda], d,  R, rc, n]*T*(\[Del]Fn[\[Omega], m])^2

If I provide units and evaluate the functions separately, everything works fine, so there is no issue with any definition here.


I two variables: $r_c$ (rc) and $\lambda$ (\[Lambda]), all other inputs are parameters. I want to find all tuples $(r_c, \lambda)$ for which the following inequality holds:

coherence[\[Omega]_, d_, R_, rc_, \[Lambda]_, n_, T_, m_] <= 1

In particular, I want to use Manipulate to find these tuples for different parameter values, i.e. I want to have slides where I can change $\omega, d, R$ given the right physical units and then get a plot for the tuples.

My solution

I use RegionPlot for the inequality and Manipulate for the parameter slide. That leads to

  coherence[\[Omega]Paper Quantity[40, "Terahertz"], 
    d  Quantity[0.25, "Millimeters"], R Quantity[3.6, "Micrometers"], 
    rc, \[Lambda], n Quantity[176.2 * 10^(27), "Meters"^(-3)], 
    T Quantity[350, "Femtoseconds"], 
    m Quantity[6, "AtomicMassUnit"]] < 1, {rc , 
   Quantity[10^-9, "Meters"], Quantity[10^-1, "Meters"]}, {\[Lambda], 
   Quantity[10^-10, ("Seconds")^-1], 
   Quantity[10^-1, ("Seconds")^-1]}, 
  PlotRange -> All], {\[Omega]Paper, 20, 60}, {d, 0.01, 1}, {R, 0.01, 
  10}, {n, 100, 200}, {T, 300, 400}, {m, 1, 10}]

Here, I specified my variables $r_c$ and $\lambda$ and provided slides for all other parameters. However, I get an empty plot not showing anything.

Can anybody help me with this? I suspect that the problem has something to do with the Quantities package and my misuse thereof.

Additional Input

This is the list of hard-coded parameters I use:

\[Omega] = Quantity[40, "Terahertz"]
R = Quantity[3.6, "Micrometers"]
d = Quantity[0.25, "Millimeters"]
n = Quantity[176.2 * 10^(27), "Meters"^(-3)]
rc = Quantity[10^(-7), "Meters"]
\[Lambda] = Quantity[10^(-17), "Seconds"^(-1)]
T = Quantity[350, "Femtoseconds"]
m = Quantity[6, "AtomicMassUnit"]

This shows that coherence works:

testCoherence = coherence[\[Omega], d, R, rc, \[Lambda], n, T, m]
Out: 2.28009*10^-15

This is an example of RegionPlot without Manipulate, which also shows an empty plot:

 Evaluate[coherence[\[Omega]Paper, d, R, rc, \[Lambda], n, T, m]] < 
  1, {rc , Quantity[10^-9, "Meters"], 
  Quantity[10^-1, "Meters"]}, {\[Lambda], 
  Quantity[10^-10, ("Seconds")^-1], Quantity[10^-1, ("Seconds")^-1]}, 
 PlotRange -> Automatic]
  • $\begingroup$ Can you present an example of a working RegionPlot for your data with hard coded parameters of your choice as a starting point? $\endgroup$
    – Syed
    Mar 26, 2022 at 14:58
  • $\begingroup$ You have not provided the definition of the function volume $\endgroup$
    – Bob Hanlon
    Mar 26, 2022 at 15:13
  • $\begingroup$ @BobHanlon Fixed both issues $\endgroup$ Mar 26, 2022 at 15:18

1 Answer 1


You should also add some precision control.


volume[R_, d_] := Pi R^2 d

atomNumber[R_, d_, n_] := 2*n*volume[R, d]

γFn[x_] := 2/x^2 (1 - Exp[-x^2] (BesselI[0, x^2] + BesselI[1, x^2]))

ηEval[λ_, d_, R_, rc_, n_] := 
 UnitConvert[λ*((atomNumber[R, d, n])^2/d^2) γFn[
    R/(Sqrt[2] rc)] (1 - Exp[-d^2/(4*rc^2)])]

∇Fn[ω_, m_] := 

coherence[ω_, d_, R_, rc_, λ_, n_, T_, m_] := 
 UnitConvert[4*ηEval[λ, d, R, rc, n]*T*(∇Fn[ω, m])^2]


     ωPaper Quantity[40, "Terahertz"],
     d Quantity[1/4, "Millimeters"],
     R Quantity[18/5, "Micrometers"],
     rc, λ,
     n Quantity[1762*10^(26), "Meters"^(-3)],
     T Quantity[350, "Femtoseconds"],
     m Quantity[6, "AtomicMassUnit"]] < 1],
  {rc, Quantity[10^-9, "Meters"], Quantity[10^-1, "Meters"]},
  {λ, Quantity[10^-10, ("Seconds")^-1], 
   Quantity[10^-1, ("Seconds")^-1]},
  PlotPoints -> 75,
  MaxRecursion -> 3,
  PlotRange -> All],
 {{ωPaper, 30}, 20, 60, 0.25, Appearance -> "Labeled"},
 {{d, 0.25}, 0.01, 1, 0.01, Appearance -> "Labeled"},
 {{R, 2}, 0.01, 10, 0.05, Appearance -> "Labeled"},
 {{n, 125}, 100, 200, 1, Appearance -> "Labeled"},
 {{T, 325}, 300, 400, 1, Appearance -> "Labeled"},
 {{m, 3}, 1, 10, 1, Appearance -> "Labeled"},
 SynchronousUpdating -> False,
 TrackedSymbols :> All]
  • $\begingroup$ Thanks. I see that you specified {n, 125}. In what units is now n? Is it like in n Quantity[1762*10^(26), "Meters"^(-3)], i.e. we have $125 * 10^{26} m^{-3}$? $\endgroup$ Mar 26, 2022 at 16:17
  • $\begingroup$ Like all of your control variables, n is a number (dimensionless) which is multiplied by a Quantity inside the arguments of coherence. The arguments to coherence have the units of the associated quantities. You originally defined n as a number in the interval {100, 200}. I gave it an initial value (125) within that interval and specified a step size of 1. $\endgroup$
    – Bob Hanlon
    Mar 26, 2022 at 16:36
  • $\begingroup$ Got it. I now changed {{n, 125}, 100, 200, 1, Appearance -> "Labeled"} to {{n, 176.2 * 10^(27)}, 170 * 10^(27), 180 * 10^(27), 10^26, Appearance -> "Labeled"} and again see no plot. But I will look into it, thanks. $\endgroup$ Mar 26, 2022 at 16:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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