I have some trouble understanding the RegionPlot command. I have defined the following input


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

atomNumber[R_, d_, n_] := 2*n*volume[R, 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_] := 
 UnitConvert[\[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_] := 

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

with units

Subscript[m, 0] = Quantity[1, "AtomicMassUnit"]
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"]
\[Omega]Paper = Quantity[40, "Terahertz"]

When I evaluate my functions separately, everything works as intended:

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

coherence[\[Omega]Paper, d, R, rc, \[Lambda], n, T, m] < 1
Out: True

I want to evaluate the above inequality for a range of values for $(\lambda, r_c)$. So: Define an input range for $\lambda$ and for $r_c$ $\Rightarrow$ evaluate coherence(input) < 1 $\Rightarrow$ plot all points $(\lambda, r_c)$ for which this is true.

For that, I use the following RegionPlot command:

 Evaluate[coherence[\[Omega]Paper Quantity[40, "Terahertz"], 
    d Quantity[1/4, "Millimeters"], R Quantity[18/5, "Micrometers"], 
    rc, \[Lambda], 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"]}, {\[Lambda], 
  Quantity[10^-10, ("Seconds")^-1], Quantity[10^-1, ("Seconds")^-1]}]

However, the code doesn't even run and the notebook stops after a while. What is the issue with this code?

My thoughts:

  • the numbers are too small (I can't change the numbers -> are there alternative programs to plot?)
  • there is some issue with my Mathematica version strong text
  • 1
    $\begingroup$ On my machine, v12.2Win7-x64, the RegionPlot command results in a kernel crash. $\endgroup$
    – Syed
    Commented Apr 25, 2022 at 7:28
  • $\begingroup$ Ok but why? As outlined, when I evaluate the functions separately, everything works as intended. What other sites can I consult? Are there alternatives? $\endgroup$ Commented Apr 25, 2022 at 7:38
  • $\begingroup$ @Syed Put differently: Is there a way to define only certain tuples as input? So instead of plotting the full range, I only want to plot 10 values of $\lambda$: $10^{-n}$ for $n \in \{0, 1, ..., 10\}$ $\endgroup$ Commented Apr 25, 2022 at 8:14

1 Answer 1


This is a workaround that at least produces output. It produces a few error messages and took almost 1 minute to calculate, but it may be useful as a starting point. Note that the scale is log-log.

$Version  (*  12.1.0 for Linux x86 (64-bit) (March 14, 2020)  *)

RegionPlot[coherence[ωPaper, d, R,
    10^logrc  Quantity["Meters"],
    10^logλ  Quantity["Seconds"^-1],
    n, T , m ] < 1,
  {logrc, -9, -1},
  {logλ, -10, -1}] // Quiet

enter image description here

All error messages have been suppressed. They are underflows caused by evaluating Exp with a large negative argument. This is considered benign, since the exponential is subtracted from 1.

Another way to visualize the region by plotting the coherence as a 3D surface. An array of about 1000 data points will produce a sufficiently smooth surface. The data points can be generated and plotted like this

data = Flatten[
    Table[{logrc, log\[Lambda], coherence[\[Omega]Paper, d, R,
       10^logrc  Quantity["Meters"],
       10^log\[Lambda]  Quantity["Seconds"^-1],
       n, T , m ]},
     {logrc, -9, -1, .2},
     {log\[Lambda], -6, -1, .2}], 1]; // Quiet

ListPlot3D[data, PlotRange -> {All, All, {0, 1}}]

enter image description here

This plot has been cut off at z=1 in order to visualize the boundary for which coherence == 1.

The logarithmic scales in the above plots may be misleading. A linear plot, such as the one below, shows that the surface is much steeper for small values of rc than it is for the larger values. The following plot is coherence vs rc for $\lambda = 0.00002/\text{sec}$

data2 = Table[{10^logrc, coherence[\[Omega]Paper, d, R,
      10^logrc  Quantity["Meters"],
      2*^-5 /Quantity["Seconds"],
      n, T , m ]},
    {logrc, -8, -3, 0.025}]; // Quiet
ListLinePlot[data2, FrameLabel -> {"rc"}]

enter image description here

  • $\begingroup$ Thanks. Upon reflection, I think I don't need to plot the full range of points to get some output. Is there a way to A) define e.g. 100 input tuples and B) evaluate the inequality for those 100 input tuples? $\endgroup$ Commented Apr 25, 2022 at 8:19
  • 1
    $\begingroup$ @LionCereals When I get a blank plot I will often change the plot command to Table so I can see the what is really going on. I did that on this one, to sort the units problem. I'll add an example of using Table to my answer. Thanks for accepting my answer, but accepting my workaround may discourage another who has the real answer. Feel free to un-accept this workaround and wait to see if you get some really good answers. $\endgroup$
    – LouisB
    Commented Apr 25, 2022 at 8:49

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.