Working with the Saha equation in Mathematica

[dthomas]

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

 

[Bill Simpson]

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.

 

[dthomas]

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

 

[Bill Simpson]

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.

 

[dthomas]

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