more notebooks

notebooks/BinaryOrbits.py

+# ---
+# jupyter:
+#   jupytext:
+#     formats: ipynb,py:light
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.13.0
+#   kernelspec:
+#     display_name: Python 3 (ipykernel)
+#     language: python
+#     name: python3
+# ---
+
+# # Orbital phase calculation and observability
+#
+# Calculate orbital phase and observability for given binary and observatory.
+
+# ## Parameters
+#
+# Parameters to be set:
+# - binary parameters: from yaml file in ../data/ directory
+# - observatory paramater: from yaml file in ../data/ directory
+# - observation period: below
+
+from astropy.time import Time
+import astropy.units as u
+import matplotlib.pyplot as plt
+import numpy as np
+import yaml
+
+# time range for observations
+obs_range = Time( ['2021-09-01T00:00:00', '2022-03-01T00:00:00'], format='isot', scale='utc')
+
+binary_data = "../data/HESSJ0632p057.yaml"
+binary_data = "../data/LSIp61303.yaml"
+
+
+# ## Utilities
+
+# binary
+def yaml_data(ymfle):
+    "return data from yaml file"
+    try:
+        stream = open(ymfle, 'r')
+    except:
+        print('file not found: %s' % (ymfle) )
+        return None
+    data = list(yaml.load_all(stream, Loader=yaml.FullLoader))
+    return data
+
+
+# # Orbital phase
+
+data = yaml_data(binary_data)
+print(data)
+
+
+def orbital_phase(mjd, mjd0, period):
+    "orbital phase calculation"
+    return (mjd-mjd0)/period - np.int_((mjd-mjd0)/period)
+
+
+time_axis=np.linspace(obs_range[0].mjd,
+                      obs_range[1].mjd,
+                      int(obs_range[1].mjd-obs_range[0].mjd)+1)
+phase_axis=orbital_phase(time_axis,
+                         data[0]['orbit']['mjd0']['val'], 
+                         data[0]['orbit']['period']['val'])
+plt.plot(time_axis,
+         phase_axis)
+plt.xlabel('MJD')
+plt.ylabel('orbital phase')
+plt.grid()
+plt.title("%s (MJD0=%.1f,period=%.1f)" %  (data[0]['name'],
+                                           data[0]['orbit']['mjd0']['val'],
+                                           data[0]['orbit']['period']['val']) )
+plt.show()
+
+# # Observability
+#
+# time vs orbital phase vs election
+
+from astroplan import FixedTarget
+from astroplan import Observer
+
+flwo = Observer(latitude=31.6716989799*u.deg, 
+                longitude=-110.951291195*u.deg,
+                elevation=1.268*u.m, 
+                name='FLWO',
+                timezone="America/Phoenix")
+
+binary_object = FixedTarget.from_name(data[0]['name'])
+
+print(binary_object)
+
+

notebooks/MassiveStars.py

+# ---
+# jupyter:
+#   jupytext:
+#     formats: ipynb,py:light
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.13.0
+#   kernelspec:
+#     display_name: Python 3 (ipykernel)
+#     language: python
+#     name: python3
+# ---
+
+# # Introduction to Be stars
+#
+# - [Chen, A. M. ; Guo, Y. D. ; Yu, Y. W. ; Takata, J., Astronomy & Astrophysics, Volume 652, id.A39, 8 pp. (2021)](https://ui.adsabs.harvard.edu/abs/2021A%26A...652A..39C/abstract)
+# - [Waters et al 1988](https://ui.adsabs.harvard.edu/abs/1988A%26A...198..200W/abstract) (stellar winds)
+#
+# Massloss in massive star occur via intensive outflows:
+# - a high-velocity low-density polar wind
+# - a low-velocity high-densitiy equatorial disk for Oe/Be stars
+#
+# Mass loss rates are of the order of $1-2 \times 10^{-7} M_{\odot}/yr$ (Waters et al 1988)
+
+import astropy.units as u
+import astropy.constants as const
+import math
+import matplotlib.pyplot as plt
+import numpy as np
+
+# ## Parameters used in examples
+#
+# Following Table 1 in [Chen et al 2021](https://ui.adsabs.harvard.edu/abs/2021A%26A...652A..39C/abstract):
+
+# +
+# general parameters
+
+# assumption from Chen et al 2021
+u_i = 0.62
+
+# PSR B1259-63/LS 2883
+object_name = "PSR B1259-63/LS 2883"
+
+orbital_period = 1236.7 * u.d
+
+# Johnston et al 1992; Shannon et al 2014)
+pulsar_period = 47.6e-3 * u.s
+pulsar_luminosity = 8.2e35 * u.erg/u.s
+
+# Negueruela et al 2011
+stellar_type = "O9.5Ve"
+stellar_mass = 30 * u.M_sun
+stellar_radius = 9.2 * u.R_sun
+stellar_T = 3.02e4 * u.K
+stellar_massloss = 2e-8 * u.M_sun/u.yr
+stellar_windvelocity = 3e8 * u.cm/u.s
+
+# LS 5039
+
+#stellar_type = "O"
+#stellar_mass = 22.9 * u.M_sun
+#stellar_radius = 9.3 * u.R_sun
+#stellar_T = 3.90e4 * u.K
+# -
+
+# ## Isotropic Wind
+#
+# Number density of the wind is given (following [Waters et al. 1988](https://ui.adsabs.harvard.edu/abs/1988A%26A...198..200W/abstract)):
+# $$
+#    n_{w,i}(r) = n_{w,0} \left( \frac{r}{r_{\star}} \right)^{-2}
+# $$
+# - $r$: radial distance from center to star
+# - $r_{\star}$ stellar surface
+# - $n_{w,0}$ base density at surface of star
+#
+# The base density depends on the mass loss rate $\dot{M}$ and wind velocity $v_{W}$:
+# $$
+#   n_{w,0} = \dot{M} / 4\pi r_{\star}^2 v_w \mu_i m_p
+# $$
+# - $\mu_i$: mean ion molecular weight
+# - $m_p$: proton mass
+
+n_w_0 = stellar_massloss/(4.*math.pi*stellar_radius**2*stellar_windvelocity*u_i*const.m_p)
+print("Base number densitity (wind): {:.2e}".format(n_w_0.to(u.cm**-3)))
+
+
+# **Note:** this value disagrees with Chen et al 2021, who cite $n_{w,0} = 4.12\times 10^8$ cm$^{-3}$
+
+# +
+def iso_wind_density(r):
+    "isotropic wind densitity"
+    return n_w_0 * (r/stellar_radius)**(-2.)
+
+distance_axis=np.linspace(1., 50., 50 )
+density_axis = iso_wind_density(distance_axis)
+
+plt.gcf().clear()
+plt.figure(figsize=(6,6))
+plt.plot(distance_axis,
+         density_axis)
+plt.xlabel("distance to stellar center (R$_{\star}$)",fontsize=14)
+plt.ylabel("iostropic wind density (1/cm$^3$)",fontsize=14)
+plt.yscale('log')
+
+plt.title("Isotropic wind density for %s" % (object_name))
+plt.show()
+# -
+
+# The temperature of the wind can be expressed as a power-law relation (Kochanek 1993; Bogomazov 2005), ignoring heating due to photoionisation.
+# $$
+#    T_{w}(r) = T_{\star} \left( \frac{r}{r_{\star}} \right)^{-\beta}
+# $$
+
+