Working with the Saha equation in Mathematica

  • Mathematica
  • Thread starter dthomas
  • Start date
  • Tags
    Mathematica
In summary, you are trying to calculate the ionization curve for hydrogen at a given matter density. You use the Saha equation and put the concentrations in terms of reaction coordinates. The equation looks something like this: Log[x[t]^2/(1 - x[t])] = 1/(k t)E - 3/2 Log[-(1/(k t))*E] + 3/2 Log[(-2pimE)/(n^(2/3) h^2)] where k is the Boltzmann constant, h is Planck's constant, and E is a given energy value. t is Temperature. The goal is to get Mathematica to plot the
  • #1
dthomas
4
0
I am trying to calculate the ionization curve for hydrogen at a given matter density (in this case it happens to be 55 atoms/cm^3)

From the Saha equation and putting the concentrations in terms of reaction coordinates. It looks something like this:

Log[x[t]^2/(1 - x[t])] = 1/(k t)E - 3/2 Log[-(1/(k t))E] + 3/2 Log[(-2pimE)/(n^(2/3) h^2)]

k is the Boltzmann constant, h is Planck's constant, and E is a given energy value
t is Temperature

Basically I need to get Mathematica to plot the reaction coordinate (x) vs. t

If anyone has some advice as to how to do this, it would be greatly appreciated
 
Physics news on Phys.org
  • #2
I don't know what version you have, older versions use ImplicitPlot, newer versions use Contour plot.

"E" is reserved by Mathematica for Euler's constant and it will be very upset if you think that is a variable. Pi is reserved by Mathematica for the constant pi. Variables being multiplied either need to be separated by a space or by "*". Plotting versus x[t] will likely be a source of unending grief for you. Given all that, try something like this

In[1]:= <<Miscellaneous`PhysicalConstants`

In[2]:= k=BoltzmannConstant;
h=PlanckConstant;
e=something;
m=something;
n=something;

and either

In[3]:= <<Graphics`ImplicitPlot`

In[4]:= ImplicitPlot[Log[x^2/(1-x)]==1/(k*t)*e-3/2 Log[-(1/(k*t))*e]+3/2 Log[(-2Pi*m*e)/(n^(2/3)*h^2)],{t,something,something}]

or

In[3]:= ContourPlot[Log[x^2/(1-x)]==1/(k*t)*e-3/2 Log[-(1/(k*t))*e]+3/2 Log[(-2*Pi*m*e)/(n^(2/3)*h^2)],{x,something,something},{t,something,something}]

where you fill in all the "something"s with appropriate values and then let us know if it worked or where it failed and what kind of error messages it reported to you.

Read the help pages or Google for

Mathematica TheMathematicaReservedWordYouAreLookingFor

to find the documentation with examples and explanations.
 
  • #3
Thanks so much for the reply it helped quite a bit. I had to tweak it a bit, so here is what I used as my code. (I am in Mathematica 7 by the way)

Clear["Global`*"]
<< PhysicalConstants`;
k = 8.617343*10^-5;(* Boltzmann constant measured in eV/K *)
h = 4.13566733*10^-15; (* Planck's constant measured in
eV s *)
energy = -13.6;
m = 1/1836;(* electron mass * proton mass/ hydrogen mass
measured in proton masses *)
n = 10^-22/1.0794;(* number density (mol H)/cm^3 *)

ContourPlot[
Log[x^2/(1 - x)] ==
1/(k*t)*energy - 3/2 Log[-(1/(k*t))*energy] +
3/2 Log[(-2*\[Pi]*m*energy)/(n^(2/3)*h^2)], {t, 1000, 4000}, {x,
0, -100}]

I am trying to verify a plot given the textbook for my course, but this gave a plot that doesn't match. I am trying to figure out how to upload the plot to this post. I couldn't use the physical constants as the equation needed everything in terms of proton masses and electron volts

http://einstein.drexel.edu/~bob/Thermodynamics/saha.pdf

I am trying to show that the curve holds true for a different matter density (assuming all matter at that density is hydrogen). The curve is on page 4 of the PDF I linked to this post. I hope that this is specific enough, thank you so much for your help
 
Last edited:
  • #4
If you are entering your own constants you can eliminate the
<<PhysicalConstants'

To try to track down the discrepancy I use your constants and then

In[6]:= Simplify[1/(k*t)*energy-3/2 Log[-(1/(k*t))*energy]+3/2 Log[(-2Pi*m*energy)/(n^(2/3)*h^2)]]

Out[6]= 127.536 - 157821./t - 1.5 Log[1/t]

Now I want to see a table of {t,x} to get some idea what is going on

In[7]:= Table[{t, x /.NSolve[Log[x^2/(1-x)] == 127.536-157821./t-1.5` Log[1/t], x]}, {t, 1000, 4000, 100}]

Note: some long strings of trailing digits snipped when put into this post.

Out[7]= {{1000, {-0.0000471618, 0.0000471596}},
{1100, {-0.0683075, 0.06394}},
{1200, {-776.204, 0.998713}},
{1300, {-2.16375*^7}},
{1400, {-1.41081*^11}},
{1500, x},
{1600, x},
{1700, {-8.23303*^19}},
...
{4000, {-4.528245*^43}}}

That isn't going to look anything like what your curves in the pdf look like. I get a stream of warnings about 1/0. Some t have two solutions, some have none, as the t gets close to 4000 the x gets close to minus infinity, not 1.

So something seems to be seriously wrong.

From this can you start to track down where the first layer of errors might be? Maybe I've made some typos you can catch. Maybe you've made some you can catch.

Perhaps once those are found and stripped away then maybe we can start to see what direction to go to be fix the next layer and start in the direction of being able to replicate the example in the text.
 
Last edited:
  • #5
ContourPlot[
Log[x^2/(1 - x)] ==
1/(k*t)*energy - 3/2 Log[-(1/(k*t))*energy] +
3/2 Log[(-2*Pi*m*energy)/(n^(2/3)*h^2)], {t, 1, 4000}, {x, 0, 1}]

I checked my constants and starting equation again this gave me the correct shape curve. I don't know why the code didn't work the first time, but it might have had to do with not clearing whatever was in my program from before this code.

Thank you so much for the help and I appreciate the responses
 

1. What is the Saha equation and why is it important?

The Saha equation is a mathematical expression that describes the ionization state of a gas at a given temperature and pressure. It is important because it allows us to understand and predict the behavior of ionized gases, which are essential components of many natural and technological processes such as stellar atmospheres and plasma physics.

2. How can I use Mathematica to work with the Saha equation?

Mathematica is a powerful computational software that includes built-in functions for solving complex equations such as the Saha equation. You can use the built-in Solve or NSolve functions to find numerical or analytical solutions, or use the Manipulate function to visualize the behavior of the equation with different parameters.

3. Can I customize the Saha equation in Mathematica for specific applications?

Yes, you can customize the Saha equation by defining your own functions and variables in Mathematica. You can also use the built-in functions to add constraints or assumptions to the equation, or use external data to incorporate real-world conditions into your calculations.

4. Are there any limitations to using Mathematica for the Saha equation?

While Mathematica is a powerful tool for working with the Saha equation, it does have some limitations. For example, it may not be able to handle extremely large or complex equations, and it may not have built-in functions for specialized applications of the Saha equation. In these cases, you may need to use alternative software or develop your own code.

5. Are there any resources available for learning how to work with the Saha equation in Mathematica?

Yes, there are many online resources such as tutorials, forums, and documentation provided by Mathematica and its community that can help you learn how to use the software for the Saha equation. Additionally, there are books and online courses available for more in-depth learning and practice.

Similar threads

  • MATLAB, Maple, Mathematica, LaTeX
Replies
2
Views
201
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
4
Views
1K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
4
Views
2K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
9
Views
3K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
12
Views
4K
Replies
0
Views
194
Replies
1
Views
692
  • MATLAB, Maple, Mathematica, LaTeX
Replies
18
Views
3K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
2K
Back
Top