<commit_before>import math as m
import numpy as nu
from scipy import integrate
from galpy import actionAngle
from galpy.potential import LogarithmicHaloPotential, PowerSphericalPotential,\
    KeplerPotential
from galpy.potential_src.Potential import evaluateRforces, evaluatezforces,\
    evaluatePotentials, evaluatephiforces, evaluateDensities
import galpy.util.bovy_plot as plot
import galpy.util.bovy_symplecticode as symplecticode
from OrbitTop import OrbitTop
class FullOrbit(OrbitTop):
    """Class that holds and integrates orbits in full 3D potentials"""
    def __init__(self,vxvv=[1.,0.,0.9,0.,0.1]):
        """
        NAME:
           __init__
        PURPOSE:
           intialize a full orbit
        INPUT:
           vxvv - initial condition [R,vR,vT,z,vz,phi]
        OUTPUT:
           (none)
        HISTORY:
           2010-08-01 - Written - Bovy (NYU)
        """
        self.vxvv= vxvv
        #For boundary-condition integration
        self._BCIntegrateFunction= _integrateFullOrbit
        return None

    def integrate(self,t,pot,method='odeint'):
        """
        NAME:
           integrate
        PURPOSE:
           integrate the orbit
        INPUT:
           t - list of times at which to output (0 has to be in this!)
           pot - potential instance or list of instances
           method= 'odeint' for scipy's odeint integration, 'leapfrog' for
                    a simple symplectic integrator
        OUTPUT:
           (none) (get the actual orbit using getOrbit()
        HISTORY:
           2010-08-01 - Written - Bovy (NYU)
        """
        if method == 'leapfrog_c': method= 'odeint'
        #Reset things that may have been defined by a previous integration
        if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp')
        if hasattr(self,'rs'): delattr(self,'rs')
        self.t= nu.array(t)
        self._pot= pot
        self.orbit= _integrateFullOrbit(self.vxvv,pot,t,method)

    def Jacobi(self,Omega,t=0.,pot=None):
        """
        NAME:
           Jacobi
        PURPOSE:
           calculate the Jacobi integral of the motion
        INPUT:
           Omega - pattern speed of rotating frame
           t= time
           pot= potential instance or list of such instances
        OUTPUT:
           Jacobi integral
        HISTORY:
           2011-04-18 - Written - Bovy (NYU)
        """
        if len(Omega) == 3:
            if isinstance(Omega,list): thisOmega= nu.array(Omega)
            else: thisOmega= Omega
            return self.E(pot=pot,t=t)-nu.dot(thisOmega,self.L())
        else:
            return self.E(pot=pot,t=t)-Omega*self.L()[2]

    def E(self,pot=None,t=0.):
        """
        NAME:
           E
        PURPOSE:
           calculate the energy
        INPUT:
           pot= potential instance or list of such instances
           t= time at which to evaluate L and t
        OUTPUT:
           energy
        HISTORY:
           2010-09-15 - Written - Bovy (NYU)
        """
        if pot is None:
            try:
                pot= self._pot
            except AttributeError:
                raise AttributeError("Integrate orbit or specify pot=")
        return evaluatePotentials(self.vxvv[0],self.vxvv[3],pot,
                                  phi=self.vxvv[5],t=t)+\
                                  self.vxvv[1]**2./2.+self.vxvv[2]**2./2.+\
                                  self.vxvv[4]**2./2.

    def e(self,analytic=False,pot=None):
        """
        NAME:
           e
        PURPOSE:
           calculate the eccentricity
        INPUT:
           analytic - compute this analytically
           pot - potential to use for analytical calculation
        OUTPUT:
           eccentricity
        HISTORY:
           2010-09-15 - Written - Bovy (NYU)
        """
        if analytic:
            if not hasattr(self,'_aA'):
                self._setupaA(pot=pot)
            (rperi,rap)= self._aA.calcRapRperi()
            return (rap-rperi)/(rap+rperi)
        if not hasattr(self,'orbit'):
            raise AttributeError("Integrate the orbit first")
        if not hasattr(self,'rs'):
            self.rs= nu.sqrt(self.orbit[:,0]**2.+self.orbit[:,3]**2.)
        return (nu.amax(self.rs)-nu.amin(self.rs))/(nu.amax(self.rs)+nu.amin(self.rs))

    def rap(self,analytic=False,pot=None):<commit_msg>change and fix E() call<commit_after>import math as m
import numpy as nu
from scipy import integrate
from galpy import actionAngle
from galpy.potential import LogarithmicHaloPotential, PowerSphericalPotential,\
    KeplerPotential
from galpy.potential_src.Potential import evaluateRforces, evaluatezforces,\
    evaluatePotentials, evaluatephiforces, evaluateDensities
import galpy.util.bovy_plot as plot
import galpy.util.bovy_symplecticode as symplecticode
from OrbitTop import OrbitTop
class FullOrbit(OrbitTop):
    """Class that holds and integrates orbits in full 3D potentials"""
    def __init__(self,vxvv=[1.,0.,0.9,0.,0.1]):
        """
        NAME:
           __init__
        PURPOSE:
           intialize a full orbit
        INPUT:
           vxvv - initial condition [R,vR,vT,z,vz,phi]
        OUTPUT:
           (none)
        HISTORY:
           2010-08-01 - Written - Bovy (NYU)
        """
        self.vxvv= vxvv
        #For boundary-condition integration
        self._BCIntegrateFunction= _integrateFullOrbit
        return None

    def integrate(self,t,pot,method='odeint'):
        """
        NAME:
           integrate
        PURPOSE:
           integrate the orbit
        INPUT:
           t - list of times at which to output (0 has to be in this!)
           pot - potential instance or list of instances
           method= 'odeint' for scipy's odeint integration, 'leapfrog' for
                    a simple symplectic integrator
        OUTPUT:
           (none) (get the actual orbit using getOrbit()
        HISTORY:
           2010-08-01 - Written - Bovy (NYU)
        """
        if method == 'leapfrog_c': method= 'odeint'
        #Reset things that may have been defined by a previous integration
        if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp')
        if hasattr(self,'rs'): delattr(self,'rs')
        self.t= nu.array(t)
        self._pot= pot
        self.orbit= _integrateFullOrbit(self.vxvv,pot,t,method)

    def Jacobi(self,Omega,t=0.,pot=None):
        """
        NAME:
           Jacobi
        PURPOSE:
           calculate the Jacobi integral of the motion
        INPUT:
           Omega - pattern speed of rotating frame
           t= time
           pot= potential instance or list of such instances
        OUTPUT:
           Jacobi integral
        HISTORY:
           2011-04-18 - Written - Bovy (NYU)
        """
        if len(Omega) == 3:
            if isinstance(Omega,list): thisOmega= nu.array(Omega)
            else: thisOmega= Omega
            return self.E(pot=pot,t=t)-nu.dot(thisOmega,self.L())
        else:
            return self.E(pot=pot,t=t)-Omega*self.L()[2]

    def E(self,*args,**kwargs):
        """
        NAME:
           E
        PURPOSE:
           calculate the energy
        INPUT:
           t - (optional) time at which to get the radius
           pot= potential instance or list of such instances
        OUTPUT:
           energy
        HISTORY:
           2010-09-15 - Written - Bovy (NYU)
        """
        if not kwargs.has_key('pot'):
            try:
                pot= self._pot
            except AttributeError:
                raise AttributeError("Integrate orbit or specify pot=")
        else:
            pot= kwargs['pot']
            kwargs.pop('pot')
        #Get orbit
        thiso= self(*args,**kwargs)
        onet= (len(thiso.shape) == 1)
        if onet:
            return evaluatePotentials(thiso[0],thiso[3],thispot,
                                      phi=thiso[5],t=t)\
                                      +thiso[1]**2./2.\
                                      +thiso[2]**2./2.\
                                      +thiso[4]**2./2.
        else:
            return nu.array([evaluatePotentials(thiso[0,ii],thiso[3,ii],
                                                pot,phi=thiso[5,ii],
                                                t=t[ii])\
                                 +thiso[1,ii]**2./2.\
                                 +thiso[2,ii]**2./2.\
                                 +thiso[4,ii]**2./2. for ii in range(len(t))])

    def e(self,analytic=False,pot=None):
        """
        NAME:
           e
        PURPOSE:
           calculate the eccentricity
        INPUT:
           analytic - compute this analytically
           pot - potential to use for analytical calculation
        OUTPUT:
           eccentricity
        HISTORY:
           2010-09-15 - Written - Bovy (NYU)
        """
        if analytic:
            if not hasattr(self,'_aA'):
                self._setupaA(pot=pot)
            (rperi,rap)= self._aA.calcRapRperi()
            return (rap-rperi)/(rap+rperi)
        if not hasattr(self,'orbit'):
            raise AttributeError("Integrate the orbit first")
        if not hasattr(self,'rs'):
            self.rs= nu.sqrt(self.orbit[:,0]**2.+self.orbit[:,3]**2.)
        return (nu.amax(self.rs)-nu.amin(self.rs))/(nu.amax(self.rs)+nu.amin(self.rs))

    def rap(self,analytic=False,pot=None):
{'fraction_non_alphanumeric': 0.10415102434838572, 'fraction_numerical': 0.016089241660409738, 'mean_word_length': 1.6579247434435576, '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': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4304876', 'n_tokens_mistral': 2890, 'n_tokens_neox': 2691, 'n_words': 680}