Programs
This is a collection of Python programs I've written, mainly about astronomy. For each one,
I have some information about the program, the mathematics behind it, and a link to the
code behind it on GitHub. The code is written in Python 3.x and is designed to be run
mainly via a command line script, but it can of course be run in an IDE, preferably one
with certain libraries built-in. Numpy, matplotlib, and SciPy are some of the more common
ones I use.
Program |
Description |
Popgen |
Popgen is the first incarnation of a stellar modeling program. It uses various
initial mass
functions (IMFs) to model the distribution of stellar masses in a population
of a size of your choice. Currently, the traditional Salpeter IMF and the newer
Kroupa IMF are supported, and the user can choose between the two.
|
Program 2 |
[Placeholder] |
Program 3 |
[Placeholder] |
Popgen
Popgen is a Python command line script designed to create stellar populations
through randomization guided by an
initial mass function (IMF),
a probability density function. At the moment, the program includes two IMFs:
the Salpeter IMF (Salpeter
(1955)) and the Kroupa IMF (Kroupa
(2002)). They are given as follows:
$$\xi(m)\Delta m=\xi_0\left(\frac{m}{M_{\odot}}\right)^{-2.35}\left(\frac{\Delta m}{M_{\odot}}\right)\tag{Salpeter IMF}$$
$$\xi(m)\Delta m=\xi_0\left(\frac{m}{M_{\odot}}\right)^{-\alpha}\left(\frac{\Delta m}{M_{\odot}}\right),\quad\alpha=\begin{cases}0.3,\quad0.01\leq m/M_{\odot}<0.08\\1.3,\quad0.08\leq m/M_\odot<0.50\\2.3,\quad0.50\leq m/M_{\odot}<100\\\end{cases}\tag{Kroupa IMF}$$
The total stars of masses between $m_a$ and $m_b$ is
$$N(m_a,m_b)=\xi_0\int_{m_a}^{m_b}\xi(m)\mathrm{d}m$$
for some to-be-determined normalization constant $\xi_0$.
The program requires the user to import five modules:
matplotlib,
NumPy, SymPy,
random
, and
argparse
:
import matplotlib.pyplot as plt
import numpy as np
from sympy import integrate
from sympy.solvers import solve
from sympy import Symbol
import random
import argparse
There are three parameters you can specify, besides the IMF of your choice. They are
the minimum mass of a star (mmin
), the maximum mass of a star
(mmax
), and the number of stars in the population (N
). The
standard choices I use are $m_{\text{min}}=0.08M_{\odot}$, $m_{\text{max}}=100M_{\odot}$,
and $N=100$. We also have to create the strings R
and
S
, which will denote the number and mass of the stars.
mmin = 0.08
mmax = 100
N = 100
R = []
S = []
We next have to set up the argument parser. It is currently set up to run the
Salpeter IMF or Kroupa IMF, which have shortcuts of -sp
and
-kr
, respectively. This is a simple use of argparse
.
parser = argparse.ArgumentParser()
parser.add_argument("-sp", "--salpeter", action="store_true", help="Salpeter IMF")
parser.add_argument("-kr", "--kroupa", action="store_true", help="Kroupa IMF")
args = parser.parse_args()
Now, we move on to defining the IMFs. An IMF is generally written in the form
While it might seem simpler to define each of the two IMFs in the form $\xi(m)\Delta m$,
and then write just one integrator, this is not possible, because the Kroupa IMF, like
many that are accurate in the low-mass regime, is a piecewise function, and
certain integrators fail in the case of discontinuous functions. We therefore must
define the Kroupa IMF as three separate functions, and integrate them each
individually. The relevant code is
def salpeter(m):
return m**(-2.35)
def salpeterint(a, b):
m = Symbol('m')
result = integrate(salpeter(m), (m, a, b))
return result
def kroupa1(m):
return m**(-0.3)
def kroupa2(m):
return m**(-1.3)
def kroupa3(m):
return m**(-2.3)
def kroupa1int(a, b):
m = Symbol('m')
result = integrate(kroupa1(m), (m, a, b))
return result
def kroupa2int(a, b):
m = Symbol('m')
result = integrate(kroupa2(m), (m, a, b))
return result
def kroupa3int(a, b):
m = Symbol('m')
result = integrate(kroupa3(m), (m, a, b))
return result
We also want to normalize these functions, because it is clear that in this form,
integrating any of them from $m_{\text{min}}$ to $m_{\text{max}}$ and multiplying
that by $N$ will not produce the total number of stars in the population, as it
should. Therefore, we define the following normalization constants:
J = salpeterint(mmin, mmax)
K = (J)**(-1)
a = kroupa1int(mmin, 0.08)
b = kroupa2int(0.08, 0.5)
c = kroupa3int(0.5, mmax)
E = a + b + c
F = (E)**(-1)
A = N*F
B = N*F
C = N*F
I've made A
, B
and C
separate variables for clarity,
even though they're all the same.
Finally, we are ready to define how the IMFs are sliced up to produce a mass for each
star. For the Salpeter IMF, this is relatively symbol, using SymPy. Be careful to use
Symbol('x')
before solving.
def salpsolve(i):
x = Symbol('x')
au = solve(N * K * salpeterint(mmin, x) - i, x)
al = solve(N * K * salpeterint(mmin, x) - (i - 1), x)
p = random.uniform(al[0], au[0])
R.append(i)
S.append(p)
This code lets us create a "bin" within which each star resides. au
and al
are the upper and lower bounds, respectively, but do not
determine the mass of the star. That's where random.uniform
comes
in, to introduce an element of randomness and ensure that each new run does
not produce the exact same group of stars.
We have to be slightly more careful with the Kroupa IMF. Rather than naively
creating one solve
pair of equations to solve, we need three, to cover
all three parts of the piecewise IMF. Therefore, we have
def kroupsolve(i):
x = Symbol('x')
if i <= a * A:
au = solve(A * kroupa1int(mmin, x) - i, x)
al = solve(A * kroupa1int(mmin, x) - (i - 1), x)
elif a * A < i and i <= b * B:
au = solve(A * a + B * kroupa2int(0.08, x) - i, x)
al = solve(A * a + B * kroupa2int(0.08, x) - (i - 1), x)
elif b * B < i:
au = solve(A * a + B * b + C * kroupa3int(0.5, x) - i, x)
al = solve(A * a + B * b + C * kroupa3int(0.5, x) - (i - 1), x)
p = random.uniform(al[0], au[0])
R.append(i)
S.append(p)
At this point, all of the heavy lifting has been done. We just need to determine which
IMF is being used, and then output the data in some form. The first is fairly simple,
going back to our argument parser:
if args.salpeter:
for i in range(1, N + 1):
salpsolve(i)
elif args.kroupa:
for i in range(1, N + 1):
kroupsolve(i)
Choosing how to output the data is a bit more complicated. You can likely modify this
section to suit your needs - it's rudimentary at the moment - but it currently creates
a plot of the number of stars against the mass of the $i$th star. The code is fairly
straightforward:
plt.semilogx(S, R, 'bx')
plt.xlabel('$\log m\:(M_{\odot})$')
plt.ylabel('$N(m\leq i)$')
plt.xlim(-1, 1.1*mmax)
plt.ylim(-1, 1.1*N)
plt.show()
This simply uses S
and R
as the $x$- and $y$- axes, respectively.
Logarithmic axes are necessary in order to resolve distinct stars at the lower end of the
mass spectrum.
The second output is as a .txt
file, unimaginatively titled "Data". It has
crude columns, with the mass of each star in solar masses given next to its number:
file = open('Data.txt', 'w')
for index in range(len(R)):
file.write("Number" + " " + "Mass" + "\n")
file.write(" " + str(R[index]) + " " + str(S[index]) + "\n")
file.close()
I do have plans to expand Popgen into something bigger. Obviously, I can come up with
better ways of graphically exploring the population, in the short term, but I can also
use it as the input into stellar population models. There are some awesome isochrone
and evolutionary track generators online, and it might be possible to use one to better
understand how a slice of a population will live and die. This version of Popgen, though,
can stand on its own.
Program 2
Program 3
All content on this page is licensed under the
Creative Commons Attribution
4.0 International License.