Python Orbit.flip方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import flip [as 别名]
def get_initial_condition(to):
    """Find an initial condition near apocenter that ends up today near the end of the stream and is close to a gyrfalcon stepsize
    to= Find an apocenter that is just more recent than this time"""
    vo, ro= 220., 8.
    # Center of the stream
    o= Orbit([0.10035165,-0.81302488,0.80986668,0.58024425,0.92753945,
               0.88763126],ro=ro,vo=vo)
    #Flip for backwards integration
    of= o.flip()
    ts= numpy.linspace(0.,to/bovy_conversion.time_in_Gyr(vo,ro),10001)
    of.integrate(ts,MWPotential2014)
    # Find the closest apocenter to the time to
    rs= numpy.sqrt(of.R(ts)**2.+of.z(ts)**2.)
    drs= rs-numpy.roll(rs,1)
    nearApo= (drs < 0.)*(numpy.roll(drs,1) > 0.)
    tsNearApo= ts[nearApo]
    tsNearApo*= bovy_conversion.time_in_Gyr(vo,ro)
    tfgf= numpy.amax(tsNearApo)
    #Round to nearest 2.**-8.
    tfgf= round(tfgf/2.**-8.)*2.**-8
    tf= tfgf/bovy_conversion.time_in_Gyr(vo,ro)
    print 'Current: %s,%s,%s,%s,%s,%s' % (of.x()[0], of.y()[0], of.z(), -of.vx()[0], -of.vy()[0], -of.vz())
    print 'At %g Gyr: %s,%s,%s,%s,%s,%s' % (tf*bovy_conversion.time_in_Gyr(vo,ro),of.x(tf)[0], of.y(tf)[0], of.z(tf), -of.vx(tf,use_physical=False)[0]*bovy_conversion.velocity_in_kpcGyr(vo,ro), -of.vy(tf,use_physical=False)[0]*bovy_conversion.velocity_in_kpcGyr(vo,ro), -of.vz(tf,use_physical=False)*bovy_conversion.velocity_in_kpcGyr(vo,ro))
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import flip [as 别名]
def impulse_deltav_general_orbitintegration(v,x,b,w,x0,v0,pot,tmax,galpot,
                                            tmaxfac=10.,nsamp=1000,
                                            integrate_method='symplec4_c'):
    """
    NAME:

       impulse_deltav_general_orbitintegration

    PURPOSE:

       calculate the delta velocity to due an encounter with a general spherical potential NOT in the impulse approximation by integrating each particle in the underlying galactic potential; allows for arbitrary velocity vectors and arbitrary shaped streams.

    INPUT:

       v - velocity of the stream (nstar,3)

       x - position along the stream (nstar,3)

       b - impact parameter

       w - velocity of the subhalo (3)

       x0 - position of closest approach (3)

       v0 - velocity of stream at closest approach (3)

       pot - Potential object or list thereof (should be spherical)

       tmax - maximum integration time

       galpot - galpy Potential object or list thereof

       nsamp(1000) - number of forward integration points

       integrate_method= ('symplec4_c') orbit integrator to use (see Orbit.integrate)

    OUTPUT:

       deltav (nstar,3)

    HISTORY:

       2015-08-17 - SANDERS

    """
    if len(v.shape) == 1: v= numpy.reshape(v,(1,3))
    if len(x.shape) == 1: x= numpy.reshape(x,(1,3))
    nstar,ndim=numpy.shape(v)
    b0 = numpy.cross(w,v0)
    b0 *= b/numpy.sqrt(numpy.sum(b0**2))
    times = numpy.linspace(0.,tmax,nsamp)
    xres = numpy.zeros(shape=(len(x),nsamp*2-1,3))
    R, phi, z= bovy_coords.rect_to_cyl(x[:,0],x[:,1],x[:,2])
    vR, vp, vz= bovy_coords.rect_to_cyl_vec(v[:,0],v[:,1],v[:,2],
                                            R,phi,z,cyl=True)
    for i in range(nstar):
        o = Orbit([R[i],vR[i],vp[i],z[i],vz[i],phi[i]])
        o.integrate(times,galpot,method=integrate_method)
        xres[i,nsamp:,0]=o.x(times)[1:]
        xres[i,nsamp:,1]=o.y(times)[1:]
        xres[i,nsamp:,2]=o.z(times)[1:]
        oreverse = o.flip()
        oreverse.integrate(times,galpot,method=integrate_method)
        xres[i,:nsamp,0]=oreverse.x(times)[::-1]
        xres[i,:nsamp,1]=oreverse.y(times)[::-1]
        xres[i,:nsamp,2]=oreverse.z(times)[::-1]
    times = numpy.concatenate((-times[::-1],times[1:]))
    nsamp = len(times)
    X = b0+xres-x0-numpy.outer(times,w)
    r = numpy.sqrt(numpy.sum(X**2,axis=-1))
    acc = (numpy.reshape(evaluateRforces(r.flatten(),0.,pot),(nstar,nsamp))/r)[:,:,numpy.newaxis]*X
    return integrate.simps(acc,x=times,axis=1)

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import flip [as 别名]
    def _determine_impact_coordtransform(self,deltaAngleTrackImpact,
                                         nTrackChunksImpact,
                                         timpact,impact_angle):
        """Function that sets up the transformation between (x,v) and (O,theta)"""
        # Integrate the progenitor backward to the time of impact
        self._gap_progenitor_setup()
        # Sign of delta angle tells us whether the impact happens to the
        # leading or trailing arm, self._sigMeanSign contains this info
        if impact_angle > 0.:
            self._gap_leading= True
        else:
            self._gap_leading= False
        if (self._gap_leading and not self._leading) \
                or (not self._gap_leading and self._leading):
            raise ValueError('Modeling leading (trailing) impact for trailing (leading) arm; this is not allowed because it is nonsensical in this framework')
        self._impact_angle= numpy.fabs(impact_angle)
        self._gap_sigMeanSign= 1.
        if (self._gap_leading and self._progenitor_Omega_along_dOmega/self._sigMeanSign < 0.) \
                or (not self._gap_leading and self._progenitor_Omega_along_dOmega/self._sigMeanSign > 0.):
            self._gap_sigMeanSign= -1.
        # Determine how much orbital time is necessary for the progenitor's orbit at the time of impact to cover the part of the stream near the impact; we cover the whole leading (or trailing) part of the stream
        if nTrackChunksImpact is None:
            #default is floor(self._deltaAngleTrackImpact/0.15)+1
            self._nTrackChunksImpact= int(numpy.floor(self._deltaAngleTrackImpact/0.15))+1
        else:
            self._nTrackChunksImpact= nTrackChunksImpact
        dt= self._deltaAngleTrackImpact\
            /self._progenitor_Omega_along_dOmega\
            /self._sigMeanSign*self._gap_sigMeanSign
        self._gap_trackts= numpy.linspace(0.,2*dt,2*self._nTrackChunksImpact-1) #to be sure that we cover it
        #Instantiate an auxiliaryTrack, which is an Orbit instance at the mean frequency of the stream, and zero angle separation wrt the progenitor; prog_stream_offset is the offset between this track and the progenitor at zero angle (same as in streamdf, but just done at the time of impact rather than the current time)
        prog_stream_offset=\
            _determine_stream_track_single(self._aA,
                                           self._gap_progenitor,
                                           self._timpact, # around the t of imp
                                           self._progenitor_angle-self._timpact*self._progenitor_Omega,
                                           self._gap_sigMeanSign,
                                           self._dsigomeanProgDirection,
                                           lambda da: super(streamgapdf,self).meanOmega(da,offset_sign=self._gap_sigMeanSign,tdisrupt=self._tdisrupt-self._timpact),
                                           0.) #angle = 0
        auxiliaryTrack= Orbit(prog_stream_offset[3])
        if dt < 0.:
            self._gap_trackts= numpy.linspace(0.,-2.*dt,2.*self._nTrackChunksImpact-1)
            #Flip velocities before integrating
            auxiliaryTrack= auxiliaryTrack.flip()
        auxiliaryTrack.integrate(self._gap_trackts,self._pot)
        if dt < 0.:
            #Flip velocities again
            auxiliaryTrack._orb.orbit[:,1]= -auxiliaryTrack._orb.orbit[:,1]
            auxiliaryTrack._orb.orbit[:,2]= -auxiliaryTrack._orb.orbit[:,2]
            auxiliaryTrack._orb.orbit[:,4]= -auxiliaryTrack._orb.orbit[:,4]
        #Calculate the actions, frequencies, and angle for this auxiliary orbit
        acfs= self._aA.actionsFreqs(auxiliaryTrack(0.),maxn=3)
        auxiliary_Omega= numpy.array([acfs[3],acfs[4],acfs[5]]).reshape(3\
)
        auxiliary_Omega_along_dOmega= \
            numpy.dot(auxiliary_Omega,self._dsigomeanProgDirection)
        # compute the transformation using _determine_stream_track_single
        allAcfsTrack= numpy.empty((self._nTrackChunksImpact,9))
        alljacsTrack= numpy.empty((self._nTrackChunksImpact,6,6))
        allinvjacsTrack= numpy.empty((self._nTrackChunksImpact,6,6))
        thetasTrack= numpy.linspace(0.,self._deltaAngleTrackImpact,
                                    self._nTrackChunksImpact)
        ObsTrack= numpy.empty((self._nTrackChunksImpact,6))
        ObsTrackAA= numpy.empty((self._nTrackChunksImpact,6))
        detdOdJps= numpy.empty((self._nTrackChunksImpact))
        if self._multi is None:
            for ii in range(self._nTrackChunksImpact):
                multiOut= _determine_stream_track_single(self._aA,
                                           auxiliaryTrack,
                                           self._gap_trackts[ii]*numpy.fabs(self._progenitor_Omega_along_dOmega/auxiliary_Omega_along_dOmega), #this factor accounts for the difference in frequency between the progenitor and the auxiliary track, no timpact bc gap_tracks is relative to timpact
                                           self._progenitor_angle-self._timpact*self._progenitor_Omega,
                                           self._gap_sigMeanSign,
                                           self._dsigomeanProgDirection,
                                           lambda da: super(streamgapdf,self).meanOmega(da,offset_sign=self._gap_sigMeanSign,tdisrupt=self._tdisrupt-self._timpact),
                                           thetasTrack[ii])
                allAcfsTrack[ii,:]= multiOut[0]
                alljacsTrack[ii,:,:]= multiOut[1]
                allinvjacsTrack[ii,:,:]= multiOut[2]
                ObsTrack[ii,:]= multiOut[3]
                ObsTrackAA[ii,:]= multiOut[4]
                detdOdJps[ii]= multiOut[5]
        else:
            multiOut= multi.parallel_map(\
                (lambda x: _determine_stream_track_single(self._aA,
                                           auxiliaryTrack,
                                           self._gap_trackts[x]*numpy.fabs(self._progenitor_Omega_along_dOmega/auxiliary_Omega_along_dOmega), #this factor accounts for the difference in frequency between the progenitor and the auxiliary track, no timpact bc gap_tracks is relative to timpact
                                           self._progenitor_angle-self._timpact*self._progenitor_Omega,
                                           self._gap_sigMeanSign,
                                           self._dsigomeanProgDirection,
                                           lambda da: super(streamgapdf,self).meanOmega(da,offset_sign=self._gap_sigMeanSign,tdisrupt=self._tdisrupt-self._timpact),
                                           thetasTrack[x])),
                range(self._nTrackChunksImpact),
                numcores=numpy.amin([self._nTrackChunksImpact,
                                     multiprocessing.cpu_count(),
                                     self._multi]))
            for ii in range(self._nTrackChunksImpact):
                allAcfsTrack[ii,:]= multiOut[ii][0]
                alljacsTrack[ii,:,:]= multiOut[ii][1]
                allinvjacsTrack[ii,:,:]= multiOut[ii][2]
#.........这里部分代码省略.........

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import flip [as 别名]
TIME = 3000*units.Myr
dt = 1 * units.Myr
div = 1.
step_size = dt/div

##Creating the Milky Way Potential
ps= potential.PowerSphericalPotentialwCutoff(alpha=1.8,rc=1.9/8.,normalize=0.05)
mn= potential.MiyamotoNagaiPotential(a=3./8.,b=0.28/8.,normalize=.6)
MWPotential= [ps,mn,scf]

##Creating the orbit
o = Orbit([229.018,-0.124,23.2,-2.296,-2.257,-58.7],radec=True,ro=8.,vo=220.,solarmotion=[-11.1,24.,7.25])
o.turn_physical_off()

##Integrate backwards in time
o = o.flip()
ts= nu.arange(0,(TIME + step_size).value,step_size.value)*units.Myr

o.integrate(ts, MWPotential, method='dopr54_c')

##Integrating Forward in time
newOrbit = Orbit([o.R(TIME), -o.vR(TIME), -o.vT(TIME), o.z(TIME), -o.vz(TIME), o.phi(TIME)],ro=8.,vo=220.)
newOrbit.turn_physical_off()
newOrbit.integrate(ts, MWPotential, method='dopr54_c')


def randomVelocity(std=.001):
    if type(std).__name__ == "Quantity":
        return nu.random.normal(scale=std.value)*std.unit
    return nu.random.normal(scale=std)