main_sequence.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import ascii\n",
    "from scipy import interpolate\n",
    "from scipy import integrate\n",
    "from astropy import constants as constants\n",
    "from astropy import units as u\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ms_table_file = 'ekert_2018_table_7_MS.dat'\n",
    "dat = ascii.read(ms_table_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Spt', 'logTeff', 'B_V', 'U_B', 'BC', 'MV', 'Teff', 'M', 'R', 'logg', 'M_over_L', 'L_over_M']\n"
     ]
    }
   ],
   "source": [
    "# Names of columns in the ascii table\n",
    "print dat.colnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Functions to interpolate Teff, luminosities, as a function of stellar mass\n",
    "Teff_vs_M = interpolate.interp1d(dat['M'], dat['Teff'], \n",
    "                                 fill_value='extrapolate')\n",
    "logTeff_vs_M = interpolate.interp1d(dat['M'], dat['logTeff'], \n",
    "                                 fill_value='extrapolate')\n",
    "L_over_M_vs_M = interpolate.interp1d(dat['M'], dat['L_over_M'], \n",
    "                                 fill_value='extrapolate')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# set up properties of the IMF\n",
    "Mmin = 0.08\n",
    "Mmax = 120.\n",
    "slope = -2.35\n",
    "\n",
    "# function to generate a power law initial mass function initialized\n",
    "# to one solar mass when integrated between Mmin and Mmax\n",
    "#\n",
    "# https://en.wikipedia.org/wiki/Initial_mass_function\n",
    "#\n",
    "def imf_power_law(M, slope=slope, Mmin=Mmin, Mmax=Mmax):\n",
    "    # normalization \n",
    "    norm = 1.0 / ((1.0/(slope+1))*(Mmax**(slope+1) - Mmin**(slope+1)))\n",
    "    imf = norm * M**slope\n",
    "    # remove stars above and below the mass range over which the IMF is normalized\n",
    "    imf[M<Mmin]=0.0\n",
    "    imf[M>Mmax]=0.0\n",
    "    return imf\n",
    "\n",
    "# check normalization\n",
    "#print integrate.quad(imf_power_law,\n",
    "#                     Mmin, Mmax, args=(slope))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set up vector of Masses, Temperatures, bolometric Luminosities, \n",
    "# and IMF weighting\n",
    "#\n",
    "# Units: Solar Masses, Kelvin, Solar Luminosities, inverse solar masses\n",
    "npoints = 200\n",
    "lgM = np.linspace(np.log10(Mmin),np.log10(Mmax),npoints)\n",
    "M = 10**lgM\n",
    "T = Teff_vs_M(M)\n",
    "Lbol = L_over_M_vs_M(M) * M\n",
    "imf = imf_power_law(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Function to return a black body spectrum for wavelength in angstroms\n",
    "# and temperature in Kelvin. Output in cgs units.\n",
    "#\n",
    "# This can be integrated using:\n",
    "# integrate.quad(blackbody_wavelength_angstroms,\n",
    "#                lambda_min, lambda_max,\n",
    "#                args=(temperature))\n",
    "# which returns a 2-element tuple where the first value contains the definite integral\n",
    "\n",
    "def blackbody_wavelength_angstroms(lam, temperature):\n",
    "    _h = constants.h.cgs.value\n",
    "    _c = constants.c.cgs.value\n",
    "    _k_B = constants.k_B.cgs.value\n",
    "    \n",
    "    I = (2.0*_h*_c**2 / lam**5 * (np.exp(_h*_c/(_k_B*temperature*(lam*1.e-8))) \n",
    "                                  - 1.)**-1)\n",
    "\n",
    "    return I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# function to return the log of the main sequence lifetime \n",
    "# as a function of stellar mass\n",
    "\n",
    "def ms_lifetime_log_yrs(M):\n",
    "    lgt = (10.09 - 3.139*np.log10(M)\n",
    "                   + 0.238238*np.log10(M)**2\n",
    "                   + 0.26163378*np.log10(M)**3)\n",
    "    # > 8 Msun\n",
    "    lgt_hi = (9.01 - 1.57*np.log10(M))\n",
    "    i = np.where(M>20)[0]\n",
    "    lgt[i]=lgt_hi[i]\n",
    "    return lgt"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}