Coding The Cosmos

Hi! Welcome to my website. It's currently under construction, so it doesn't look like much.

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.

Useful Links