<commit_before># -*- coding: utf-8 -*-
"""
Created on Thu Jul 16 15:42:56 2015

@author: ibackus
"""
import warnings
import pynbody as pb
SimArray = pb.array.SimArray
import numpy as np

from diskpy.utils import match_units
from diskpy.pdmath import binned_mean

def Q(snapshot, molecular_mass = 2.0, bins=100, use_velocity=False, \
use_omega=True):
    """Calculates the Toomre Q as a function of r, assuming radial temperature
    profile and kappa ~= omega
    
    Parameters
    ----------
    
    snapshot : tipsy snapshot
    molecular_mass : float
        Mean molecular mass (for sound speed).  Default = 2.0
    bins : int or array
        Either the number of bins or the bin edges
    use_velocity : Bool
        Determines whether to use the particles' velocities to calculate orbital
        velocity.  Useful if the circular orbital velocities are set in the
        snapshot.
    use_omega : Bool
        Default=True.  Use omega as a proxy for kappa to reduce noise

    Returns
    -------
    
    Q : array
        Toomre Q as a function of r
    r_edges : array
        Radial bin edges
    """

    # Physical constants
    kB = SimArray([1.0],'k')
    G = SimArray([1.0],'G')
    # Calculate surface density
    sig, r_edges = sigma(snapshot, bins)
    # Calculate sound speed
    m = match_units(molecular_mass,'m_p')[0]
    c_s_all = np.sqrt(kB*snapshot.g['temp']/m)
    # Bin/average sound speed
    dummy, c_s, dummy2 = binned_mean(snapshot.g['rxy'], c_s_all, binedges=r_edges)

    if use_omega:
        # Calculate keplerian angular velocity (as a proxy for the epicyclic
        # frequency, which is a noisy calculation)
        if use_velocity:
            # Calculate directly from particle's velocity
            dummy, omega, dummy2 = binned_mean(snapshot.g['rxy'], \
            snapshot.g['vt']/snapshot.g['rxy'], binedges=r_edges)

        else:
            # Estimate, from forces, using pynbody
            p = pb.analysis.profile.Profile(snapshot, bins=r_edges)
            omega = p['omega']

        kappa_calc = omega

    else:

        if use_velocity:
            # Calculate directly from particle's velocities
            kappa_calc, dummy = kappa(snapshot, r_edges)

        else:
            # Estimate, from forces, using pynbody
            p = pb.analysis.profile.Profile(snapshot, bins=r_edges)
            kappa_calc = p['kappa']

    return (kappa_calc*c_s/(np.pi*G*sig)).in_units('1'), r_edges

def Qeff(snapshot, molecular_mass = 2.0, bins=100, use_velocity=False, \
use_omega=True, alpha=0.18, beta=2.2):
    """Estimates the effective Toomre Q as a function of r, defined as:
    
    .. math:: Q_{eff} = \\beta Q (h/R)^{\\alpha}
        
    See Q and h for the estimates of Q and h
    
    Parameters
    ----------
    
    snapshot : tipsy snapshot
    molecular_mass : float
        Mean molecular mass (for sound speed).  Default = 2.0
    bins : int or array<commit_msg>Updated (slightly) default alpha & beta values for Q_eff<commit_after># -*- coding: utf-8 -*-
"""
Created on Thu Jul 16 15:42:56 2015

@author: ibackus
"""
import warnings
import pynbody as pb
SimArray = pb.array.SimArray
import numpy as np

from diskpy.utils import match_units
from diskpy.pdmath import binned_mean

def Q(snapshot, molecular_mass = 2.0, bins=100, use_velocity=False, \
use_omega=True):
    """Calculates the Toomre Q as a function of r, assuming radial temperature
    profile and kappa ~= omega
    
    Parameters
    ----------
    
    snapshot : tipsy snapshot
    molecular_mass : float
        Mean molecular mass (for sound speed).  Default = 2.0
    bins : int or array
        Either the number of bins or the bin edges
    use_velocity : Bool
        Determines whether to use the particles' velocities to calculate orbital
        velocity.  Useful if the circular orbital velocities are set in the
        snapshot.
    use_omega : Bool
        Default=True.  Use omega as a proxy for kappa to reduce noise

    Returns
    -------
    
    Q : array
        Toomre Q as a function of r
    r_edges : array
        Radial bin edges
    """

    # Physical constants
    kB = SimArray([1.0],'k')
    G = SimArray([1.0],'G')
    # Calculate surface density
    sig, r_edges = sigma(snapshot, bins)
    # Calculate sound speed
    m = match_units(molecular_mass,'m_p')[0]
    c_s_all = np.sqrt(kB*snapshot.g['temp']/m)
    # Bin/average sound speed
    dummy, c_s, dummy2 = binned_mean(snapshot.g['rxy'], c_s_all, binedges=r_edges)

    if use_omega:
        # Calculate keplerian angular velocity (as a proxy for the epicyclic
        # frequency, which is a noisy calculation)
        if use_velocity:
            # Calculate directly from particle's velocity
            dummy, omega, dummy2 = binned_mean(snapshot.g['rxy'], \
            snapshot.g['vt']/snapshot.g['rxy'], binedges=r_edges)

        else:
            # Estimate, from forces, using pynbody
            p = pb.analysis.profile.Profile(snapshot, bins=r_edges)
            omega = p['omega']

        kappa_calc = omega

    else:

        if use_velocity:
            # Calculate directly from particle's velocities
            kappa_calc, dummy = kappa(snapshot, r_edges)

        else:
            # Estimate, from forces, using pynbody
            p = pb.analysis.profile.Profile(snapshot, bins=r_edges)
            kappa_calc = p['kappa']

    return (kappa_calc*c_s/(np.pi*G*sig)).in_units('1'), r_edges

def Qeff(snapshot, molecular_mass = 2.0, bins=100, use_velocity=False, \
use_omega=True, alpha=0.1766, beta=2.121):
    """Estimates the effective Toomre Q as a function of r, defined as:
    
    .. math:: Q_{eff} = \\beta Q (h/R)^{\\alpha}
        
    See Q and h for the estimates of Q and h.  
    
    The values for alpha and beta are from Backus & Quinn, 2016 and are 
    appropriate for an isothermal disk.
    
    Parameters
    ----------
    
    snapshot : tipsy snapshot
    molecular_mass : float
        Mean molecular mass (for sound speed).  Default = 2.0
    bins : int or array
{'fraction_non_alphanumeric': 0.09871031746031746, 'fraction_numerical': 0.01455026455026455, 'mean_word_length': 2.3793296089385474, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '737844', 'n_tokens_mistral': 1950, 'n_tokens_neox': 1829, 'n_words': 600}