<commit_before># coding: utf-8

""" MW Potential used in Law & Majewski 2010 """

from __future__ import division, print_function

__author__ = "adrn <<EMAIL>>"

# Standard library
import os, sys
import math

# Third-party
import numpy as np
import astropy.units as u
from astropy.constants import G
from scipy.signal import argrelmin,argrelmax

from .core import CartesianPotential, CompositePotential, PotentialParameter
from .common import MiyamotoNagaiPotential, HernquistPotential, LogarithmicPotentialLJ
from ._lm10_acceleration import lm10_acceleration
from .. import usys
from ..util import _parse_quantity

from ..inference.parameter import ModelParameter
from ..inference.prior import LogUniformPrior

class LawMajewski2010(CompositePotential):

    def __init__(self, **parameters):
        """ Represents the functional form of the Galaxy potential used by
            Law and Majewski 2010.

            Miyamoto-Nagai disk
            Hernquist bulge
            Logarithmic halo

            Model parameters: q1, qz, phi, v_halo

            Parameters
            ----------
            parameters : dict
                A dictionary of parameters for the potential definition.
        """

        params = dict()
        params['q1'] = ModelParameter(name='q1', truth=1.38,
                            prior=LogUniformPrior(1.0,1.7))
        params['q2'] = ModelParameter(name='q2', truth=1.,
                            prior=LogUniformPrior(0.71,1.5))
        params['qz'] = ModelParameter(name='qz', truth=1.36,
                            prior=LogUniformPrior(1.0,1.7))
        params['phi'] = ModelParameter(name='phi',
                            truth=(97.*u.deg).decompose(usys).value,
                            prior=LogUniformPrior((80.*u.deg).decompose(usys).value,
                                                  (120.*u.deg).decompose(usys).value))
        params['v_halo'] = ModelParameter(name='v_halo',
                            truth=(121.858*u.km/u.s).decompose(usys).value,
                            prior=LogUniformPrior((100.*u.km/u.s).decompose(usys).value,
                                                  (150.*u.km/u.s).decompose(usys).value))
        params['R_halo'] = ModelParameter(name='R_halo',
                            truth=(12.*u.kpc).decompose(usys).value,
                            prior=LogUniformPrior((8.*u.kpc).decompose(usys).value,
                                                  (20*u.kpc).decompose(usys).value))

        self._parameter_dict = dict()
        self.parameters = dict(params)
        for k,v in self.parameters.items():
            #self.__dict__[k] = v
            self._parameter_dict[k] = parameters.get(k, params[k].truth)

        ####################
        # bulge = HernquistPotential(usys,
        #                            m=3.4E10*u.M_sun,
        #                            c=0.7*u.kpc)

        # disk = MiyamotoNagaiPotential(usys,
        #                               m=1.E11*u.M_sun,
        #                               a=6.5*u.kpc,
        #                               b=0.26*u.kpc)

        # p_dict = dict([(k,v) for k,v in self._parameter_dict.items()])
        # halo = LogarithmicPotentialLJ(usys, **p_dict)

        # super(LawMajewski2010, self).__init__(usys,
        #                                       bulge=bulge,
        #                                       disk=disk,
        #                                       halo=halo)
        ####################
        self.units = usys
        self._G = G.decompose(bases=usys).value

    def _acceleration_at(self, r, n_particles, acc):
        return lm10_acceleration(r, n_particles, acc, **self._parameter_dict)

    def _enclosed_mass(self, R):
        """ Compute the enclosed mass at the position r. Assumes it's far from
            the disk and bulge, and is a spherically averaged approximation.
        """

        #m_halo_enc = self["halo"]._parameters["v_halo"]**2 * R/self._G
        #m_enc = self["disk"]._parameters["m"] + \
        #        self["bulge"]._parameters["m"] + \
        #        m_halo_enc

        # TODO: HACK!!!
        Rh = self._parameter_dict['R_halo']
        vh = self._parameter_dict['v_halo']
        G = self._G
        m_halo_enc = (2*R**3*vh**2) / (G*(R**2 + Rh**2))
        m_enc = 1.E11 + 3.4E10 + m_halo_enc

        return m_enc

    def _tidal_radius(self, m, r):
        """ Compute the tidal radius of a massive particle at the specified
            position(s). Assumes position and mass are in the same unit
            system as the potential.

            Parameters
            ----------
            m : numeric
                Mass.
            r : array_like
                Position.
        """

        # Radius of Sgr center relative to galactic center
        R = np.sqrt(np.sum(r**2., axis=-1))
        m_enc = self._enclosed_mass(R)
        Rh = self._parameter_dict['R_halo']

        dlnM_dlnR = (3*Rh**2 + R**2)/(Rh**2 + R**2)
        f = (1 - dlnM_dlnR/3.)**(-0.3333333333333)

        return R * (m / m_enc)**(0.3333333333333) * f * 1.4

    def _read_tidal_radius(self, m, r, alpha=1.):
        """ Compute the tidal radius of a massive particle at the specified
            position(s). Assumes position and mass are in the same unit
            system as the potential.

            Parameters
            ----------
            m : numeric
                Mass.
            r : array_like
                Position.
            alpha : numeric
                1 for prograde, 0 for radial, -1 for retrograde orbits.
        """

        # Radius of Sgr center relative to galactic center
        R = np.sqrt(np.sum(r**2., axis=-1))
<commit_msg>define as classical tidal radius<commit_after># coding: utf-8

""" MW Potential used in Law & Majewski 2010 """

from __future__ import division, print_function

__author__ = "adrn <<EMAIL>>"

# Standard library
import os, sys
import math

# Third-party
import numpy as np
import astropy.units as u
from astropy.constants import G
from scipy.signal import argrelmin,argrelmax

from .core import CartesianPotential, CompositePotential, PotentialParameter
from .common import MiyamotoNagaiPotential, HernquistPotential, LogarithmicPotentialLJ
from ._lm10_acceleration import lm10_acceleration
from .. import usys
from ..util import _parse_quantity

from ..inference.parameter import ModelParameter
from ..inference.prior import LogUniformPrior

class LawMajewski2010(CompositePotential):

    def __init__(self, **parameters):
        """ Represents the functional form of the Galaxy potential used by
            Law and Majewski 2010.

            Miyamoto-Nagai disk
            Hernquist bulge
            Logarithmic halo

            Model parameters: q1, qz, phi, v_halo

            Parameters
            ----------
            parameters : dict
                A dictionary of parameters for the potential definition.
        """

        params = dict()
        params['q1'] = ModelParameter(name='q1', truth=1.38,
                            prior=LogUniformPrior(1.0,1.7))
        params['q2'] = ModelParameter(name='q2', truth=1.,
                            prior=LogUniformPrior(0.71,1.5))
        params['qz'] = ModelParameter(name='qz', truth=1.36,
                            prior=LogUniformPrior(1.0,1.7))
        params['phi'] = ModelParameter(name='phi',
                            truth=(97.*u.deg).decompose(usys).value,
                            prior=LogUniformPrior((80.*u.deg).decompose(usys).value,
                                                  (120.*u.deg).decompose(usys).value))
        params['v_halo'] = ModelParameter(name='v_halo',
                            truth=(121.858*u.km/u.s).decompose(usys).value,
                            prior=LogUniformPrior((100.*u.km/u.s).decompose(usys).value,
                                                  (150.*u.km/u.s).decompose(usys).value))
        params['R_halo'] = ModelParameter(name='R_halo',
                            truth=(12.*u.kpc).decompose(usys).value,
                            prior=LogUniformPrior((8.*u.kpc).decompose(usys).value,
                                                  (20*u.kpc).decompose(usys).value))

        self._parameter_dict = dict()
        self.parameters = dict(params)
        for k,v in self.parameters.items():
            #self.__dict__[k] = v
            self._parameter_dict[k] = parameters.get(k, params[k].truth)

        ####################
        # bulge = HernquistPotential(usys,
        #                            m=3.4E10*u.M_sun,
        #                            c=0.7*u.kpc)

        # disk = MiyamotoNagaiPotential(usys,
        #                               m=1.E11*u.M_sun,
        #                               a=6.5*u.kpc,
        #                               b=0.26*u.kpc)

        # p_dict = dict([(k,v) for k,v in self._parameter_dict.items()])
        # halo = LogarithmicPotentialLJ(usys, **p_dict)

        # super(LawMajewski2010, self).__init__(usys,
        #                                       bulge=bulge,
        #                                       disk=disk,
        #                                       halo=halo)
        ####################
        self.units = usys
        self._G = G.decompose(bases=usys).value

    def _acceleration_at(self, r, n_particles, acc):
        return lm10_acceleration(r, n_particles, acc, **self._parameter_dict)

    def _enclosed_mass(self, R):
        """ Compute the enclosed mass at the position r. Assumes it's far from
            the disk and bulge, and is a spherically averaged approximation.
        """

        #m_halo_enc = self["halo"]._parameters["v_halo"]**2 * R/self._G
        #m_enc = self["disk"]._parameters["m"] + \
        #        self["bulge"]._parameters["m"] + \
        #        m_halo_enc

        # TODO: HACK!!!
        Rh = self._parameter_dict['R_halo']
        vh = self._parameter_dict['v_halo']
        G = self._G
        m_halo_enc = (2*R**3*vh**2) / (G*(R**2 + Rh**2))
        m_enc = 1.E11 + 3.4E10 + m_halo_enc

        return m_enc

    def _tidal_radius(self, m, r):
        """ Compute the tidal radius of a massive particle at the specified
            position(s). Assumes position and mass are in the same unit
            system as the potential.

            Parameters
            ----------
            m : numeric
                Mass.
            r : array_like
                Position.
        """

        # Radius of Sgr center relative to galactic center
        R = np.sqrt(np.sum(r**2., axis=-1))
        m_enc = self._enclosed_mass(R)
        Rh = self._parameter_dict['R_halo']

        dlnM_dlnR = (3*Rh**2 + R**2)/(Rh**2 + R**2)
        f = (1 - dlnM_dlnR/3.)**(-0.3333333333333)

        return R * (m / (3*m_enc))**(0.3333333333333) * f

    def _read_tidal_radius(self, m, r, alpha=1.):
        """ Compute the tidal radius of a massive particle at the specified
            position(s). Assumes position and mass are in the same unit
            system as the potential.

            Parameters
            ----------
            m : numeric
                Mass.
            r : array_like
                Position.
            alpha : numeric
                1 for prograde, 0 for radial, -1 for retrograde orbits.
        """

        # Radius of Sgr center relative to galactic center
        R = np.sqrt(np.sum(r**2., axis=-1))

{'fraction_non_alphanumeric': 0.12103469369920475, 'fraction_numerical': 0.02508083544525037, 'mean_word_length': 1.6151736745886653, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '88563', 'n_tokens_mistral': 3631, 'n_tokens_neox': 3315, 'n_words': 739}