skirt-stab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "700ab242",
   "metadata": {},
   "source": [
    "### FSPS implementation in SKIRT: stab file\n",
    "This script will take FSPS (https://github.com/cconroy20/fsps) data output and make an interpolation to increase the resolution ($R= \\frac{ \\Delta \\lambda}{\\lambda}$) as stated in the document by Baes, Camps, Kapoor & Rivas-Sánchez 2022.\n",
    "\n",
    "J. Fritz & M Rivas-Sánchez:\n",
    "\n",
    "\"When we set up a library, one sensible option could be to assume a fixed spectral resolution for the lines,\n",
    "meaning that we choose a value R, and determine the width $\\Delta \\lambda_{k}$ of each line according to\n",
    "\n",
    "$\\Delta \\lambda_{k} = \\frac{2\\lambda_{k}}{R}$ or $R = \\frac{2\\lambda_{k}}{\\Delta \\lambda_{k}} $.\"\n",
    "\n",
    "$FWHM = 2 \\sqrt{2ln(2)} \\sigma$\n",
    "\n",
    "\"If we know that the area under the Gaussian curve is ∼ 99.993% when we reach 4$\\sigma$ from the center, we can safely truncate the Gaussian profile at 4$\\sigma$. For a given intrinsic wavelength resolution R, we can then determine the standard deviation of the k’th line $\\sigma_{k}$ as:\n",
    "\n",
    "$\\sigma_{k} = \\frac{\\Delta \\lambda_{k}}{8} = \\frac{\\lambda_{k}}{4R}$.\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "70eca529",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Importing libraries\n",
    "import fsps\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "import math\n",
    "import pts.storedtable.io as stab\n",
    "from astropy.visualization import astropy_mpl_style\n",
    "import astropy.io.ascii as ascii\n",
    "from __future__ import print_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "aa47b44c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Technical specifications\n",
    "nint = 37 #number of new points per emission lines\n",
    "fwhm2sig = 2.*np.sqrt(np.log(2)) #relationship between fwhm and sigma\n",
    "#fsps resolution (fwhm)\n",
    "fwhm_hires, fwhm_lores = 2.3, 30.\n",
    "#solving for sigma  \n",
    "sig_hires, sig_lores = fwhm_hires/fwhm2sig, fwhm_lores/fwhm2sig\n",
    "#delta resolution\n",
    "Delta_hires, Delta_lores = 4.*sig_hires, 4.*sig_lores\n",
    "#FSPS [3600,7400] Ang resolution FWHM=2.3 Ang\n",
    "l1, l2 = 3600., 7400.\n",
    "\n",
    "#Reading files\n",
    "llines = ascii.read('/Users/.../fsps_lumlines.txt')\n",
    "nlines = len(llines)\n",
    "lc = np.sort(llines) #increasing wavelength lines \n",
    "\n",
    "spectrum = ascii.read('/Users/.../data_fsps.txt')\n",
    "nwav = len(spectrum) #number of wavelengths\n",
    "wavel = spectrum[0][:] #wavelength in Ang"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c029b673",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Creating new array for interpolated wavelengths\n",
    "int_lin_array = []\n",
    "new_array = []\n",
    "llim = 0\n",
    "\n",
    "for _i in range(nlines):\n",
    "    if lc[_i][0] >= l1 and lc[_i][0] <= l2: #Miles resolution\n",
    "        delta = Delta_hires\n",
    "    else:\n",
    "        delta = Delta_lores\n",
    "    li = lc[_i][0] - delta #lower\n",
    "    lu = lc[_i][0] + delta #upper\n",
    "    dum_wave = wavel[np.where(wavel < li)]\n",
    "    wave = wavel[np.where(dum_wave > llim)]\n",
    "    llim = lu\n",
    "    new_array = np.append(new_array, wave) #\n",
    "    dl = (lu-li)/float(nint-1)\n",
    "    dum_arr = [(li + dl*n) for n in range(nint)] #new wavelength points\n",
    "    new_array = np.append(new_array, dum_arr)\n",
    "\n",
    "new_wavel = np.sort(new_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "667361f9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initializing FSPS...\n",
      "(b'mist', b'miles', b'DL07  ')\n",
      "Determining grids...\n",
      "Calculating spectra for logZ=-2.5 and logU=-4...\n",
      "Calculating spectra for logZ=-2.5 and logU=-3.5...\n",
      "Calculating spectra for logZ=-2.5 and logU=-3...\n",
      "Calculating spectra for logZ=-2.5 and logU=-2.5...\n",
      "Calculating spectra for logZ=-2.5 and logU=-2...\n",
      "Calculating spectra for logZ=-2.5 and logU=-1.5...\n",
      "Calculating spectra for logZ=-2.5 and logU=-1...\n",
      "Calculating spectra for logZ=-2 and logU=-4...\n",
      "Calculating spectra for logZ=-2 and logU=-3.5...\n",
      "Calculating spectra for logZ=-2 and logU=-3...\n",
      "Calculating spectra for logZ=-2 and logU=-2.5...\n",
      "Calculating spectra for logZ=-2 and logU=-2...\n",
      "Calculating spectra for logZ=-2 and logU=-1.5...\n",
      "Calculating spectra for logZ=-2 and logU=-1...\n",
      "Calculating spectra for logZ=-1.5 and logU=-4...\n",
      "Calculating spectra for logZ=-1.5 and logU=-3.5...\n",
      "Calculating spectra for logZ=-1.5 and logU=-3...\n",
      "Calculating spectra for logZ=-1.5 and logU=-2.5...\n",
      "Calculating spectra for logZ=-1.5 and logU=-2...\n",
      "Calculating spectra for logZ=-1.5 and logU=-1.5...\n",
      "Calculating spectra for logZ=-1.5 and logU=-1...\n",
      "Calculating spectra for logZ=-1 and logU=-4...\n",
      "Calculating spectra for logZ=-1 and logU=-3.5...\n",
      "Calculating spectra for logZ=-1 and logU=-3...\n",
      "Calculating spectra for logZ=-1 and logU=-2.5...\n",
      "Calculating spectra for logZ=-1 and logU=-2...\n",
      "Calculating spectra for logZ=-1 and logU=-1.5...\n",
      "Calculating spectra for logZ=-1 and logU=-1...\n",
      "Calculating spectra for logZ=-0.6 and logU=-4...\n",
      "Calculating spectra for logZ=-0.6 and logU=-3.5...\n",
      "Calculating spectra for logZ=-0.6 and logU=-3...\n",
      "Calculating spectra for logZ=-0.6 and logU=-2.5...\n",
      "Calculating spectra for logZ=-0.6 and logU=-2...\n",
      "Calculating spectra for logZ=-0.6 and logU=-1.5...\n",
      "Calculating spectra for logZ=-0.6 and logU=-1...\n",
      "Calculating spectra for logZ=-0.4 and logU=-4...\n",
      "Calculating spectra for logZ=-0.4 and logU=-3.5...\n",
      "Calculating spectra for logZ=-0.4 and logU=-3...\n",
      "Calculating spectra for logZ=-0.4 and logU=-2.5...\n",
      "Calculating spectra for logZ=-0.4 and logU=-2...\n",
      "Calculating spectra for logZ=-0.4 and logU=-1.5...\n",
      "Calculating spectra for logZ=-0.4 and logU=-1...\n",
      "Calculating spectra for logZ=-0.3 and logU=-4...\n",
      "Calculating spectra for logZ=-0.3 and logU=-3.5...\n",
      "Calculating spectra for logZ=-0.3 and logU=-3...\n",
      "Calculating spectra for logZ=-0.3 and logU=-2.5...\n",
      "Calculating spectra for logZ=-0.3 and logU=-2...\n",
      "Calculating spectra for logZ=-0.3 and logU=-1.5...\n",
      "Calculating spectra for logZ=-0.3 and logU=-1...\n",
      "Calculating spectra for logZ=-0.2 and logU=-4...\n",
      "Calculating spectra for logZ=-0.2 and logU=-3.5...\n",
      "Calculating spectra for logZ=-0.2 and logU=-3...\n",
      "Calculating spectra for logZ=-0.2 and logU=-2.5...\n",
      "Calculating spectra for logZ=-0.2 and logU=-2...\n",
      "Calculating spectra for logZ=-0.2 and logU=-1.5...\n",
      "Calculating spectra for logZ=-0.2 and logU=-1...\n",
      "Calculating spectra for logZ=-0.1 and logU=-4...\n",
      "Calculating spectra for logZ=-0.1 and logU=-3.5...\n",
      "Calculating spectra for logZ=-0.1 and logU=-3...\n",
      "Calculating spectra for logZ=-0.1 and logU=-2.5...\n",
      "Calculating spectra for logZ=-0.1 and logU=-2...\n",
      "Calculating spectra for logZ=-0.1 and logU=-1.5...\n",
      "Calculating spectra for logZ=-0.1 and logU=-1...\n",
      "Calculating spectra for logZ=0 and logU=-4...\n",
      "Calculating spectra for logZ=0 and logU=-3.5...\n",
      "Calculating spectra for logZ=0 and logU=-3...\n",
      "Calculating spectra for logZ=0 and logU=-2.5...\n",
      "Calculating spectra for logZ=0 and logU=-2...\n",
      "Calculating spectra for logZ=0 and logU=-1.5...\n",
      "Calculating spectra for logZ=0 and logU=-1...\n",
      "Calculating spectra for logZ=0.1 and logU=-4...\n",
      "Calculating spectra for logZ=0.1 and logU=-3.5...\n",
      "Calculating spectra for logZ=0.1 and logU=-3...\n",
      "Calculating spectra for logZ=0.1 and logU=-2.5...\n",
      "Calculating spectra for logZ=0.1 and logU=-2...\n",
      "Calculating spectra for logZ=0.1 and logU=-1.5...\n",
      "Calculating spectra for logZ=0.1 and logU=-1...\n",
      "Calculating spectra for logZ=0.2 and logU=-4...\n",
      "Calculating spectra for logZ=0.2 and logU=-3.5...\n",
      "Calculating spectra for logZ=0.2 and logU=-3...\n",
      "Calculating spectra for logZ=0.2 and logU=-2.5...\n",
      "Calculating spectra for logZ=0.2 and logU=-2...\n",
      "Calculating spectra for logZ=0.2 and logU=-1.5...\n",
      "Calculating spectra for logZ=0.2 and logU=-1...\n",
      "Calculating spectra for logZ=0.5 and logU=-4...\n",
      "Calculating spectra for logZ=0.5 and logU=-3.5...\n",
      "Calculating spectra for logZ=0.5 and logU=-3...\n",
      "Calculating spectra for logZ=0.5 and logU=-2.5...\n",
      "Calculating spectra for logZ=0.5 and logU=-2...\n",
      "Calculating spectra for logZ=0.5 and logU=-1.5...\n",
      "Calculating spectra for logZ=0.5 and logU=-1...\n",
      "Writing stored table file FSPS-emissionlines.stab...\n",
      "\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "This fragment was taken from:\n",
    "https://skirt.ugent.be/pts9/namespacepts_1_1storedtable_1_1do_1_1fsps__with__lines__to__stored__table.html\n",
    "\n",
    "# *****************************************************************\n",
    "# **       PTS -- Python Toolkit for working with SKIRT          **\n",
    "# **       © Astronomical Observatory, Ghent University          **\n",
    "# *****************************************************************\n",
    "\"\"\"\n",
    "# define the script body inside a function to avoid execution at the top level\n",
    "def do() -> \"Create FSPS-generated SED family including emission lines as stored table\":\n",
    "    # the file path for the .stab output file (here placed in the current directory)\n",
    "    stabPath = \"FSPS-emissionlines.stab\"\n",
    "\n",
    "    # the metallicity and ionization parameter grids\n",
    "    Zsun = 0.0142 #Padova = 0.0190\n",
    "    logZv = [-2.5, -2, -1.5, -1, -0.6, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.5]\n",
    "    logUv = [-4, -3.5, -3, -2.5, -2, -1.5, -1]\n",
    "\n",
    "    # initialize the ssp object (it is advised to do this only once)\n",
    "    print(\"Initializing FSPS...\", flush=True)\n",
    "    sp = fsps.StellarPopulation(compute_vega_mags = False,\n",
    "                                zcontinuous=1,               # Type of interpolation for metallicity\n",
    "                                sfh=0,                       # SSPs\n",
    "                                imf_type=2,                  # Kroupa IMF\n",
    "                                add_dust_emission=False,     # Turn off dust emission\n",
    "                                add_neb_emission=True,       # Turn on nebular emission\n",
    "                                add_neb_continuum=True,\n",
    "                                nebemlineinspec=True,        # Removes lines from spectra\n",
    "                                cloudy_dust=False,           # Exclude dust from nebular emission\n",
    "                                logzsol=-2,                  # Stellar metallicity\n",
    "                                gas_logz=-2,                 # Gas metallicity (must be the same as stellar)\n",
    "                                gas_logu=-4)                 # Ionization parameter\n",
    "    \n",
    "    print (sp.libraries)\n",
    "\n",
    "    # perform a dummy calculation to get the wavelength and age grids\n",
    "    # (if we don't specify the age, a data cube with SSPs for a range of ages is created)\n",
    "    print(\"Determining grids...\", flush=True)\n",
    "    wv, _ = sp.get_spectrum(peraa=True)     # wavelengths in Angstrom\n",
    "    logtv = sp.log_age\n",
    "    tv = 10.**logtv                         # ages in year\n",
    "    Zv = Zsun * 10.**np.array(logZv)        # metallicities from log(Z/Zsun)\n",
    "    Uv = 10.**np.array(logUv)\n",
    "\n",
    "    # allocate hypercube for the luminosities\n",
    "    L = np.zeros((len(new_wavel), len(Zv), len(tv), len(Uv)))\n",
    "\n",
    "    # perform calculations for all metallicities and ionization parameters\n",
    "    for iZ, logZ in enumerate(logZv):\n",
    "        for iU, logU in enumerate(logUv):\n",
    "            print(\"Calculating spectra for logZ={} and logU={}...\".format(logZ, logU), flush=True)\n",
    "\n",
    "            # modify parameters in the sp object\n",
    "            sp.params['logzsol'] = logZ\n",
    "            sp.params['gas_logz'] = logZ\n",
    "            sp.params['gas_logu'] = logU\n",
    "\n",
    "            # get specific luminosities in Lsun/Angstrom and put them into the hypercube\n",
    "            _, final_spectra = sp.get_spectrum(peraa=True)\n",
    "            \n",
    "            Lv = np.zeros((final_spectra.shape[0],new_wavel.shape[0]))\n",
    "\n",
    "            for i in range(Lv.shape[0]):\n",
    "                Lv[i,:] = np.interp(new_wavel, _, final_spectra[i,:])\n",
    "            \n",
    "            L[:, iZ, :, iU] = Lv.T\n",
    "\n",
    "    # convert from input units Lsun/Angstrom to output units W/m\n",
    "    L *= 3.848e26 * 1e10\n",
    "    wv *= 1e-10\n",
    "\n",
    "    # write the stored table\n",
    "    print(\"Writing stored table file {}...\".format(stabPath), flush=True)\n",
    "    stab.writeStoredTable(stabPath,\n",
    "            ['lambda','Z','t','U'], ['m','1','yr','1'], ['log','log','log','log'], [new_wavel,Zv,tv,Uv],\n",
    "            ['Llambda'], ['W/m'], ['log'], [L])\n",
    "    print()\n",
    "    print(\"Done.\", flush=True)\n",
    "\n",
    "\n",
    "# remove the \"#\" in front of the next line to perform the script body\n",
    "do()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "156fc421",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}