diff --git a/.github/workflows/docker-image.yml b/.github/workflows/docker-image.yml old mode 100644 new mode 100755 diff --git a/.gitignore b/.gitignore new file mode 100755 index 0000000..7f3c645 --- /dev/null +++ b/.gitignore @@ -0,0 +1,33 @@ +# ---- Python bytecode / caches ---- +__pycache__/ +*.py[cod] +*$py.class + +# ---- Jupyter ---- +.ipynb_checkpoints/ +*.ipynb + +# ---- Build / packaging ---- +*.egg-info/ +dist/ +build/ +.eggs/ + +# ---- Virtual envs ---- +.venv/ +venv/ +env/ + +# ---- Temporary / local output / local test folders ---- +temp/ +tmp/ +*.log +/*.dat +/alpha_disc/tests_magn.py + +# ---- OS/editor junk ---- +.DS_Store +Thumbs.db +.vscode/ +.idea/ + diff --git a/Dockerfile b/Dockerfile old mode 100644 new mode 100755 diff --git a/LICENSE.md b/LICENSE.md old mode 100644 new mode 100755 diff --git a/README.md b/README.md old mode 100644 new mode 100755 index 6fd5142..b2e5749 --- a/README.md +++ b/README.md @@ -311,16 +311,64 @@ print(varkappa_C, rho_C, T_C, P_C) # Opacity, bulk density, temperature and gas print(Sigma0) # Surface density of disc print(vertstr.tau()) # optical thickness of disc print(vertstr.C_irr, vertstr.T_irr) # irradiation parameter and temperature + + +``` + +## Magnetic field treatment +To calculate the structure of an accretion disc around a magnetized object with its own magnetic +field, we added the argument `magn_args` containing all info related to the magnetic field. This is a dictionary +that can contain the following fields: +```python +magn_args = {'model': "KuzinLipunova", # model of the magnetic torque and truncation radius, + 'r_star': 1.2e6, # the star raduis; here - neutron star, + 'mu_magn': 1e28, # magnetic dipole moment of the star, + 'freq':200, # frequency of the star spin, + 'chi_deg':45, # the angle between the star rotation axis and the dipole axis, + 'r_out_magn': 'rlight', # for r > r_out_magn, the magnetic field is set to zero everywhere, + 'r_lower_limit': 1.3e6, 'r_upper_limit': 'rcor', # hard lower/upper limits for r_inner + ########### and some other model-specific keywords: ############### + 'eta': 1, 'kappa': 0.5, 'kappa_prime':0.5, 'delta0_in':0.05, + "ra_coef":0.5, 'kappa_alfven':0.1} +``` +The magnetic field support is added for all cases without an external irradiation. For every class of the vertical +structure with the name `SomeName`, there is a class `SomeNameMagnetic` supporting the argument `magn_args`. +An example of the vetical structure calculation with a magnetic field with tabular opacicities and EoS: +```python +from alpha_disc import mesa_vs + +M = 2e33 +alpha = 0.01 +r = 2e9 +F = 2e34 + +magn_args = {'model': 'KuzinLipunova', 'mu_magn': 1e27, 'eta': 0.8, 'r_star': 1.2e6, + 'freq':200, 'kappa': 0.02, 'chi_deg':30, 'kappa_prime':0.03, 'delta0_in':0.05, + 'r_out_magn': 'rlight', 'r_lower_limit': 1.3e6, 'r_upper_limit': 'rcor'} +vertstr = mesa_vs.MesaVerticalStructureRadConvMagnetic(Mx=M, alpha=alpha, r=r, F=F, + magn_args=magn_args) +z0r, result = vertstr.fit() +varkappa_C, rho_C, T_C, P_C, Sigma0 = vertstr.parameters_C() +print(z0r) +print(varkappa_C, rho_C, T_C, P_C) +print(Sigma0) +print(vertstr.tau()) + ``` + ## Structure choice Module `profiles`, containing `StructureChoice()` function, serves as interface for creating the right structure object in a simpler way. The input parameter of all structure classes is `F` - viscous torque. One can use other input parameters instead viscous torque `F` (such as effective temperature $T_{\rm eff}$ and accretion rate $\dot{M}$) using `input` parameter, and choose the structure type using `structure` parameter. The relation between $T_{\rm eff}, \dot{M}$ and $F$: ```math - \sigma_{\rm SB}T_{\rm eff}^4 = \frac{3}{8\pi} \frac{F\omega_{\rm K}}{r^2}, \quad F = \dot{M}h \left(1 - \sqrt{\frac{r_{\rm in}}{r}}\right) + F_{\rm in}, + \sigma_{\rm SB}T_{\rm eff}^4 = \frac{3}{8\pi} \frac{F\omega_{\rm K}}{r^2}, \quad F = \dot{M}h \left(1 - \sqrt{\frac{r_{\rm in}}{r}}\right) + F_{\rm in} + F_\mathrm{magn}, ``` -where $r_{\rm in} = 3r_{\rm g}=6GM/c^2$. The default value of viscous torque at the inner boundary of the disc $F_{\rm in}=0$ (it corresponds to Schwarzschild black hole as central source). If $F_{\rm in}\neq0$ you should set the non-zero value of $F_{\rm in}$ manually (`F_in` parameter) for correct calculation of the relation above. +where $r_{\rm in} = 3r_{\rm g}=6GM/c^2$ without the magnetic field or $r_\mathrm{in}=r_\mathrm{in}(\dot{M})$ with the magnetic field. The default value of viscous torque at the inner boundary of the disc +$F_{\rm in}=0$ (it corresponds to Schwarzschild black hole as central source). If $F_{\rm in}\neq0$ you should set the non-zero value of $F_{\rm in}$ manually (`F_in` parameter) for correct calculation of the relation above. + +**Note**: The examples below are given for the non-magnetic case, but all used methods: `StructureChoice`, `Vertical_Profile`, +`S_curve`, `Radial_Profile` support the magnetic argument `magn_args` described above. Usage: ``` python3 @@ -413,6 +461,20 @@ column density, temperature and energy flux in the disc and $\alpha$ is Shakura- parameter ([Shakura & Sunyaev 1973](https://ui.adsabs.harvard.edu/abs/1973A&A....24..337S)). By default, the inner viscous torque $F_{\rm in}=0$ (this case corresponds to a black hole as an accretor). +In the magnetic case, the first equation and the boundary conditions $P(z\_0)$ and $Q\_0$ are modified. +Namely, the equation of the hydromagnetostatics with the modified pressure boundary is +```math +\begin{split} +\frac{{\rm d}P}{{\rm d}z} &= -\rho\,\omega^2_{K} z - \frac{{\rm d}P_\mathrm{magn,~induced}}{{\rm d}z} \qquad\quad\,\, P_{\rm gas}(z_0) = P' + \Delta P_\mathrm{magn,~photosphere}; \\ +\end{split} +``` +where both magnetic modifications are usually known functions of $r$ and $z$, and the $Q\_0$ naturally changes: +```math +\begin{equation} + Q_0 = \frac{3}{8 \pi} \frac{F \omega_K}{r^2} = \frac{3}{8 \pi} \frac{\omega_K}{r^2} \left(F_\mathrm{no~field} + F_\mathrm{magn}\right). +\end{equation} +``` + The temperature gradient $\nabla\equiv\frac{{\rm d}\ln T}{{\rm d}\ln P}$ is defined according to the Schwarzschild criterion: ```math @@ -456,7 +518,7 @@ After the normalizing $P_{\rm gas}, Q, T, \Sigma$ on their characteristic values and replacing $z$ on $\hat{z} = 1 - z/z_0$ (in code it is the `t` variable), one has: ```math \begin{split} -\frac{{\rm d}\hat{P}}{{\rm d}\hat{z}} &= \frac{z_0^2}{P_0}\,\omega^2_{\rm K} \,\rho (1-\hat{z}) - \frac{4aT_0^3}{3 P_0} \frac{{\rm d}\hat{T}}{{\rm d}\hat{z}} \qquad\qquad \hat{P}(0) = P'/P_0; \\ +\frac{{\rm d}\hat{P}}{{\rm d}\hat{z}} &= \frac{z_0^2}{P_0}\,\omega^2_{\rm K} \,\rho (1-\hat{z}) - \frac{4aT_0^3}{3 P_0} \frac{{\rm d}\hat{T}}{{\rm d}\hat{z}} + \frac{z_0}{P_0} \frac{{\rm d}P_\mathrm{magn,~induced}}{{\rm d}z} \qquad\qquad \hat{P}(0) = (P' + \Delta P_\mathrm{magn})/P_0; \\ \frac{{\rm d}\hat{\Sigma}}{{\rm d}\hat{z}} &= 2\,\frac{z_0}{\Sigma_{00}}\,\rho \qquad\qquad\qquad\qquad\qquad\qquad\qquad \hat{\Sigma}(0) = 0; \\ \frac{{\rm d}\hat{T}}{{\rm d}\hat{z}} &= \nabla \frac{\hat{T}}{P_{\rm tot}} z_0^2\,\omega^2_{\rm K} \,\rho (1-\hat{z}) \quad\qquad\qquad\qquad \hat{T}(0) = T_{\rm eff}/T_0; \\ \frac{{\rm d}\hat{Q}}{{\rm d}\hat{z}} &= -\frac32\,\frac{z_0}{Q_0}\,\omega_{\rm K} \alpha P_{\rm tot} \qquad\qquad\qquad \hat{Q}(0) = 1, \quad \hat{Q}(1) = 0; \\ diff --git a/alpha_disc/__init__.py b/alpha_disc/__init__.py old mode 100644 new mode 100755 diff --git a/alpha_disc/magn.py b/alpha_disc/magn.py new file mode 100755 index 0000000..66bf84a --- /dev/null +++ b/alpha_disc/magn.py @@ -0,0 +1,632 @@ +#!/usr/bin/env python3 +import numpy as np +from astropy import constants as const +from scipy.optimize import brentq +from numpy import sin, cos, pi +import inspect + +G = const.G.cgs.value +R_GAS = const.R.cgs.value +K_BOLTZ = const.k_B.cgs.value +H_PLANCK = const.h.cgs.value +C_LIGHT = const.c.cgs.value +SIGMA_STEFBOLTZ = const.sigma_sb.cgs.value +R_SOLAR = const.R_sun.cgs.value +M_SOLAR = const.M_sun.cgs.value + +ALLOWED_R0_MODELS = ['alfven', 'kuzinlipunova', 'kluzniakrappaport07_i', + 'kluzniakrappaport07_ii', 'wang97', 'standard'] +ALLOWED_R_KEYWORDS = ['rlight', 'rcor'] + + +def _call_with_filtered_dict(func, data_dict, explicit_args=None): + sig = inspect.signature(func).parameters.keys() + if explicit_args is not None: + combined_args = data_dict | explicit_args + else: + combined_args = data_dict + final_kwargs = {k: v for k, v in combined_args.items() if k in sig} + return func(**final_kwargs) + + +def fn(dist0, dist, n = 1): + return 1 - (dist0 / dist)**(n/2) + +def r_cor(Mx, freq): + """ Corotation radius.""" + GM = G * Mx + return np.cbrt(GM / 4 / pi**2 / freq**2) + +def r_alfven(Mx, mu_magn, Mdot): + """Alfven radius (without any numerical coefs)""" + GM = G * Mx + return (mu_magn**4 / Mdot**2 / GM)**(1./7.) + +def r_isco(Mx): + """innermost stable circular orbit radius""" + return 6 * G * Mx / C_LIGHT**2 + +def r_light(freq): + """ Light cylinder radius""" + return C_LIGHT / 2 / pi / freq + +def f_p(eta, delta0_in, chi_deg): + """F_P function in Kuzin Lipunova model""" + chi = np.deg2rad(chi_deg) + aux_= (1 - eta)**2 + pi * np.sqrt(delta0_in*0.5) * eta * (1 - eta) + return cos(chi)**2 * aux_ + pi**2/16 * delta0_in * sin(chi)**2 + +def _require_not_none(**kwargs): + missing = [k for k, v in kwargs.items() if v is None] + if missing: + raise ValueError(f"Required params must not be None: {', '.join(missing)}") + + +def radius_inner(Mx, model, mu_magn=None, Mdot=None,r_star=None, chi_deg=None, freq=None, + r_lower_limit=None, r_upper_limit=None, + ra_coef=None, eta=None, kappa=None, delta0_in=None, + ): + """ + Compute the inner disk radius r0 for a magnetized central pbject under several + literature prescriptions, with optional limiting by lower/upper bounds. + + This routine evaluates one of the supported models for the truncation radius + (Alfvenic or torque-balance based), solves the corresponding algebraic equation + for the dimensionless fastness parameter w when needed, and returns + r0 = r_c * w**(2/3) or r0 = r_Alfven as appropriate, then enforces physically + motivated limits (ISCO and stellar surface) and any user-provided bounds. + + Parameters + ---------- + model : {'Alfven', 'KuzinLipunova', 'KluzniakRappaport07_I', + 'KluzniakRappaport07_II', 'Wang97', 'Standard'} + Choice of inner-radius prescription (case-insensitive): + - 'Standard': maximum of r_star and R_ISCO. + - 'Alfven' : r0 = ra_coef * r_A (Alfven radius defined without any numerical + coefficients, multiplied by ra_coef). + - 'kuzinlipunova' : Kuzin & Lipunova torque-balance form. + - 'KluzniakRappaport07_I' : Kluźniak & Rappaport (2007) case I + (ArXiv 0709.2361, eq. 17). + - 'KluzniakRappaport07_II' : Kluźniak & Rappaport (2007) case II + (ArXiv 0709.2361, eq. B7). + - 'Wang97' : Wang (1997) prescription + (https://ui.adsabs.harvard.edu/abs/1997ApJ...475L.135W/abstract, + eq. (7) with \Gamma === 1) + Mx : float + Stellar Mx [g]. + mu_magn : float, optional + Magnetic dipole moment of the star [G cm^3]. Default is None. + Mdot : float, optional + Mx accretion rate [g s^-1]. Default is None. + r_star : float or None, optional + Stellar radius (cm). Used as a hard lower limit together with ISCO. + If None, set to R_ISCO. Default is None. + chi_deg : float, optional + Angle chi between the magnetic dipole axis and the disc symmetry axis, + in degrees. Default is None. + freq : float, optional + Stellar spin frequency (Hz). Used for the corotation radius and light cylinder. + Default is None. + r_lower_limit : {float, 'rcor', 'rlight', None}, optional + Additional lower bound to apply to r0. If a string: + - 'rcor', use r_c (corotation radius). + - 'rlight', use r_L (light-cylinder radius). + If None, no extra lower bound is applied. Default is None. + r_upper_limit : {float, 'rcor', 'rlight', None}, optional + Additional upper bound to apply to r0. Same string semantics as + r_lower_limit. If None, no extra upper bound is applied. Default is None. + ra_coef : float, optional + Coefficient for the Alfven radius. Default is None. + eta : float, optional + Coupling/efficiency/penetration factor appearing in Kuzin-Lipunova and + Wang models. Default is None. + kappa : float, optional + Parameter entering the Kuzin–Lipunova forms. Default is None. + delta0_in : float, optional + Boundary layer thickness parameter (dimensionless) used in the models + of Kuzin-Lipunova and Wang97. + Default is None. + + Returns + ------- + float + The inner radius r0 in cm, after applying: + (1) the model formula, + (2) the hard lower limit max(r_star, r_ISCO), + (3) the optional r_lower_limit and r_upper_limit if provided. + + + """ + + if model.lower() not in ALLOWED_R0_MODELS: + raise ValueError(f"model of the inner radius should be one of: {ALLOWED_R0_MODELS}.") + + risco = r_isco(Mx=Mx) + if r_star is None: + r_star = risco + rmax = max(r_star, risco) + if model.lower() == 'standard': + return rmax + _require_not_none( mu_magn=mu_magn, Mdot=Mdot, r_star=r_star) + _ra = r_alfven(Mdot=Mdot, Mx=Mx, mu_magn=mu_magn) + if model.lower() != 'alfven': + _require_not_none(freq=freq) + _rc = r_cor(Mx=Mx, freq=freq) + _ra_rc = _ra / _rc # the same as \xi in Kluzniak, Rappaport + + if model.lower() == 'alfven': + _require_not_none(ra_coef=ra_coef) + _r0 = _ra * ra_coef + + if model.lower() == 'kuzinlipunova': + _require_not_none(eta=eta, chi_deg=chi_deg, delta0_in=delta0_in, + kappa=kappa) + _chi = np.deg2rad(chi_deg) + _s2 = sin(_chi)**2 + _c2 = 1 - _s2 + _f_p = f_p(eta=eta, delta0_in=delta0_in, chi_deg=chi_deg) + + if eta < 1e-3: ### special case eta=0 + w = _ra_rc**1.5 * (2 * kappa * _f_p)**(3/7) + else: + one = lambda _w: (1-_w) * eta**2 + two = lambda _w: (1-_w) * eta**2 * (1 - 3 * delta0_in) * _s2 + three = lambda _w: kappa * _f_p + to_solve = lambda _w: 0.5 - _ra_rc**3.5 * _w**(-7/3) * (one(_w) + two(_w) + three(_w)) + try: + w = brentq(to_solve, 1e-10, 10000) + except: + w = 0 + _r0 = _rc * w**(2/3) + + if model.lower() == 'kluzniakrappaport07_i': + to_solve = lambda _w: 0.5 - _ra_rc**3.5 * _w**(-10/3) * (1 - _w) + try: + w = brentq(to_solve, 1e-10, 10000) + except: + w = 0 + _r0 = _rc * w**(2/3) + + if model.lower() == 'kluzniakrappaport07_ii': + to_solve = lambda _w: 0.5 - _ra_rc**3.5 * _w**(-7/3) * (1 - _w) + try: + w = brentq(to_solve, 1e-10, 10000) + except: + w = 0 + _r0 = _rc * w**(2/3) + + if model.lower() == 'wang97': + _require_not_none(eta=eta, chi_deg=chi_deg, + delta0_in=delta0_in) + _chi = np.deg2rad(chi_deg) + _s2 = sin(_chi)**2 + _c2 = 1 - _s2 + one = lambda _w: (1 - _w) * _c2 + two = lambda _w: (8 - 5 * _w) * delta0_in * _s2 + + to_solve = lambda _w: 0.5 - eta**2 * _ra_rc**3.5 * _w**(-7/3) * (one(_w) + two(_w)) + try: + w = brentq(to_solve, 1e-10, 10000) + except: + w = 0 + _r0 = _rc * w**(2/3) + + _r0 = np.maximum(rmax, _r0) + + if r_lower_limit is not None: + if isinstance(r_lower_limit, str) and r_lower_limit not in ALLOWED_R_KEYWORDS: + raise ValueError(f"If r_lower_limit is a string, it should be one of {ALLOWED_R_KEYWORDS}") + if isinstance(r_lower_limit, str) and r_lower_limit == 'rcor': + r_lower_limit = _rc + if isinstance(r_lower_limit, str) and r_lower_limit == 'rlight': + r_lower_limit = r_light(freq=freq) + _r0 = np.maximum(r_lower_limit, _r0) + + if r_upper_limit is not None: + if isinstance(r_upper_limit, str) and r_upper_limit not in ALLOWED_R_KEYWORDS: + raise ValueError(f"If r_upper_limit is a string, it should be one of {ALLOWED_R_KEYWORDS}") + if isinstance(r_upper_limit, str) and r_upper_limit == 'rcor': + r_upper_limit = _rc + if isinstance(r_upper_limit, str) and r_upper_limit == 'rlight': + r_upper_limit = r_light(freq=freq) + _r0 = np.minimum(r_upper_limit, _r0) + + return _r0 + +def magnetic_torque(r, mu_magn=None, Mx=None, r0=None, model=None, freq=None, + r_out_magn=None, eta=None, kappa_prime=None, + delta0_in=None, chi_deg=None, kappa_alfven=None): + """ + Compute the magnetic torque (or magnetic stress term) as a function of radius + for several commonly used magnetosphere–disk coupling prescriptions. + + The function returns a radial profile :math:`F_\\mathrm{magn}(r)` evaluated for + the chosen `model`. For the `'standard'` model the magnetic torque is set to + zero. For all other models the input radius (or radii) must satisfy + r >= r_0. + + Parameters + ---------- + r : float or array_like + Radius (or radii) at which to evaluate the magnetic torque. Must be + greater than or equal to `r0`. Scalars return a scalar, arrays return an + array of matching shape. + mu_magn : float + Magnetic dipole moment of the central object. + Mx : float + Mass of the central object (passed to the corotation-radius helper + `r_cor` when needed). + r0 : float + The inner disk radius. + model : {'Standard', 'Alfven', 'KuzinLipunova', 'KluzniakRappaport07_I', \ + 'KluzniakRappaport07_II', 'Wang97'} + Choice of magnetic-torque prescription. Matching is case-insensitive. + + - 'standard': returns identically zero torque. + - 'alfven': constant torque level set by `kappa_alfven`. + - 'kuzinlipunova': Kuzin & Lipunova-type expression using + `eta`, `kappa_prime`, `delta0_in`, and `chi_deg`. + - 'kluzniakrappaport07_i': Kluzniak & Rappaport (2007), case I. + (ArXiv 0709.2361, eq. 24). + - 'kluzniakrappaport07_ii': Kluzniak & Rappaport (2007), case II + (piecewise below/above corotation). (ArXiv 0709.2361, eq. B5, B6). + - 'wang97': Wang (1997)-type expression using `eta`, `delta0_in`, + and `chi_deg`. + (https://ui.adsabs.harvard.edu/abs/1997ApJ...475L.135W/abstract, + eqs. 8, 9.) + freq : float, optional + Spin frequency of the central object. Required for all models except + 'standard' and 'alfven'. Used to compute the corotation radius + and the fastness parameter. + r_out_magn : None or float or {'rcor', 'rlight'}, optional + Optional outer radius for the magnetic interaction region. + If provided, the torque is evaluated with an effective distance + dist_ = min(r, r_out_magn). If a string keyword is given: + + - 'rcor': uses the corotation radius computed from `freq`. + - 'rlight': uses the light-cylinder radius computed from `freq`. + + If not provided, ``dist_ = r``. + eta : float, optional + Dimensionless parameter used in some prescriptions ( + ``'kuzinlipunova'``, ``'wang97'``). Required for those models. + kappa_prime : float, optional + Dimensionless coupling parameter used in the 'kuzinlipunova' model. + delta0_in : float, optional + The effective z0/r at the inner edge f=of the disc. + Required for 'kuzinlipunova' and 'wang97'. + chi_deg : float, optional + Magnetic obliquity angle in degrees. Required for models that include an + inclination dependence ('kuzinlipunova' and 'wang97'). + kappa_alfven : float, optional + Dimensionless normalization for the 'alfven' model. Required if + model='alfven'. + + Returns + ------- + F_magn : float or numpy.ndarray + Magnetic torque (or magnetic stress term) evaluated at `r`. If `r` is a + scalar, a Python float is returned. If `r` is array-like, a NumPy array + with the same shape as `r` is returned. + + """ + if model.lower() not in ALLOWED_R0_MODELS: + raise ValueError(f"model of the magnetic torque should be one of: {ALLOWED_R0_MODELS}.") + if model.lower() == 'standard': + return 0 * r + _require_not_none(Mx=Mx, mu_magn=mu_magn, r0=r0, model=model) + r = np.asarray(r) + if np.ndim(r) == 0: + if r < r0: + raise ValueError(f"r = {r} in F_magn should be > r0 = {r0}.") + else: + if any([r_i < r0 for r_i in r]): + raise ValueError("all r in F_magn should be > r0.") + + if model.lower() != 'alfven' or isinstance(r_out_magn, str): + _require_not_none(freq=freq) + _rc = r_cor(Mx=Mx, freq=freq) + w_fastness = (r0 / _rc)**1.5 + + dist_ = r + + if r_out_magn is not None: + if isinstance(r_out_magn, str) and r_out_magn not in ALLOWED_R_KEYWORDS: + raise ValueError(f"If r_out_magn is a string, it should be one of: {ALLOWED_R_KEYWORDS}") + if isinstance(r_out_magn, str) and r_out_magn == 'rcor': + r_out_magn = _rc + if isinstance(r_out_magn, str) and r_out_magn == 'rlight': + r_out_magn = r_light(freq=freq) + # print('r out magn = ', r_out_magn) + dist_ = np.minimum(r, r_out_magn) + + + if model.lower() == 'alfven': + # print('alf') + _require_not_none(kappa_alfven = kappa_alfven) + _F_magn = kappa_alfven * mu_magn**2 / r0**3 * np.ones(r.size) + + if model.lower() == 'kuzinlipunova': + # print('kl') + _require_not_none(eta=eta, kappa_prime=kappa_prime, chi=chi_deg, + delta0_in=delta0_in) + chi = np.deg2rad(chi_deg) + f3_ = fn(r0, dist_, 3) + f6_ = fn(r0, dist_, 6) + _g_func = (eta**2 / 3 * (2 * w_fastness * f3_ - f6_) * cos(chi)**2 + + 2 * kappa_prime * f_p(eta=eta, delta0_in=delta0_in, chi_deg=chi_deg)) + + _F_magn = mu_magn**2 / r0**3 * _g_func + + if model.lower() == 'kluzniakrappaport07_i': + _F_magn = mu_magn**2/9 * (3 / r0**3 + - 3 / dist_**3 + - 2 * _rc**1.5 / r0**4.5 + + 2 * _rc**1.5 / dist_**4.5) + + if model.lower() == 'kluzniakrappaport07_ii': + def _F_magn_under_cor(_r): + return mu_magn**2 / 3 / _r**3 * (2 * (_r / _rc)**1.5 * ( + (_r / r0)**1.5 - 1) + + 1 + - (_r / r0)**3) + def _F_magn_over_cor(_r): + return mu_magn**2 / 9 / _r**3 * ( - 2 * (_r / _rc)**3 + + 2 * (_rc / _r)**1.5 + + 6 * (_r**2 / r0 / _rc)**1.5 + - 3 * (_r / r0)**3 + - 3 + ) + _F_magn = np.where(dist_ < _rc, _F_magn_under_cor(dist_), _F_magn_over_cor(dist_)) + + if model.lower() == 'wang97': + _require_not_none(chi_deg=chi_deg, eta=eta, delta0_in=delta0_in) + _s2 = sin(np.deg2rad(chi_deg))**2 + _F_magn = mu_magn**2 / r0**3 * ( (1 - _s2) * 1/3 * ( + 2 * w_fastness * fn(r0, dist_, 3) - fn(r0, dist_, 6)) + + 2 * _s2 * delta0_in * (w_fastness - 1) + ) + + return _F_magn.item() if np.ndim(r)==0 else _F_magn + + +class MagneticField: + """ + A mixin for the magnetic field treatment in the vertical structure + calculation. + + Takes the dictionary `magn_args` as an argument. If it's None or if + it contains an entry ['model'] == 'standard', the behavior falls to a + fully non-magnetic case. Essentially, the behavior is non-magnetic + in all cases except model == 'KuzinLipunova'. The dictionary, if provided, + should always contain 'model' and 'r_star'. + + It overwrites the two BaseVerticalStructure functions: + - additional_photospheric_pressure() --- for the aditional + pressure the discontinuous magnetic field can exert on a + photosphere. Set == 0 in all cases except model == 'KuzinLipunova'. + + - dP_magn_induced_dz --- for the additional source in the equation + of vertical hydtrostatics. It should be a d/dz(b_phi^2/8pi), + where b_phi is a phi-component of the induced magnetic field inside + the disc. In case model=='KuzinLipunova', set using the approximation + b_phi(r, z) = b_phi(r, z=z0) * z/z0. + + Inside, it calls the ViscousTorque class to calculate the inner radius + of the disk given F. I guess this should be fixed in future when we have the + AccretionDisc class containing all info... + + """ + def __init__(self, *args, **kwargs): + magn_args = kwargs.pop("magn_args", None) + super().__init__(*args, **kwargs) + + if magn_args is None: + self.model = 'standard' + magn_args = {'model': 'standard'} + else: + if magn_args['model'].lower() not in ALLOWED_R0_MODELS: + raise ValueError(f"""The model in magn_args in the class MagneticField + \nshould be one of:{ALLOWED_R0_MODELS}.""") + self.model = magn_args['model'] + self.magn_args = magn_args + + if self.model != 'standard': + self._set_magn_pars(magn_args=magn_args) + + def _set_magn_pars(self, magn_args): + self.mu_magn = magn_args['mu_magn'] + self.model = magn_args['model'] + self.r_star = magn_args['r_star'] + + if self.model.lower() == 'alfven': + self.ra_coef = magn_args['ra_coef'] + if self.model.lower() == 'kuzinlipunova': + self.freq = magn_args['freq'] + self.chi_deg = magn_args['chi_deg'] + self.chi = np.deg2rad(self.chi_deg) + self.eta = magn_args['eta'] + self.kappa = magn_args['kappa'] + self.kappa_prime = magn_args['kappa_prime'] + self.delta0_in = magn_args['delta0_in'] + + + + def _f_p(self): + return f_p(eta=self.eta, delta0_in=self.delta0_in, chi_deg=self.chi_deg) + + def _r_out_magn(self, Mx): + """An analogue of r_out_magn, but since it depends on Mx which is + not known at the init time of this mixin, I made it a funciton of Mx.""" + _r_out = self.magn_args['r_out_magn'] + if _r_out is not None: + if isinstance(_r_out, str) and _r_out not in ALLOWED_R_KEYWORDS: + raise ValueError(f"If _r_out is a string, it should be one of {ALLOWED_R_KEYWORDS}") + if isinstance(_r_out, str) and _r_out == 'rcor': + _r_out = r_cor(Mx=Mx, freq=self.freq) + if isinstance(_r_out, str) and _r_out == 'rlight': + _r_out = r_light(freq=self.freq) + return _r_out + + @property + def _r0(self): + _Torques = ViscousTorque(Mx=self.Mx, magn_args=self.magn_args) + _Mdot = _Torques.Mdot_from_F(F=self.F, r=self.r) + return _Torques.r0(_Mdot) + + def additional_photospheric_pressure(self): + if self.model.lower() == 'kuzinlipunova': + _s2 = sin(self.chi)**2 + _c2 = 1 - _s2 # cos^2(chi) + D = max(( (2 * self.delta0_in)**0.5, + (self.r**2 / self._r0**2 - 1)**0.5 )) # Lai1999 + P_aly = self.mu_magn ** 2 / self.r ** 6 / 8 / pi * ( + 2.5 * _s2 + (4 / pi / D *(1 - self.eta)) ** 2 * _c2 + ) + return P_aly + else: + return 0 + + def _b_induced_surface(self, _r): + _r = np.asarray(_r) + if self.model.lower() == 'kuzinlipunova': + _rc = r_cor(Mx=self.Mx, freq=self.freq) + R = np.where( + (_r < self._r_out_magn(Mx=self.Mx)), + (self.mu_magn / _r**3) * (1.0 - (_r / _rc)**1.5), + 0.0 + ) * self.eta * cos(self.chi) + return R.item() if np.isscalar(_r) else R + else: + return _r * 0 + + def dP_magn_induced_dz(self, z): + if self.model.lower() == 'kuzinlipunova': + ### uses approximation: b(r, z) = b_surface(r) * z/z0, + ### so d(b^2)/dz = b_surface^2 * 2z/z0^2 + db2_dz = self._b_induced_surface(_r=self.r)**2 * 2 * z / self.z0**2 + return db2_dz / 8 / pi + else: + return z*0 + + +class ViscousTorque: + """ + A class sumarizing the information about the visocus and magnetic torques + distribution over radius, all that can be found prior to soling for the + vertical structure. All quantities are in CGS. + + All methods are vectorized over r and not vectorized over Mdot. + + Attributes + ---------- + Mx : double + Mass of the central object. + F_in : double + Inner torque. Default 0. Used only if it was impossible to solve for + Mdot given F or T_eff, so the simplified standard disc expressions are + used. + GM : double + Yeah. + rg : double + Grav radius 2GM/c^2. + magn_args : dict or None + A dictionary of arguments for finding the inenr radius and the magnetic + distributed torque. + + Methods + ------- + r0(Mdot): + Finds the inner radius of the disc by calling the function radius_inner. + See the documentation for this function. + F_standard(Mdot, r): + The Shalura-Sunyaev viscous torque: Mdot sqrt(GMr) (1 - sqrt(r0/r)). + F_magn(Mdot, r): + The magnetic torque. The function magnetic_torque if called, see its + documentation. + F_vis(Mdot, r): + The total torque, the sum of F_standard, F_magn, and F_in. + Q_vis(Mdot, r): + Viscous heat flux from the unit of one surface. + Q_vis = 3/4 W_rphi w_kepl. + T_eff(Mdot, r): + The effective temperature of a disc. + sigma_sb T_eff^4 = Q_vis. + Mdot_from_F(F, r): + Finds such an Mdot as to give the F at the radius r. Tries to solve the + equation for Mdot, and it fails, returns the analytical standard expression. + Mdot_from_Teff(Teff, r): + Finds such an Mdot as to give the Teff at the radius r. Tries to solve the + equation for Mdot, and it fails, returns the analytical standard expression. + + """ + def __init__(self, Mx, F_in=0, magn_args=None): + self.Mx = Mx + self.GM = G * Mx + self.rg = 2 * self.GM / C_LIGHT**2 + self.F_in = F_in + if magn_args is None: + self.model = 'standard' + magn_args = {'model': 'standard'} + else: + self.model = magn_args['model'] + if self.model.lower() not in ALLOWED_R0_MODELS: + raise ValueError(f"""The model in magn_args in the class ViscousTorque + \nshould be one of: {ALLOWED_R0_MODELS}.""") + self.magn_args = magn_args + + def r0(self, Mdot): + return _call_with_filtered_dict(func=radius_inner, + explicit_args = {'Mx': self.Mx, 'Mdot':Mdot}, + data_dict=self.magn_args) + + def _fn(self, Mdot, r, n=1): + return fn(dist0 = self.r0(Mdot), dist=r, n=n) + + def _w_kepl(self, r): + return np.sqrt(self.GM / r**3) + + def _h(self, r): + return np.sqrt(self.GM * r) + + def F_standard(self, Mdot, r): + return self._h(r) * self._fn(r=r, Mdot=Mdot, n=1) * Mdot + + def F_magn(self, Mdot, r): + return _call_with_filtered_dict(func=magnetic_torque, + explicit_args = {'r':r,'Mx':self.Mx, 'r0':self.r0(Mdot)}, + data_dict=self.magn_args + ) + + def F_vis(self, Mdot, r): + return self.F_standard(Mdot, r) + self.F_magn(Mdot, r) + self.F_in + + def Q_vis(self, Mdot, r): + ### 3/4 W_rphi w_kepl = 3/4 (F / 2 pi r^2) w_kepl + return 0.75 * (self.F_vis(Mdot, r) / 2 / pi / r**2) * self._w_kepl(r) + + def T_eff(self, Mdot, r): + return (self.Q_vis(Mdot, r) / SIGMA_STEFBOLTZ)**(1 / 4) + + def Mdot_from_F(self, F, r): + try: + return brentq(f = lambda _mdot: F - self.F_standard(_mdot, r), + a=1.0, b=1e22) + except ValueError: + print("Could not solve for Mdot given F; returning Mdot from simple standard expression.") + return (F - self.F_in) / self.F_standard(Mdot=1.0, r=r) + + def Mdot_from_Teff(self, Teff, r): + try: + return brentq(f = lambda _mdot: Teff - self.T_eff(_mdot, r), + a=1.0, b=1e22) + except ValueError: + print("Could not solve for Mdot given Teff; returning Mdot from simple standard expression.") + _F = 8 * pi / 3 * self._h(r) ** 7 / self.GM ** 4 * SIGMA_STEFBOLTZ * Teff ** 4 + return (_F - self.F_in) / self.F_standard(Mdot=1.0, r=r) + + + + + + \ No newline at end of file diff --git a/alpha_disc/mesa_vs.py b/alpha_disc/mesa_vs.py old mode 100644 new mode 100755 index bbf274d..dc1265b --- a/alpha_disc/mesa_vs.py +++ b/alpha_disc/mesa_vs.py @@ -36,7 +36,7 @@ from scipy.optimize import root, least_squares, brentq from alpha_disc.vs import BaseVerticalStructure, Vars, IdealGasMixin, RadiativeTempGradient, PgasPradNotConvergeError - +from alpha_disc.magn import MagneticField sigmaSB = const.sigma_sb.cgs.value sigmaT = const.sigma_T.cgs.value atomic_mass = const.u.cgs.value @@ -85,8 +85,8 @@ def __str__(self): class BaseMesaVerticalStructure(BaseVerticalStructure): - def __init__(self, Mx, alpha, r, F, eps=1e-5, mu=0.6, abundance='solar'): - super().__init__(Mx, alpha, r, F, eps=eps, mu=mu) + def __init__(self, Mx, alpha, r, F, eps=1e-5, mu=0.6, abundance='solar', **kwargs): + super().__init__(Mx, alpha, r, F, eps=eps, mu=mu, **kwargs) self.mesaop = Opac(abundance) @@ -220,7 +220,8 @@ class RadiativeAdiabaticGradient: def dlnTdlnP(self, y, t): rho, eos = self.rho(y, True) varkappa = self.opacity(y, lnfree_e=eos.lnfree_e) - P_full = y[Vars.P] * self.P_norm + 4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4 + # P_full = y[Vars.P] * self.P_norm + 4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4 + P_full = y[Vars.P] * self.P_norm + self.rad_pressure(y[Vars.T] * self.T_norm) if t == 1: dTdz_der = (self.dQdz(y, t) / y[Vars.T] ** 3) * 3 * varkappa * rho * self.z0 * self.Q_norm / ( @@ -250,7 +251,7 @@ class RadConvTempGradient: def dlnTdlnP(self, y, t): rho, eos = self.rho(y, True) varkappa = self.opacity(y, lnfree_e=eos.lnfree_e) - P_full = y[Vars.P] * self.P_norm + 4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4 + P_full = y[Vars.P] * self.P_norm + self.rad_pressure(y[Vars.T] * self.T_norm) if t == 1: dTdz_der = (self.dQdz(y, t) / y[Vars.T] ** 3) * 3 * varkappa * rho * self.z0 * self.Q_norm / ( @@ -709,6 +710,16 @@ class MesaVerticalStructure(MesaGasMixin, MesaOpacityMixin, RadiativeTempGradien """ pass +# Classes with MESA opacity without Irradiation +class MesaVerticalStructureMagnetic(MagneticField, MesaGasMixin, MesaOpacityMixin, + RadiativeTempGradient, BaseMesaVerticalStructure): + """ + Vertical structure class for MESA opacities and EoS with radiative energy transport. + Includes the magnetic field. + + """ + pass + class MesaIdealGasVerticalStructure(IdealGasMixin, MesaOpacityMixin, RadiativeTempGradient, BaseMesaVerticalStructure): """ @@ -717,6 +728,14 @@ class MesaIdealGasVerticalStructure(IdealGasMixin, MesaOpacityMixin, RadiativeTe """ pass +class MesaIdealGasVerticalStructureMagnetic(MagneticField, IdealGasMixin, + MesaOpacityMixin, RadiativeTempGradient, BaseMesaVerticalStructure): + """ + Vertical structure class for MESA opacities and ideal gas EoS with radiative energy transport. + Includes the magnetic field. + + """ + pass class MesaVerticalStructureAd(MesaGasMixin, MesaOpacityMixin, AdiabaticTempGradient, BaseMesaVerticalStructure): """ @@ -725,6 +744,16 @@ class MesaVerticalStructureAd(MesaGasMixin, MesaOpacityMixin, AdiabaticTempGradi """ pass +class MesaVerticalStructureAdMagnetic(MagneticField, MesaGasMixin, MesaOpacityMixin, + AdiabaticTempGradient, BaseMesaVerticalStructure): + """ + Vertical structure class for MESA opacities and EoS with adiabatic energy transport. + Includes the magnetic field. + + """ + pass + + class MesaVerticalStructureRadAd(MesaGasMixin, MesaOpacityMixin, RadiativeAdiabaticGradient, BaseMesaVerticalStructure): """ @@ -733,6 +762,16 @@ class MesaVerticalStructureRadAd(MesaGasMixin, MesaOpacityMixin, RadiativeAdiaba """ pass +class MesaVerticalStructureRadAdMagnetic(MagneticField, MesaGasMixin, MesaOpacityMixin, + RadiativeAdiabaticGradient, BaseMesaVerticalStructure): + """ + Vertical structure class for MESA opacities and EoS with radiative+adiabatic energy transport. + Includes the magnetic field. + + + """ + pass + class MesaVerticalStructureRadConv(MesaGasMixin, MesaOpacityMixin, RadConvTempGradient, BaseMesaVerticalStructure): """ @@ -741,6 +780,15 @@ class MesaVerticalStructureRadConv(MesaGasMixin, MesaOpacityMixin, RadConvTempGr """ pass +class MesaVerticalStructureRadConvMagnetic(MagneticField, MesaGasMixin, + MesaOpacityMixin, RadConvTempGradient, BaseMesaVerticalStructure): + """ + Vertical structure class for MESA opacities and EoS with radiative+convective energy transport. + Includes the magnetic field. + + """ + pass + # Classes with MESA opacity and with advanced Irradiation (Mescheryakov et al. 2011) class MesaVerticalStructureExternalIrradiation(MesaGasMixin, MesaOpacityMixin, RadiativeTempGradient, diff --git a/alpha_disc/profiles.py b/alpha_disc/profiles.py old mode 100644 new mode 100755 index f49a3b5..82cb022 --- a/alpha_disc/profiles.py +++ b/alpha_disc/profiles.py @@ -22,7 +22,8 @@ from astropy import constants as const from scipy.integrate import simpson from scipy.interpolate import InterpolatedUnivariateSpline - +from collections import namedtuple +from alpha_disc.magn import ViscousTorque import alpha_disc.vs as vert try: @@ -38,7 +39,8 @@ def StructureChoice(M, alpha, r, Par, input, structure, mu=0.6, abundance='solar', F_in=0, nu_irr=None, L_X_irr=None, spectrum_irr=None, spectrum_irr_par=None, args_spectrum_irr=(), kwargs_spectrum_irr={}, - C_irr=None, T_irr=None, cos_theta_irr=None, cos_theta_irr_exp=1 / 12): + C_irr=None, T_irr=None, cos_theta_irr=None, cos_theta_irr_exp=1 / 12, + magn_args=None, rad_pressure_include=True): """ Initialise the chosen vertical structure object. @@ -114,6 +116,11 @@ def StructureChoice(M, alpha, r, Par, input, structure, mu=0.6, abundance='solar T_irr : double If structure in ['MesaIrrZero', 'MesaRadAdIrrZero', 'MesaRadConvIrrZero'], T_irr is the irradiation temperature. + magn_args : dict + If the magnetic field treatment is to be included, this dicionary + should comtain the parameters necessary for the calculation of the + magneic torque, inner radius, etc. Only possible for models without + irratiation. Default None. Returns ------- @@ -127,74 +134,119 @@ def StructureChoice(M, alpha, r, Par, input, structure, mu=0.6, abundance='solar Accretion rate in g/s. """ - h = np.sqrt(G * M * r) - rg = 2 * G * M / c ** 2 - r_in = 3 * rg - if r <= r_in: - raise Exception(f'Radius r should be greater than inner radius r_in = 3*rg. ' + Torques = ViscousTorque(Mx=M, F_in=F_in, magn_args=magn_args) + rg = Torques.rg + # r_in = Torques.r0 + if r <= 3*rg: + raise Exception(f'Radius r should be greater than R_ISCO = 3*rg. ' f'Actual radius r = {r / rg:g} rg.') - func = 1 - np.sqrt(r_in / r) if input == 'Teff': Teff = Par - F = 8 * np.pi / 3 * h ** 7 / (G * M) ** 4 * sigmaSB * Par ** 4 - Mdot = (F - F_in) / (h * func) + Mdot = Torques.Mdot_from_Teff(Teff=Teff, r=r) if Mdot < 0: - raise Exception(f'Mdot = {Mdot:g} g/s < 0, incorrect F_in = {F_in:g} g*cm^2/s^2.') - elif input == 'Mdot': - Mdot = Par - F = Par * h * func + F_in - Teff = (3 / (8 * np.pi) * (G * M) ** 4 * F / (sigmaSB * h ** 7)) ** (1 / 4) - elif input == 'Mdot_Mdot_edd': - Mdot = Par * 1.39e18 * M / M_sun - F = Mdot * h * func + F_in - Teff = (3 / (8 * np.pi) * (G * M) ** 4 * F / (sigmaSB * h ** 7)) ** (1 / 4) - elif input == 'Mdot_Msun_yr': - Mdot = Par * M_sun / 31557600.0 - F = Mdot * h * func + F_in - Teff = (3 / (8 * np.pi) * (G * M) ** 4 * F / (sigmaSB * h ** 7)) ** (1 / 4) + raise Exception(f'Mdot = {Mdot:g} g/s < 0, incorrect F_in = {F_in:g} g*cm^2/s^2.') + F = Torques.F_vis(Mdot, r) + elif input in ('Mdot', 'Mdot_Mdot_edd', 'Mdot_Msun_yr'): + if input == 'Mdot': Mdot = Par + if input == 'Mdot_Mdot_edd': Mdot = Par * 1.39e18 * M / M_sun + if input == 'Mdot_Msun_yr': Mdot = Par * M_sun / 31557600.0 + F = Torques.F_vis(Mdot, r) + Teff = Torques.T_eff(Mdot, r) elif input == 'F': F = Par - Mdot = (Par - F_in) / (h * func) + Mdot = Torques.Mdot_from_F(F=F, r=r) if Mdot < 0: raise Exception(f'Mdot = {Mdot:g} g/s < 0, incorrect F_in = {F_in:g} g*cm^2/s^2.') - Teff = (3 / (8 * np.pi) * (G * M) ** 4 * Par / (sigmaSB * h ** 7)) ** (1 / 4) + Teff = Torques.T_eff(Mdot, r) else: raise Exception("Incorrect input, try 'Teff', 'F', 'Mdot', 'Mdot_Mdot_edd' or 'Mdot_Msun_yr'.") + + if magn_args is not None and structure in ('MesaIrrZero', 'MesaRadAdIrrZero', + 'MesaRadConvIrrZero', 'MesaIrr', + 'MesaRadAdIrr', 'MesaRadConvIrr'): + print("""A structure with irradiation was chosen, but the magnetic field is not supported for this case. + \nMagnetic field is set to zero.""") + r_in = Torques.r0(Mdot) + if r <= r_in: + raise Exception(f'Radius r should be greater than the inner radius = {r_in:g} cm. ' + f'Actual radius r = {r / r_in:g} r_in.') if structure == 'Kramers': - vs = vert.IdealKramersVerticalStructure(M, alpha, r, F, mu=mu) + if magn_args is None: + vs = vert.IdealKramersVerticalStructure(M, alpha, r, F, mu=mu, + rad_pressure_include=rad_pressure_include) + else: + vs = vert.IdealKramersVerticalStructureMagnetic(M, alpha, r, F, mu=mu, + magn_args=magn_args, + rad_pressure_include=rad_pressure_include) elif structure == 'BellLin': - vs = vert.IdealBellLin1994VerticalStructure(M, alpha, r, F, mu=mu) + if magn_args is None: + vs = vert.IdealBellLin1994VerticalStructure(M, alpha, r, F, mu=mu, + rad_pressure_include=rad_pressure_include) + else: + vs= vert.IdealBellLin1994VerticalStructureMagnetic(M, alpha, r, F, mu=mu, + magn_args=magn_args, + rad_pressure_include=rad_pressure_include) elif structure == 'Mesa': try: if np.isnan(mesa_vs): raise ModuleNotFoundError('Mesa2py is not installed') except TypeError: - vs = mesa_vs.MesaVerticalStructure(M, alpha, r, F, abundance=abundance) + if magn_args is None: + vs = mesa_vs.MesaVerticalStructure(M, alpha, r, F, abundance=abundance, + rad_pressure_include=rad_pressure_include) + else: + vs = mesa_vs.MesaVerticalStructureMagnetic(M, alpha, r, F, abundance=abundance, + magn_args=magn_args, + rad_pressure_include=rad_pressure_include) elif structure == 'MesaIdealGas': try: if np.isnan(mesa_vs): raise ModuleNotFoundError('Mesa2py is not installed') except TypeError: - vs = mesa_vs.MesaIdealGasVerticalStructure(M, alpha, r, F, mu=mu, abundance=abundance) + if magn_args is None: + vs = mesa_vs.MesaIdealGasVerticalStructure(M, alpha, r, F, mu=mu, abundance=abundance, + rad_pressure_include=rad_pressure_include) + else: + vs = mesa_vs.MesaIdealGasVerticalStructureMagnetic(M, alpha, r, F, mu=mu, abundance=abundance, + magn_args = magn_args, + rad_pressure_include=rad_pressure_include) elif structure == 'MesaAd': try: if np.isnan(mesa_vs): raise ModuleNotFoundError('Mesa2py is not installed') except TypeError: - vs = mesa_vs.MesaVerticalStructureAd(M, alpha, r, F, abundance=abundance) + if magn_args is None: + vs = mesa_vs.MesaVerticalStructureAd(M, alpha, r, F, abundance=abundance, + rad_pressure_include=rad_pressure_include) + else: + vs = mesa_vs.MesaVerticalStructureAdMagnetic(M, alpha, r, F, abundance=abundance, + magn_args=magn_args, + rad_pressure_include=rad_pressure_include) elif structure == 'MesaRadAd': try: if np.isnan(mesa_vs): raise ModuleNotFoundError('Mesa2py is not installed') except TypeError: - vs = mesa_vs.MesaVerticalStructureRadAd(M, alpha, r, F, abundance=abundance) + if magn_args is None: + vs = mesa_vs.MesaVerticalStructureRadAd(M, alpha, r, F, abundance=abundance, + rad_pressure_include=rad_pressure_include) + else: + vs = mesa_vs.MesaVerticalStructureRadAdMagnetic(M, alpha, r, F, abundance=abundance, + magn_args=magn_args, + rad_pressure_include=rad_pressure_include) elif structure == 'MesaRadConv': try: if np.isnan(mesa_vs): raise ModuleNotFoundError('Mesa2py is not installed') except TypeError: - vs = mesa_vs.MesaVerticalStructureRadConv(M, alpha, r, F, abundance=abundance) + if magn_args is None: + vs = mesa_vs.MesaVerticalStructureRadConv(M, alpha, r, F, abundance=abundance, + rad_pressure_include=rad_pressure_include) + else: + vs = mesa_vs.MesaVerticalStructureRadConvMagnetic(M, alpha, r, F, abundance=abundance, + magn_args=magn_args, + rad_pressure_include=rad_pressure_include) elif structure == 'MesaIrrZero': try: if np.isnan(mesa_vs): @@ -316,7 +368,9 @@ def Vertical_Profile(M, alpha, r, Par, input, structure, mu=0.6, abundance='sola args_spectrum_irr=(), kwargs_spectrum_irr={}, cos_theta_irr=None, cos_theta_irr_exp=1 / 12, C_irr=None, T_irr=None, z0r_estimation=None, Sigma0_estimation=None, P_ph_0=None, verbose=False, - n=100, add_Pi_values=True, path_dots=None): + n=100, add_Pi_values=True, path_dots=None, + to_return=False, to_save=True, magn_args=None, + rad_pressure_include=True): """ Calculates vertical structure and makes table with disc parameters as functions of vertical coordinate. Table also contains input parameters of structure, parameters in the symmetry plane and parameter normalisations. @@ -413,9 +467,23 @@ def Vertical_Profile(M, alpha, r, Par, input, structure, mu=0.6, abundance='sola Whether to write Pi-parameters (see Ketsaris & Shakura, 1998) to the output file header. path_dots : str Where to save data table. + to_return : bool + Whether to return a namedtuple containing the calculated quantities and + debug data. Default is False. + to_save : bool + Whether to save the calculated data in a file. Default is True. + magn_args : dict or None, optional + The dictionary containing the parameters necessary for the magentic field + tretment. If None, the behaviour switches to a non-magnetic case. + Default is None. + rad_pressure_include : bool, optional + Whether to include the radiation pressure contribution in the equation + of state. Default is True. Returns ------- + -------------------------- if to_save == True: ---------------------------- + Table with calculated Vertical disc profile will save to path_dots. Table contains: 1) input parameters of the system -- M in Msun, alpha, r in cm and in rg, effective temperature Teff (viscous temperature Tvis in case of irradiation), accretion rate Mdot, viscous torque F, @@ -442,16 +510,27 @@ def Vertical_Profile(M, alpha, r, Par, input, structure, mu=0.6, abundance='sola In case of tabular EoS: 5) grad_ad, free_e -- adiabatic temperature gradient dlnT/dlnP and mean number of free electrons per nucleon. + + ------------------------- if to_return == True: --------------------------- + A namedtuple containing:[ + "t", "S", "P", 'T', "Q", "rho", "varkappa", "tau", "grad", + "r", "r_rg", "Sigma0", "Mdot", "Teff", "F", "z0r", + "rho_C", "T_C", "P_C", "structure", 'mu', + "PradPgas_C", "varkappa_C", + "tau_C", + "Pi1", "Pi2", "Pi3", "Pi4", + ] """ - if path_dots is None: + if path_dots is None and to_save: raise Exception("ATTENTION: the data wil not be saved, since 'path_dots' is None.") vs, F, Teff, Mdot = StructureChoice(M, alpha, r, Par, input, structure, mu, abundance, nu_irr=nu_irr, L_X_irr=L_X_irr, spectrum_irr=spectrum_irr, spectrum_irr_par=spectrum_irr_par, args_spectrum_irr=args_spectrum_irr, kwargs_spectrum_irr=kwargs_spectrum_irr, cos_theta_irr=cos_theta_irr, cos_theta_irr_exp=cos_theta_irr_exp, - C_irr=C_irr, T_irr=T_irr, F_in=F_in) + C_irr=C_irr, T_irr=T_irr, F_in=F_in, magn_args=magn_args, + rad_pressure_include=rad_pressure_include) if structure in ['MesaIrr', 'MesaRadAdIrr', 'MesaRadConvIrr']: kwargs_fit = {'z0r_estimation': z0r_estimation, 'Sigma0_estimation': Sigma0_estimation, 'verbose': verbose, 'P_ph_0': P_ph_0} @@ -512,10 +591,28 @@ def Vertical_Profile(M, alpha, r, Par, input, structure, mu=0.6, abundance='sola header_norm = f'\nSigma_norm = {vs.sigma_norm:e}, P_norm = {vs.P_norm:e}, ' \ f'T_norm = {vs.T_norm:e}, Q_norm = {vs.Q_norm:e}' if add_Pi_values: + Pi1, Pi2, Pi3, Pi4 = vs.Pi_finder() header_Pi = '\nPi1 = {:f}, Pi2 = {:f}, Pi3 = {:f}, Pi4 = {:f}'.format(*vs.Pi_finder()) header = header + header_input + header_C + header_norm + header_conv + header_Pi + header_input_irr - if path_dots is not None: + if path_dots is not None and to_save: np.savetxt(path_dots, dots_arr, header=header) + if to_return: + VertStructureResults = namedtuple("VertStructureResults", [ + "t", "S", "P", 'T', "Q", "rho", "varkappa", "tau", "grad", + "r", "r_rg", "Sigma0", "Mdot", "Teff", "F", "z0r", + "rho_C", "T_C", "P_C", "structure", 'mu', + "PradPgas_C", "varkappa_C", + "tau_C", + "Pi1", "Pi2", "Pi3", "Pi4", + ]) + _res = VertStructureResults(t=t, S=S, P=P, T=T, Q=Q, rho=rho, + varkappa=varkappa, tau=tau_arr, grad=grad_plot(np.log(P_full)), + r=r, r_rg=r/rg, Sigma0=Sigma0, + Mdot=Mdot, Teff=Teff, F=F, z0r=z0r, rho_C=rho_C, T_C=T_C, + P_C=P_C, structure=structure, mu=mu, PradPgas_C=delta, + varkappa_C=varkappa_C, tau_C=tau, + Pi1=Pi1, Pi2=Pi2, Pi3=Pi3, Pi4=Pi4) + return _res return @@ -523,7 +620,9 @@ def S_curve(Par_min, Par_max, M, alpha, r, input, structure, mu=0.6, abundance=' L_X_irr=None, spectrum_irr=None, spectrum_irr_par=None, args_spectrum_irr=(), kwargs_spectrum_irr={}, cos_theta_irr=None, cos_theta_irr_exp=1 / 12, C_irr=None, T_irr=None, z0r_start_estimation=None, Sigma0_start_estimation=None, P_ph_0=None, verbose=False, - n=100, tau_break=True, add_Pi_values=True, path_dots=None): + verbose_scurve=True, + n=100, tau_break=True, add_Pi_values=True, path_dots=None, + magn_args=None, rad_pressure_include=True, to_save=True, to_return=False): """ Calculates S-curve and makes table with disc parameters on the S-curve. Table contains input parameters of system, surface density Sigma0, viscous torque F, @@ -630,9 +729,27 @@ def S_curve(Par_min, Par_max, M, alpha, r, input, structure, mu=0.6, abundance=' Whether to write Pi-parameters (see Ketsaris & Shakura, 1998) to the output file. path_dots : str Where to save data table. - + verbose_scurve : bool + Whether to print info during the S-curve calculation. + Default is True. + to_return : bool + Whether to return a namedtuple containing the calculated quantities and + debug data. Default is False, and returned are: (Sigma_minus_index, + Sigma_minus_index). If True, returned are: (Sigma_minus_index, + Sigma_minus_index, SCurveResults). + to_save : bool + Whether to save the calculated data in a file. Default is True. + magn_args : dict or None, optional + The dictionary containing the parameters necessary for the magentic field + tretment. If None, the behaviour switches to a non-magnetic case. + Default is None. + rad_pressure_include : bool, optional + Whether to include the radiation pressure contribution in the equation + of state. Default is True. Returns - ------- + ------- + -------------------------- if to_save == True: ---------------------------- + Table with calculated S-curve will save to path_dots. Table contains: 1) input parameters of the system -- M in Msun, alpha, r in cm and in rg, structure type, mu (in case of analytical EoS) or abundance (in case of tabular EoS); @@ -655,9 +772,16 @@ def S_curve(Par_min, Par_max, M, alpha, r, input, structure, mu=0.6, abundance=' Also table contains Sigma_plus_index, Sigma_minus_index -- turn point indexes of the S-curve. Finally, table contains 'Non-converged_fits' -- number of unsuccessfully fitted structures. + + ------------------------- if to_return == True: --------------------------- + If to_return == False, returned are: (Sigma_minus_index, + Sigma_minus_index). If True, returned are: (Sigma_minus_index, + Sigma_minus_index, SCurveResults). SCurveResults contain: + [Sigma0, Mdot, Teff, F, z0r, T_C, + rho_C, P_C, tau, PradPgas_C, varkappa_C]. """ - if path_dots is None: + if path_dots is None and to_return: raise Exception("ATTENTION: the data wil not be saved, since 'path_dots' is None.") Sigma_minus_index = 0 @@ -697,17 +821,26 @@ def S_curve(Par_min, Par_max, M, alpha, r, input, structure, mu=0.6, abundance=' if structure in ['MesaIrr', 'MesaRadAdIrr', 'MesaRadConvIrr']: header += ' \tcost \tSigma_ph' header = header + '\nAll values are in CGS units.' + header_end - np.savetxt(path_dots, [], header=header) + if to_save: + np.savetxt(path_dots, [], header=header) + + SCurveResults = namedtuple("SCurveResults", [ + "Sigma0", "Mdot", "Teff", "F", "z0r", "T_C", + "rho_C", "P_C", "tau", "PradPgas_C", "varkappa_C", + ]) + acc = {k: [] for k in SCurveResults._fields} for i, Par in enumerate(np.geomspace(Par_max, Par_min, n)): - print(i) + if verbose_scurve: + print(i) vs, F, Teff, Mdot = StructureChoice(M, alpha, r, Par, input, structure, mu, abundance, F_in=F_in, nu_irr=nu_irr, L_X_irr=L_X_irr, spectrum_irr=spectrum_irr, spectrum_irr_par=spectrum_irr_par, args_spectrum_irr=args_spectrum_irr, kwargs_spectrum_irr=kwargs_spectrum_irr, cos_theta_irr=cos_theta_irr, cos_theta_irr_exp=cos_theta_irr_exp, - C_irr=C_irr, T_irr=T_irr) + C_irr=C_irr, T_irr=T_irr, magn_args=magn_args, + rad_pressure_include=rad_pressure_include) if structure in ['MesaIrr', 'MesaRadAdIrr', 'MesaRadConvIrr']: kwargs_fit = {'z0r_estimation': z0r_estimation, 'Sigma0_estimation': sigma_par_estimation, 'verbose': verbose, 'P_ph_0': P_ph_0} @@ -718,8 +851,9 @@ def S_curve(Par_min, Par_max, M, alpha, r, input, structure, mu=0.6, abundance=' try: result = vs.fit(**kwargs_fit) except Exception as e: - print(e) - print('Non-converged fit') + if verbose_scurve: + print(e) + print('Non-converged fit') Non_converged_fits += 1 continue try: @@ -734,19 +868,21 @@ def S_curve(Par_min, Par_max, M, alpha, r, input, structure, mu=0.6, abundance=' pass tau = vs.tau() - print(f'Mdot = {Mdot:1.3e} g/s, {Teff_string} = {Teff:g} K, tau = {tau:g}, z0r = {z0r:g}') + if verbose_scurve: + print(f'Mdot = {Mdot:1.3e} g/s, {Teff_string} = {Teff:g} K, tau = {tau:g}, z0r = {z0r:g}') if tau < 1 and tau_break: - print('Note: tau<1, tau_break=True. Cycle ends, when tau<1.') + if verbose_scurve: + print('Note: tau<1, tau_break=True. Cycle ends, when tau<1.') break varkappa_C, rho_C, T_C, P_C, Sigma0 = vs.parameters_C() PradPgas_C = (4 * sigmaSB) / (3 * c) * T_C ** 4 / P_C output_string = [Sigma0, Teff, Mdot, F, z0r, rho_C, T_C, P_C, tau, PradPgas_C, varkappa_C] - - print(f'Sigma0 = {Sigma0:g} g/cm^2') - print('Prad/Pgas_C = ', PradPgas_C) + if verbose_scurve: + print(f'Sigma0 = {Sigma0:g} g/cm^2') + print('Prad/Pgas_C = ', PradPgas_C) rho, eos = vs.law_of_rho(P_C, T_C, full_output=True) try: @@ -765,7 +901,8 @@ def S_curve(Par_min, Par_max, M, alpha, r, input, structure, mu=0.6, abundance=' QirrQvis = vs.Q_irr / vs.Q0 T_irr_, C_irr_ = vs.T_irr, vs.C_irr output_string.extend([QirrQvis, T_irr_, C_irr_]) - print(f'T_irr, C_irr = {T_irr_:g} K, {C_irr_:g}') + if verbose_scurve: + print(f'T_irr, C_irr = {T_irr_:g} K, {C_irr_:g}') if structure in ['MesaIrr', 'MesaRadAdIrr', 'MesaRadConvIrr']: output_string.extend([result.cost, vs.Sigma_ph]) @@ -792,20 +929,40 @@ def S_curve(Par_min, Par_max, M, alpha, r, input, structure, mu=0.6, abundance=' Sigma_plus_key = False output_string = np.array(output_string) + if to_save: + with open(path_dots, 'a') as file: + np.savetxt(file, output_string, newline=' ') + file.write('\n') + acc["Sigma0"].append(Sigma0) + acc["Teff"].append(Teff) + acc["z0r"].append(z0r) + acc["Mdot"].append(Mdot) + acc["F"].append(F) + acc["rho_C"].append(rho_C) + acc["P_C"].append(P_C) + acc["varkappa_C"].append(varkappa_C) + acc["tau"].append(tau) + acc["PradPgas_C"].append(PradPgas_C) + + if to_save: with open(path_dots, 'a') as file: - np.savetxt(file, output_string, newline=' ') - file.write('\n') - with open(path_dots, 'a') as file: - file.write(f'# Sigma_plus_index = {Sigma_plus_index:d} Sigma_minus_index = {Sigma_minus_index:d}') - file.write(f'\n# Non-converged_fits = {Non_converged_fits}') + file.write(f'# Sigma_plus_index = {Sigma_plus_index:d} Sigma_minus_index = {Sigma_minus_index:d}') + file.write(f'\n# Non-converged_fits = {Non_converged_fits}') + if to_return: + acc_np = {k: np.asarray(v, dtype=float) for k, v in acc.items()} + results_tuple = SCurveResults(**acc_np) + return Sigma_minus_index, Sigma_plus_index, results_tuple return Sigma_minus_index, Sigma_plus_index -def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abundance='solar', F_in=0, +def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abundance='solar', F_in=0, nu_irr=None, L_X_irr=None, spectrum_irr=None, spectrum_irr_par=None, args_spectrum_irr=(), kwargs_spectrum_irr={}, cos_theta_irr=None, cos_theta_irr_exp=1 / 12, C_irr=None, T_irr=None, z0r_start_estimation=None, Sigma0_start_estimation=None, P_ph_0=None, verbose=False, - n=100, tau_break=True, add_Pi_values=True, path_dots=None): + verbose_rad=False, to_return=False, to_save=True, + n=100, tau_break=True, add_Pi_values=True, path_dots=None, + magn_args=None, + rad_pressure_include=True): """ Calculates radial structure of disc. Return table, which contains input parameters of the system, surface density Sigma0, viscous torque F, accretion rate Mdot, effective temperature Teff @@ -916,6 +1073,23 @@ def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abun verbose : bool Whether to print values of free parameters at each iteration during fitting. Default is False, the fitting process performs silently. + verbose_rad : bool + Whether to print info during radial structure calculation at every radius. + Default is False, the radial calculation structure calculation process + performs silently. + to_return : bool + Whether to return a namedtuple containing the calculated quantities and + debug data. Default is False. + to_save : bool + Whether to save the calculated data in a file. Default is True. + magn_args : dict or None, optional + The dictionary containing the parameters necessary for the magentic field + tretment. If None, the behaviour switches to a non-magnetic case. + Default is None. + rad_pressure_include : bool, optional + Whether to include the radiation pressure contribution in the equation + of state. Default is True. + n : int Number of dots to calculate. tau_break : bool @@ -927,6 +1101,8 @@ def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abun Returns ------- + -------------------------- if to_save == True: ---------------------------- + Table with calculated Radial disc profile will save to path_dots. Table contains: 1) input parameters of the system -- M in Msun, alpha, structure type, mu (in case of analytical EoS) or abundance (in case of tabular EoS); @@ -947,9 +1123,16 @@ def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abun 8) cost, Sigma_ph -- cost function and column density of the layers above the photosphere. If structure is fitted successfully, cost must be less than 1e-16. Table also contains 'Non-converged_fits' -- number of unsuccessfully fitted structures. - + + ------------------------- if to_return == True: -------------------------- + RadStructureResults : namedtuple + Contains the calculated quantities as numpy arrays. The attributes + have the same names as described above, + e.g., the dependence of the Teff on radius may be plotted as + plt.plot(res.r, res.Teff). + """ - if path_dots is None: + if (path_dots is None) and to_save: raise Exception("ATTENTION: the data wil not be saved, since 'path_dots' is None.") z0r_estimation = z0r_start_estimation @@ -980,15 +1163,26 @@ def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abun if structure in ['MesaIrr', 'MesaRadAdIrr', 'MesaRadConvIrr']: header += ' \tcost \tSigma_ph' header = header + '\nAll values are in CGS units.' + header_end - np.savetxt(path_dots, [], header=header) + if to_save: + np.savetxt(path_dots, [], header=header) try: input_broadcast = np.broadcast(r_arr, Par, cos_theta_irr, cos_theta_irr_exp, C_irr, T_irr) except ValueError as e: raise ValueError(f'Array-like input parameters must have the same size n = {n}.') from e + RadStructureResults = namedtuple("RadStructureResults", [ + "r", "r_rg", "Sigma0", "Mdot", "Teff", "F", "z0r", + "rho_C", "T_C", "P_C", "tau", "PradPgas_C", "varkappa_C", + "free_e", "conv_param_z", "conv_param_sigma", + "Pi1", "Pi2", "Pi3", "Pi4", + "QirrQvis", "T_irr_", "C_irr_", + "result_cost", "vs_Sigma_ph", + ]) + acc = {k: [] for k in RadStructureResults._fields} for i, input_pars in enumerate(input_broadcast): - print(i) + if verbose_rad: + print(i) r = input_pars[0] vs, F, Teff, Mdot = StructureChoice(M=M, alpha=alpha, r=r, Par=input_pars[1], input=input, structure=structure, mu=mu, abundance=abundance, F_in=F_in, @@ -997,7 +1191,9 @@ def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abun args_spectrum_irr=args_spectrum_irr, kwargs_spectrum_irr=kwargs_spectrum_irr, cos_theta_irr=input_pars[2], cos_theta_irr_exp=input_pars[3], - C_irr=input_pars[4], T_irr=input_pars[5]) + C_irr=input_pars[4], T_irr=input_pars[5], + magn_args=magn_args, + rad_pressure_include=rad_pressure_include) if structure in ['MesaIrr', 'MesaRadAdIrr', 'MesaRadConvIrr']: kwargs_fit = {'z0r_estimation': z0r_estimation, 'Sigma0_estimation': sigma_par_estimation, 'verbose': verbose, 'P_ph_0': P_ph_0} @@ -1008,8 +1204,9 @@ def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abun try: result = vs.fit(**kwargs_fit) except Exception as e: - print(e) - print('Non-converged fit') + if verbose_rad: + print(e) + print('Non-converged fit') Non_converged_fits += 1 continue try: @@ -1025,21 +1222,25 @@ def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abun tau = vs.tau() rg = 2 * G * M / c ** 2 - print(f'r = {r:1.3e} cm = {r / rg:g} rg, Mdot = {Mdot:1.3e} g/s, ' + if verbose_rad: + print(f'r = {r:1.3e} cm = {r / rg:g} rg, Mdot = {Mdot:1.3e} g/s, ' f'{Teff_string} = {Teff:g} K, tau = {tau:g}, z0r = {z0r:g}') if tau < 1 and tau_break: - print('Note: tau<1, tau_break=True. Cycle ends, when tau<1.') + if verbose_rad: + print('Note: tau<1, tau_break=True. Cycle ends, when tau<1.') break varkappa_C, rho_C, T_C, P_C, Sigma0 = vs.parameters_C() PradPgas_C = (4 * sigmaSB) / (3 * c) * T_C ** 4 / P_C - output_string = [r, r / rg, Sigma0, Mdot, Teff, F, z0r, rho_C, T_C, P_C, tau, PradPgas_C, varkappa_C] - print(f'Sigma0 = {Sigma0:g} g/cm^2') - print('Prad/Pgas_C = ', PradPgas_C) + output_string = [r, r / rg, Sigma0, Mdot, Teff, F, z0r, rho_C, T_C, P_C, tau, PradPgas_C, varkappa_C] + if verbose_rad: + print(f'Sigma0 = {Sigma0:g} g/cm^2') + print('Prad/Pgas_C = ', PradPgas_C) rho, eos = vs.law_of_rho(P_C, T_C, full_output=True) + free_e, conv_param_z, conv_param_sigma = np.nan, np.nan, np.nan try: _ = eos.c_p free_e = np.exp(eos.lnfree_e) @@ -1047,26 +1248,61 @@ def Radial_Profile(M, alpha, r_start, r_end, Par, input, structure, mu=0.6, abun output_string.extend([free_e, conv_param_z, conv_param_sigma]) except AttributeError: pass - + Pi1, Pi2, Pi3, Pi4 = np.nan, np.nan, np.nan, np.nan if add_Pi_values: - output_string.extend(vs.Pi_finder()) - + Pi1, Pi2, Pi3, Pi4 = vs.Pi_finder() + output_string.extend([Pi1, Pi2, Pi3, Pi4]) + QirrQvis, T_irr_, C_irr_ = np.nan, np.nan, np.nan if structure in ['MesaIrr', 'MesaRadAdIrr', 'MesaRadConvIrr', 'MesaIrrZero', 'MesaRadAdIrrZero', 'MesaRadConvIrrZero']: QirrQvis = vs.Q_irr / vs.Q0 T_irr_, C_irr_ = vs.T_irr, vs.C_irr output_string.extend([QirrQvis, T_irr_, C_irr_]) - print(f'T_irr, C_irr = {T_irr_:g} K, {C_irr_:g}') - + if verbose_rad: + print(f'T_irr, C_irr = {T_irr_:g} K, {C_irr_:g}') + + result_cost, vs_Sigma_ph = np.nan, np.nan if structure in ['MesaIrr', 'MesaRadAdIrr', 'MesaRadConvIrr']: - output_string.extend([result.cost, vs.Sigma_ph]) + result_cost, vs_Sigma_ph = result.cost, vs.Sigma_ph + output_string.extend([result_cost, vs_Sigma_ph]) output_string = np.array(output_string) + if to_save: + with open(path_dots, 'a') as file: + np.savetxt(file, output_string, newline=' ') + file.write('\n') + acc["r"].append(r) + acc["r_rg"].append(r / rg) + acc["Sigma0"].append(Sigma0) + acc["Mdot"].append(Mdot) + acc["Teff"].append(Teff) + acc["F"].append(F) + acc["z0r"].append(z0r) + acc["rho_C"].append(rho_C) + acc["T_C"].append(T_C) + acc["P_C"].append(P_C) + acc["tau"].append(tau) + acc["PradPgas_C"].append(PradPgas_C) + acc["varkappa_C"].append(varkappa_C) + acc["free_e"].append(free_e) + acc["conv_param_z"].append(conv_param_z) + acc["conv_param_sigma"].append(conv_param_sigma) + acc["Pi1"].append(Pi1) + acc["Pi2"].append(Pi2) + acc["Pi3"].append(Pi3) + acc["Pi4"].append(Pi4) + acc["QirrQvis"].append(QirrQvis) + acc["T_irr_"].append(T_irr_) + acc["C_irr_"].append(C_irr_) + acc["result_cost"].append(result_cost) + acc["vs_Sigma_ph"].append(vs_Sigma_ph) + if to_save: with open(path_dots, 'a') as file: - np.savetxt(file, output_string, newline=' ') - file.write('\n') - with open(path_dots, 'a') as file: - file.write(f'# Non-converged_fits = {Non_converged_fits}') + file.write(f'# Non-converged_fits = {Non_converged_fits}') + acc_np = {k: np.asarray(v, dtype=float) for k, v in acc.items()} + results_tuple = RadStructureResults(**acc_np) + if to_return: + return results_tuple return diff --git a/alpha_disc/vs.py b/alpha_disc/vs.py index bcfa06e..93aa579 100755 --- a/alpha_disc/vs.py +++ b/alpha_disc/vs.py @@ -15,7 +15,7 @@ from astropy import constants as const from scipy.integrate import solve_ivp, simpson from scipy.optimize import brentq - +from alpha_disc.magn import MagneticField sigmaSB = const.sigma_sb.cgs.value R_gas = const.R.cgs.value G = const.G.cgs.value @@ -87,9 +87,12 @@ class BaseVerticalStructure: F : double Moment of viscosity forces in g*cm^2/s^2. eps : double, optional - Accuracy of vertical structure calculation. + Accuracy of vertical structure calculation. Default 1e-5. mu : double, optional - Molecular weight for ideal gas equation of state. + Molecular weight for ideal gas equation of state. Default 0.6. + rad_pressure_include : bool + Whether to consider the radiation pressure contribution to the total + pressure. Default True. Methods ------- @@ -106,13 +109,19 @@ class BaseVerticalStructure: """ - def __init__(self, Mx, alpha, r, F, eps=1e-5, mu=0.6): + def __init__(self, Mx, alpha, r, F, eps=1e-5, mu=0.6, + rad_pressure_include=True): self.mu = mu self.Mx = Mx self.GM = G * Mx self.alpha = alpha self.r = r self.F = F + if rad_pressure_include: + self.eta_rad = 1.0 + else: + self.eta_rad = 0.0 + self.omegaK = np.sqrt(self.GM / self.r ** 3) self.eps = eps @@ -132,34 +141,62 @@ def z0(self, z0): self.P_norm = (4 / 3) * self.Q_norm / (self.alpha * z0 * self.omegaK) self.T_norm = self.omegaK ** 2 * self.mu * z0 ** 2 / R_gas self.sigma_norm = 28 * self.Q_norm / (3 * self.alpha * z0 ** 2 * self.omegaK ** 3) + + def rad_pressure(self, T): + return 4 * sigmaSB / (3 * c) * T ** 4 * self.eta_rad def law_of_viscosity(self, P): return self.alpha * P def law_of_rho(self, P, T, full_output): + """If full_output==False, should return rho(P, T). + + Else, should return rho(P, T) and namedtuple containing: + dlnRho_dlnPgas_const_T, + dlnRho_dlnT_const_Pgas, + mu, + lnfree_e, + grad_ad + """ raise NotImplementedError def law_of_opacity(self, rho, T, lnfree_e): + """Should return kappa(rho, T)""" raise NotImplementedError def viscosity(self, y): return self.law_of_viscosity(y[Vars.P] * self.P_norm + - 4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4) + self.rad_pressure(y[Vars.T] * self.T_norm)) def rho(self, y, full_output): - return self.law_of_rho(y[Vars.P] * self.P_norm, y[Vars.T] * self.T_norm, full_output=full_output) + return self.law_of_rho(P=y[Vars.P] * self.P_norm, T=y[Vars.T] * self.T_norm, full_output=full_output) def opacity(self, y, lnfree_e): rho = self.rho(y, full_output=False) - return self.law_of_opacity(rho, y[Vars.T] * self.T_norm, lnfree_e=lnfree_e) + return self.law_of_opacity(rho=rho, T=y[Vars.T] * self.T_norm, lnfree_e=lnfree_e) + + def additional_photospheric_pressure(self): + """Additional pressure [cgs] to be added to the photospheric pressure + in CGS. + Zero in absence of the field; in case of field present it should be overwritten by mixins. + """ + return 0.0 + + def dP_magn_induced_dz(self, z): + """Returns d(P magn induced)/dz === d(b_induced_phi**2/8pi)/dz in CGS. + Used as an additional right-hand side for the dPtot/dz: + dPtot/dz = (dPtot/dz)_no_field - dP_magn_induced_dz. + Zero in absence of the field; in case of field present it should be overwritten by mixins. + """ + return z*0 - def photospheric_pressure_equation(self, tau, P): # P = P_total + def photospheric_pressure_equation(self, tau, P_total): # P = P_total T = self.Teff * (1 / 2 + 3 * tau / 4) ** (1 / 4) - P_rad = 4 * sigmaSB / (3 * c) * T ** 4 - if P - P_rad < 0 or np.isnan(P - P_rad): - raise PgasPradNotConvergeError(*(P - P_rad), P_rad=P_rad, t=0.0, z0r=self.z0 / self.r) - rho, eos = self.law_of_rho(P - P_rad, T, True) - varkappa = self.law_of_opacity(rho, T, lnfree_e=eos.lnfree_e) + P_rad = self.rad_pressure(T) + if P_total - P_rad < 0 or np.isnan(P_total - P_rad): + raise PgasPradNotConvergeError(*(P_total - P_rad), P_rad=P_rad, t=0.0, z0r=self.z0 / self.r) + rho, eos = self.law_of_rho(P = P_total - P_rad, T=T, full_output=True) + varkappa = self.law_of_opacity(rho=rho, T=T, lnfree_e=eos.lnfree_e) return self.z0 * self.omegaK ** 2 / varkappa def P_ph(self): @@ -167,10 +204,12 @@ def P_ph(self): solution = solve_ivp( self.photospheric_pressure_equation, [0, 2 / 3], - [1e-7 * self.P_norm + 4 * sigmaSB / (3 * c) * self.Teff ** 4 / 2], rtol=self.eps, method='RK23' + [1e-7 * self.P_norm + 0.5 * self.rad_pressure(self.Teff)], + rtol=self.eps, method='RK23' ) - P_rad = 4 * sigmaSB / (3 * c) * self.Teff ** 4 + P_rad = self.rad_pressure(self.Teff) result = solution.y[0][-1] - P_rad # P_gas = P_tot - P_rad + result += self.additional_photospheric_pressure() if result < 0 or np.isnan(result): raise PgasPradNotConvergeError(P_gas=result, P_rad=P_rad, t=0.0, z0r=self.z0 / self.r) return result @@ -197,6 +236,7 @@ def initial(self): return y def dlnTdlnP(self, y, t): + """Returns d(lnT)/d(lnPtotal) """ raise NotImplementedError def dQdz(self, y, t): @@ -205,12 +245,14 @@ def dQdz(self, y, t): def dydt(self, t, y): """ - The right side of ODEs system. - + The right side of ODEs system. Magnetic fields modify this system only + as an additional source in a hydrostatic equation: + dPtot/dz = -rho wk^2 z - dPmagn/dz, + where dPmagn/dz is a known function of z. Parameters ---------- t : array-like - Modified vertical coordinate (t = 1 - z/z0). + Modified vertical coordinate (t = 1 - z/z0, so z = z0(1-t)). y : Current values of (dimensionless) unknown functions. @@ -223,24 +265,25 @@ def dydt(self, t, y): if y[Vars.P] < 0 or np.isnan(y[Vars.P]): try: raise PgasPradNotConvergeError(P_gas=y[Vars.P] * self.P_norm, - P_rad=4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4, + P_rad=self.rad_pressure(y[Vars.T] * self.T_norm), t=t, z0r=self.z0 / self.r, Sigma0_par=self.Sigma0_par) except AttributeError: raise PgasPradNotConvergeError(P_gas=y[Vars.P] * self.P_norm, - P_rad=4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4, + P_rad=self.rad_pressure(y[Vars.T] * self.T_norm), t=t, z0r=self.z0 / self.r) from None rho, eos = self.rho(y, full_output=True) - A = rho * (1 - t) * self.omegaK ** 2 * self.z0 ** 2 / self.P_norm - B = 16 * sigmaSB / (3 * c) * self.T_norm ** 4 * y[Vars.T] ** 3 / self.P_norm - - P_full = y[Vars.P] * self.P_norm + 4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4 - dP_full = A * self.P_norm + dp_full_dt = rho * (1 - t) * self.omegaK ** 2 * self.z0 ** 2 / self.P_norm # d (P total dimentionless)/dt + dp_full_dt += self.dP_magn_induced_dz(z=self.z0 * (1 - t)) * self.z0 / self.P_norm + B = 4 * self.rad_pressure(self.T_norm * y[Vars.T]) / y[Vars.T] / self.P_norm + P_full = y[Vars.P] * self.P_norm + self.rad_pressure(y[Vars.T] * self.T_norm) + + dP_full_dt = dp_full_dt * self.P_norm # dP_full/dt - grad = self.dlnTdlnP(y, t) - dTdz = grad * dP_full * y[Vars.T] / P_full + grad = self.dlnTdlnP(y, t) # d(ln T) / d(ln P_full) === (P_full/T) * (dT/dP_full) + dTdz = grad * dP_full_dt * y[Vars.T] / P_full dy[Vars.S] = 2 * rho * self.z0 / self.sigma_norm - dy[Vars.P] = A - B * dTdz + dy[Vars.P] = dp_full_dt - B * dTdz dy[Vars.Q] = self.dQdz(y, t) dy[Vars.T] = dTdz return dy @@ -326,10 +369,10 @@ def Pi_finder(self): """ varkappa_C, rho_C, T_C, P_C, Sigma0 = self.parameters_C() - - Pi_1 = (self.omegaK ** 2 * self.z0 ** 2 * rho_C) / (P_C + 4 * sigmaSB / (3 * c) * T_C ** 4) + P_C_tot = P_C + self.rad_pressure(T_C) + Pi_1 = (self.omegaK ** 2 * self.z0 ** 2 * rho_C) / P_C_tot Pi_2 = Sigma0 / (2 * self.z0 * rho_C) - Pi_3 = (3 / 4) * (self.alpha * self.omegaK * (P_C + 4 * sigmaSB / (3 * c) * T_C ** 4) * Sigma0) / ( + Pi_3 = (3 / 4) * (self.alpha * self.omegaK * P_C_tot * Sigma0) / ( self.Q0 * rho_C) Pi_4 = (3 / 32) * (self.Teff / T_C) ** 4 * (Sigma0 * varkappa_C) @@ -401,13 +444,13 @@ def dlnTdlnP(self, y, t): if t == 1: dTdz_der = (self.dQdz(y, t) / y[Vars.T] ** 3) * 3 * varkappa * rho * self.z0 * self.Q_norm / ( 16 * sigmaSB * self.T_norm ** 4) - P_full = y[Vars.P] * self.P_norm + 4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4 + P_full = y[Vars.P] * self.P_norm + self.rad_pressure(y[Vars.T]*self.T_norm) dP_full_der = - rho * self.omegaK ** 2 * self.z0 ** 2 dlnTdlnP_rad = (P_full / y[Vars.T]) * (dTdz_der / dP_full_der) else: dTdz = (abs(y[Vars.Q]) / y[Vars.T] ** 3) * 3 * varkappa * rho * self.z0 * self.Q_norm / ( 16 * sigmaSB * self.T_norm ** 4) - P_full = y[Vars.P] * self.P_norm + 4 * sigmaSB / (3 * c) * y[Vars.T] ** 4 * self.T_norm ** 4 + P_full = y[Vars.P] * self.P_norm + self.rad_pressure(y[Vars.T]*self.T_norm) dP_full = rho * (1 - t) * self.omegaK ** 2 * self.z0 ** 2 dlnTdlnP_rad = (P_full / y[Vars.T]) * (dTdz / dP_full) return dlnTdlnP_rad @@ -469,6 +512,28 @@ class IdealBellLin1994VerticalStructure(IdealGasMixin, BellLin1994TwoComponentOp """ pass +class IdealKramersVerticalStructureMagnetic(MagneticField, IdealGasMixin, KramersOpacityMixin, + RadiativeTempGradient, BaseVerticalStructure): + """ + Vertical structure class for Kramers opacity law and ideal gas EOS. + Includes the magnetic field. + + """ + pass + + +class IdealBellLin1994VerticalStructureMagnetic(MagneticField, IdealGasMixin, + BellLin1994TwoComponentOpacityMixin, RadiativeTempGradient, + BaseVerticalStructure): + """ + Vertical structure class for opacity laws from (Bell & Lin, 1994) and ideal gas EOS. + Includes the magnetic field. + + """ + pass + + + def main(): M = 10 * M_sun diff --git a/pyproject.toml b/pyproject.toml old mode 100644 new mode 100755 index 951500f..dd4af12 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "alpha-disc" -version = "1.0" +version = "1.1" authors = [ {email = "tavleev.as15@physics.msu.ru"}, {name = "Andrey Tavleev"} diff --git a/requirements.txt b/requirements.txt old mode 100644 new mode 100755 diff --git a/setup.py b/setup.py old mode 100644 new mode 100755 diff --git a/tests/test_data/S-curve.dat b/tests/test_data/S-curve.dat old mode 100644 new mode 100755 diff --git a/tests/test_data/S-curve_mesa.dat b/tests/test_data/S-curve_mesa.dat old mode 100644 new mode 100755 diff --git a/tests/test_data/radial_struct.dat b/tests/test_data/radial_struct.dat old mode 100644 new mode 100755 diff --git a/tests/test_data/radial_struct_mesa.dat b/tests/test_data/radial_struct_mesa.dat old mode 100644 new mode 100755 diff --git a/tests/test_data/vs.dat b/tests/test_data/vs.dat old mode 100644 new mode 100755 diff --git a/tests/test_data/vs_mesa.dat b/tests/test_data/vs_mesa.dat old mode 100644 new mode 100755 diff --git a/tests/test_mesa_vs.py b/tests/test_mesa_vs.py old mode 100644 new mode 100755 diff --git a/tests/test_vs.py b/tests/test_vs.py old mode 100644 new mode 100755