Python Orbit.y方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
def evolveorbit(icon, ti, tau, pot):
    global x
    global y
    o = Orbit(vxvv=icon) # [R,vR,vT,z,vz,phi]
    tf = ti+tau
    ts = np.linspace(ti,tf,100)
    o.integrate(ts, pot, method = 'leapfrog')
    x.append(o.x(ts[0]))
    y.append(o.y(ts[0]))
    return [o.R(tf),o.vR(tf),o.vT(tf),o.z(tf),o.vz(tf),o.phi(tf)]

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
def integrate_xyzuvw(params,ts,lsr_orbit,MWPotential2014):
    """Convenience function"""
    vxvv = params.copy()
    vxvv[2]=1.0/params[2]
    o = Orbit(vxvv=vxvv, radec=True, solarmotion='schoenrich')
    o.integrate(ts,MWPotential2014,method='odeint')
    xyzuvw = np.zeros( (len(ts),6) )
    xyzuvw[:,0] = 1e3*(o.x(ts)-lsr_orbit.x(ts))
    xyzuvw[:,1] = 1e3*(o.y(ts)-lsr_orbit.y(ts))
    xyzuvw[:,2] = 1e3*(o.z(ts))
    xyzuvw[:,3] = o.U(ts) - lsr_orbit.U(ts)
    xyzuvw[:,4] = o.V(ts) - lsr_orbit.V(ts)
    xyzuvw[:,5] = o.W(ts)
    return xyzuvw

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
def integrate_xyzuvw(params,times,lsr_orbit=None,Potential=MWPotential2014):
    """Convenience function. Integrates the motion of a star backwards in time.
    
    Parameters
    ----------
    params: numpy array 
        Kinematic parameters (RAdeg,DEdeg,Plx,pmRA,pmDE,RV)
        
    ts: times for back propagation, in Myr.
        
    lsr_orbit:
        WARNING: messy...
        lsr_orbit= Orbit(vxvv=[1.,0,1,0,0.,0],vo=220,ro=8, solarmotion='schoenrich')
        lsr_orbit.integrate(ts,MWPotential2014,method='odeint')
    MWPotential2014:
        WARNING: messy...
        from galpy.potential import MWPotential2014        
    """
    #We allow MWPotential and lsr_orbit to be passed for speed, but can compute/import 
    #now.
    if not lsr_orbit:
        lsr_orbit = get_lsr_orbit(times)
        
    #Convert times to galpy units:
    ts = -(times/1e3)/bovy_conversion.time_in_Gyr(220.,8.)
        
    params = np.array(params)
    vxvv = params.copy()
    #We'd prefer to pass parallax in mas, but distance in kpc is accepted. This
    #reciporical could be done elsewhere...
    vxvv[2]=1.0/params[2]   
    o = Orbit(vxvv=vxvv, radec=True, solarmotion='schoenrich')
    
    o.integrate(ts,Potential)#,method='odeint')
    xyzuvw = np.zeros( (len(ts),6) )
    xyzuvw[:,0] = 1e3*(o.x(ts)-lsr_orbit.x(ts))
    xyzuvw[:,1] = 1e3*(o.y(ts)-lsr_orbit.y(ts))
    #The lsr_orbit is zero in the z direction by definition - the local standard of 
    #rest is at the midplane of the Galaxy
    xyzuvw[:,2] = 1e3*(o.z(ts))  
    #UVW is relative to the sun. We *could* have used vx, vy and vz. Would these
    #have been relative to the LSR?
    xyzuvw[:,3] = o.U(ts) - lsr_orbit.U(ts)
    xyzuvw[:,4] = o.V(ts) - lsr_orbit.V(ts)
    xyzuvw[:,5] = o.W(ts) - lsr_orbit.W(ts) #NB This line changed !!!
    return xyzuvw

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
def plot_stream_xz(plotfilename):
    #Read stream
    data= numpy.loadtxt(os.path.join(_STREAMSNAPDIR,'gd1_evol_hitres_00800.dat'),
                        delimiter=',')
    aadata= numpy.loadtxt(os.path.join(_STREAMSNAPAADIR,
                                       'gd1_evol_hitres_aa_00800.dat'),
                        delimiter=',')
    thetar= aadata[:,6]
    thetar= (numpy.pi+(thetar-numpy.median(thetar))) % (2.*numpy.pi)
    if 'sim' in plotfilename:
        sindx= numpy.fabs(thetar-numpy.pi) > (4.*numpy.median(numpy.fabs(thetar-numpy.median(thetar))))
    else:
        sindx= numpy.fabs(thetar-numpy.pi) > (1.5*numpy.median(numpy.fabs(thetar-numpy.median(thetar))))
    includeorbit= True
    if includeorbit:
        npts= 201
        pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9)
        pts= numpy.linspace(0.,4.,npts)
        #Calculate progenitor orbit around this point
        pox= numpy.median(data[:,1])
        poy= numpy.median(data[:,3])
        poz= numpy.median(data[:,2])
        povx= numpy.median(data[:,4])
        povy= numpy.median(data[:,6])
        povz= numpy.median(data[:,5])
        pR,pphi,pZ= bovy_coords.rect_to_cyl(pox,poy,poz)
        pvR,pvT,pvZ= bovy_coords.rect_to_cyl_vec(povx,povy,povz,pR,
                                                 pphi,pZ,cyl=True)
        ppo= Orbit([pR/8.,pvR/220.,pvT/220.,pZ/8.,pvZ/220.,pphi])
        pno= Orbit([pR/8.,-pvR/220.,-pvT/220.,pZ/8.,-pvZ/220.,pphi])
        ppo.integrate(pts,pot)
        pno.integrate(pts,pot)
        pvec= numpy.zeros((3,npts*2-1))
        pvec[0,:npts-1]= pno.x(pts)[::-1][:-1]
        pvec[1,:npts-1]= pno.z(pts)[::-1][:-1]
        pvec[2,:npts-1]= pno.y(pts)[::-1][:-1]
        pvec[0,npts-1:]= ppo.x(pts)
        pvec[1,npts-1:]= ppo.z(pts)
        pvec[2,npts-1:]= ppo.y(pts)
        pvec*= 8.
    includetrack= True
    if includetrack:
        #Setup stream model
        lp= potential.LogarithmicHaloPotential(q=0.9,normalize=1.)
        aAI= actionAngleIsochroneApprox(b=0.8,pot=lp)
        obs= numpy.array([1.56148083,0.35081535,-1.15481504,
                          0.88719443,-0.47713334,0.12019596])
        #Obs is at time 1312, need to go back 2 Gyr to time 800
        obs[1]*= -1.
        obs[2]*= -1.
        obs[4]*= -1.
        o= Orbit(obs)
        ts= numpy.linspace(0.,2.*977.7922212082034/1000./bovy_conversion.time_in_Gyr(220.,8.),1001)
        o.integrate(ts,lp)
        obs= o(ts[-1])._orb.vxvv
        obs[1]*= -1.
        obs[2]*= -1.
        obs[4]*= -1.       
        tdisrupt= 4.5-2.*977.7922212082034/1000.
        sdf= streamdf(_SIGV/220.,progenitor=Orbit(obs),pot=lp,aA=aAI,
                      leading=True,nTrackChunks=_NTRACKCHUNKS,
                      tdisrupt=tdisrupt/bovy_conversion.time_in_Gyr(220.,8.),
                      deltaAngleTrack=1.5*3./5.,multi=_NTRACKCHUNKS)
        sdft= streamdf(_SIGV/220.,progenitor=Orbit(obs),pot=lp,aA=aAI,
                       leading=False,nTrackChunks=_NTRACKCHUNKS,
                       tdisrupt=tdisrupt/bovy_conversion.time_in_Gyr(220.,8.),
                       deltaAngleTrack=1.5*3./5.,multi=_NTRACKCHUNKS)
    if 'sim' in plotfilename:
        #Replace data with simulated data
        forwardXY= sdf.sample(int(round(numpy.sum(sindx)/2.)),
                              xy=True)
        backwardXY= sdft.sample(int(round(numpy.sum(sindx)/2.)),
                                xy=True)
        data= numpy.empty((forwardXY.shape[1]+backwardXY.shape[1],7))
        data[:forwardXY.shape[1],1]= forwardXY[0]*8.
        data[:forwardXY.shape[1],2]= forwardXY[2]*8.
        data[:forwardXY.shape[1],3]= forwardXY[1]*8.
        data[:forwardXY.shape[1],4]= forwardXY[3]*220.
        data[:forwardXY.shape[1],5]= forwardXY[5]*220.
        data[:forwardXY.shape[1],6]= forwardXY[4]*220.
        data[forwardXY.shape[1]:,1]= backwardXY[0]*8.
        data[forwardXY.shape[1]:,2]= backwardXY[2]*8.
        data[forwardXY.shape[1]:,3]= backwardXY[1]*8.
        data[forwardXY.shape[1]:,4]= backwardXY[3]*220.
        data[forwardXY.shape[1]:,5]= backwardXY[5]*220.
        data[forwardXY.shape[1]:,6]= backwardXY[4]*220.
        sindx= numpy.ones(data.shape[0],dtype='bool')      
    #Plot
    bovy_plot.bovy_print()
    bovy_plot.bovy_plot(data[sindx,1],data[sindx,2],'k,',ms=2.,
                        xlabel=r'$X\,(\mathrm{kpc})$',
                        ylabel=r'$Z\,(\mathrm{kpc})$',
                        xrange=[-12.5,-3.],
                        yrange=[-12.5,-7.])
    if numpy.sum(True-sindx) > 0:
        #Also plot progenitor
        pindx= copy.copy(True-sindx)
        pindx[0:9900]= False #subsample
        bovy_plot.bovy_plot(data[pindx,1],data[pindx,2],
                            'k,',overplot=True)
#.........这里部分代码省略.........

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
def make_sim_movie(proj='xz',comov=False,skippng=False,
                   includeorbit=True):
    #Directories
    savedirpng= './movies/oph/pngs/'
    timestr= '_1Gyr'
#    timestr= ''
#    massstr= '_lowmass'
    massstr= ''
    basefilename= 'oph%s_evol%s_' % (massstr,timestr)
    moviefilename= 'oph%s_evol%s' % (massstr,timestr)
    #Read data
    #datafile= 'oph_evol_hitres.dat'
    datafile= 'oph%s_evol%s.dat' % (massstr,timestr)
    print "Reading data ..."
    data= numpy.loadtxt(datafile,comments='#')
    print "Done reading data"
    if proj.lower() == 'xz':
        includeorbit= False #just to be sure
        basefilename+= 'xz_'
        moviefilename+= '_xz'
        x= data[:,1]
        y= data[:,3]
        if comov:
            basefilename+= 'comov_'
            moviefilename+= '_comov'
            xrange=[-12.,12.]
            yrange=[-10.,10.]           
        else:
            xrange=[-18.,18.]
            yrange=[-18.,18.]           
        xlabel=r'$X\,(\mathrm{kpc})$'
        ylabel=r'$Z\,(\mathrm{kpc})$'
    elif proj.lower() == 'yz':
        includeorbit= False #just to be sure
        basefilename+= 'yz_'
        moviefilename+= '_yz'
        x= data[:,2]
        y= data[:,3]
        if comov:
            basefilename+= 'comov_'
            moviefilename+= '_comov'
            xrange=[-15.,15.]
            yrange=[-15.,15.]           
        else:
            xrange=[-18.,18.]
            yrange=[-18.,18.]           
        xlabel=r'$Y\,(\mathrm{kpc})$'
        ylabel=r'$Z\,(\mathrm{kpc})$'
    elif proj.lower() == 'orbplane':
        basefilename+= 'orbplane_'
        moviefilename+= '_orbplane'
        x= numpy.zeros_like(data[:,1])
        y= numpy.zeros_like(data[:,2])
        nx= 20000
        nt= len(x)/nx
        diff= numpy.empty(nt)
        if includeorbit:
            npts= 201
            pot= potential.MWPotential2014
            pts= numpy.linspace(0.,4.,npts)
            px= numpy.zeros(nt*(2*npts-1))
            py= numpy.zeros(nt*(2*npts-1))
        ii= 0
        #Calculate median angular momentum at t=0, use this to always go to the orbital plane
        Lx= numpy.median(data[ii*nx:(ii+1)*nx,2]*data[ii*nx:(ii+1)*nx,6]\
                             -data[ii*nx:(ii+1)*nx,3]*data[ii*nx:(ii+1)*nx,5])
        Ly= numpy.median(data[ii*nx:(ii+1)*nx,3]*data[ii*nx:(ii+1)*nx,4]\
                             -data[ii*nx:(ii+1)*nx,1]*data[ii*nx:(ii+1)*nx,6])
        Lz= numpy.median(data[ii*nx:(ii+1)*nx,1]*data[ii*nx:(ii+1)*nx,5]\
                             -data[ii*nx:(ii+1)*nx,2]*data[ii*nx:(ii+1)*nx,4])
        L= numpy.sqrt(Lx**2.+Ly**2.+Lz**2.)
        Lx/= L
        Ly/= L
        Lz/= L
        Txz= numpy.zeros((3,3))
        Tz= numpy.zeros((3,3))
        Txz[0,0]= Lx/numpy.sqrt(Lx**2.+Ly**2.)
        Txz[1,1]= Lx/numpy.sqrt(Lx**2.+Ly**2.)
        Txz[1,0]= Ly/numpy.sqrt(Lx**2.+Ly**2.)
        Txz[0,1]= -Ly/numpy.sqrt(Lx**2.+Ly**2.)
        Txz[2,2]= 1.
        Tz[0,0]= Lz
        Tz[1,1]= 1.
        Tz[2,2]= Lz
        Tz[2,0]= -numpy.sqrt(Lx**2.+Ly**2.)
        Tz[0,2]= numpy.sqrt(Lx**2.+Ly**2.)
        TL= numpy.dot(Tz,Txz)
        for ii in range(nt):
            if includeorbit:
                #Calculate progenitor orbit around this point
                pox= numpy.median(data[ii*nx:(ii+1)*nx,1])
                poy= numpy.median(data[ii*nx:(ii+1)*nx,3])
                poz= numpy.median(data[ii*nx:(ii+1)*nx,2])
                povx= numpy.median(data[ii*nx:(ii+1)*nx,4])
                povy= numpy.median(data[ii*nx:(ii+1)*nx,6])
                povz= numpy.median(data[ii*nx:(ii+1)*nx,5])
                pR,pphi,pZ= bovy_coords.rect_to_cyl(pox,poy,poz)
                pvR,pvT,pvZ= bovy_coords.rect_to_cyl_vec(povx,povy,povz,pR,
                                                         pphi,pZ,cyl=True)
                ppo= Orbit([pR/8.,pvR/220.,pvT/220.,pZ/8.,pvZ/220.,pphi])
#.........这里部分代码省略.........

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [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 y [as 别名]
def trace_forward(xyzuvw, time_in_past,
                  Potential=MWPotential2014, solarmotion=None):
    """Trace forward one star in xyzuvw coords

    Parameters
    ----------
    xyzuvw: numpy float array
        xyzuvw relative to the local standard of rest at some point in the past.

    time_in_past: float
        Time in the past that we want to trace forward to the present.
    """
    if solarmotion == None:
        xyzuvw_sun = np.zeros(6)
    elif solarmotion == 'schoenrich':
        xyzuvw_sun = [0, 0, 25, 11.1, 12.24, 7.25]
    else:
        raise UserWarning
    if time_in_past == 0:
        time_in_past = 1e-5
        # print("Time in past must be nonzero: {}".format(time_in_past))
        # raise UserWarning
    times = np.linspace(0, time_in_past, 2)

    # Start off with an LSR orbit...
    lsr_orbit = get_lsr_orbit(times)

    # Convert times to galpy units:
    ts = -(times / 1e3) / bovy_conversion.time_in_Gyr(220., 8.)

    # Add on the lsr_orbit to get the times in the past.
    xyzuvw_gal = np.zeros(6)
    xyzuvw_gal[0] = xyzuvw[0] / 1e3 + lsr_orbit.x(ts[-1]) - lsr_orbit.x(ts[0])
    xyzuvw_gal[1] = xyzuvw[1] / 1e3 + lsr_orbit.y(ts[-1]) - lsr_orbit.y(ts[0])
    # Relative to a zero-height orbit, excluding the solar motion.
    xyzuvw_gal[2] = xyzuvw[2] / 1e3
    xyzuvw_gal[3] = xyzuvw[3] + lsr_orbit.U(ts[-1]) - lsr_orbit.U(ts[0])
    xyzuvw_gal[4] = xyzuvw[4] + lsr_orbit.V(ts[-1]) - lsr_orbit.V(ts[0])
    xyzuvw_gal[5] = xyzuvw[5] + lsr_orbit.W(ts[-1]) - lsr_orbit.W(ts[0])

    # Now convert to units that galpy understands by default.
    # FIXME: !!! Reverse [0] sign here.
    l = np.degrees(np.arctan2(xyzuvw_gal[1], -xyzuvw_gal[0]))
    b = np.degrees(
        np.arctan2(xyzuvw_gal[2], np.sqrt(np.sum(xyzuvw_gal[:2] ** 2))))
    dist = np.sqrt(np.sum(xyzuvw_gal[:3] ** 2))
    # As we are already in LSR co-ordinates, put in zero for solar motion and
    #  zo.
    o = Orbit([l, b, dist, xyzuvw_gal[3], xyzuvw_gal[4], xyzuvw_gal[5]],
              solarmotion=[0, 0, 0], zo=0, lb=True, uvw=True, vo=220, ro=8)

    # Integrate this orbit forwards in time.
    o.integrate(-ts, Potential)  # ,method='odeint')

    # Add in the solar z and the solar motion.
    xyzuvw_now = np.zeros(6)
    xyzuvw_now[0] = 1e3 * (o.x(-ts[-1]) - lsr_orbit.x(0))
    xyzuvw_now[1] = 1e3 * o.y(-ts[-1])
    xyzuvw_now[2] = 1e3 * o.z(-ts[-1]) - xyzuvw_sun[2]
    xyzuvw_now[3] = o.U(-ts[-1]) - xyzuvw_sun[3]
    xyzuvw_now[4] = o.V(-ts[-1]) - xyzuvw_sun[4]
    xyzuvw_now[5] = o.W(-ts[-1]) - xyzuvw_sun[5]

    # pdb.set_trace()

    return xyzuvw_now

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
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)

time1 = nu.arange(0, TIME.value, dt.value)*units.Myr
orbits_pos = nu.empty((len(time1) + 1,9,len(ts)), dtype=units.quantity.Quantity)
orbits_pos[0, :, :] = ts, newOrbit.x(ts), newOrbit.y(ts), newOrbit.z(ts), newOrbit.vx(ts), newOrbit.vy(ts), newOrbit.vz(ts), newOrbit.ra(ts), newOrbit.dec(ts)
orbits_pos[:,:,:] = orbits_pos[0,:,:]
i = 0
std = 0.004
stdR = std
stdT = std
stdz = std
for t in time1:
    print t
    dvR = randomVelocity(stdR)
    dvT = randomVelocity(stdT)
    dvz = randomVelocity(stdz)
    #dvR, dvT, dvz = 0,0,0
    tempOrbit = Orbit([newOrbit.R(t), newOrbit.vR(t) + dvR, newOrbit.vT(t) + dvT, newOrbit.z(t), newOrbit.vz(t) + dvz, newOrbit.phi(t)],ro=8.,vo=220.)
    tempOrbit.turn_physical_off()
    time = nu.arange(0,(TIME + step_size - t).value,step_size.value)*units.Myr

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
        #ra,dec,Dist,UVel,VVel,WVel,RVel,pmRA,pmDE

        orbit = Orbit(vxvv=[ra, dec, distance, pm_ra, pm_dec, velocity], radec=True, uvw=False, lb=False, vo=v0_mc, ro=r0_mc, zo=z0_mc, solarmotion='schoenrich')

        distancesJ.append(distanceJ)
        distancesK.append(distanceK)

        velocity_field = [
                ra,
                dec,
                distance,
                pm_ra,
                pm_dec,
                velocity,
                orbit.x() * r0,
                orbit.y() * r0,
                orbit.z() * r0,
                orbit.U()[0],
                orbit.V()[0],
                orbit.W()[0],
                orbit.vx()[0] * v0_mc,
                orbit.vy()[0] * v0_mc,
                orbit.vz() * v0_mc,
                orbit.vR() * v0_mc,
                orbit.vT() * v0_mc,
                orbit.vphi()[0] * v0_mc,
                orbit.L()[0][2] * v0 * r0, # Lz, 18
                orbit.E(pot=MWPotential) * pow(v0, 2)/2, # Energy
                pow(np.sum(pow(orbit.L()[0][:2], 2)), 0.5) * v0 * r0 # Lperp
            ]
            

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
def plot_stream_xz(plotfilename):
    #Read stream
    data= numpy.loadtxt(os.path.join(_STREAMSNAPDIR,'gd1_evol_hitres_01312.dat'),
                        delimiter=',')
    includeorbit= True
    if includeorbit:
        npts= 201
        pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9)
        pts= numpy.linspace(0.,4.,npts)
        #Calculate progenitor orbit around this point
        pox= numpy.median(data[:,1])
        poy= numpy.median(data[:,3])
        poz= numpy.median(data[:,2])
        povx= numpy.median(data[:,4])
        povy= numpy.median(data[:,6])
        povz= numpy.median(data[:,5])
        pR,pphi,pZ= bovy_coords.rect_to_cyl(pox,poy,poz)
        pvR,pvT,pvZ= bovy_coords.rect_to_cyl_vec(povx,povy,povz,pR,
                                                 pphi,pZ,cyl=True)
        ppo= Orbit([pR/8.,pvR/220.,pvT/220.,pZ/8.,pvZ/220.,pphi])
        pno= Orbit([pR/8.,-pvR/220.,-pvT/220.,pZ/8.,-pvZ/220.,pphi])
        ppo.integrate(pts,pot)
        pno.integrate(pts,pot)
        pvec= numpy.zeros((3,npts*2-1))
        pvec[0,:npts-1]= pno.x(pts)[::-1][:-1]
        pvec[1,:npts-1]= pno.z(pts)[::-1][:-1]
        pvec[2,:npts-1]= pno.y(pts)[::-1][:-1]
        pvec[0,npts-1:]= ppo.x(pts)
        pvec[1,npts-1:]= ppo.z(pts)
        pvec[2,npts-1:]= ppo.y(pts)
        pvec*= 8.
    includetrack= True
    if includetrack:
        #Setup stream model
        lp= potential.LogarithmicHaloPotential(q=0.9,normalize=1.)
        aAI= actionAngleIsochroneApprox(b=0.8,pot=lp)
        obs= numpy.array([1.56148083,0.35081535,-1.15481504,
                          0.88719443,-0.47713334,0.12019596])
        sdf= streamdf(_SIGV/220.,progenitor=Orbit(obs),pot=lp,aA=aAI,
                      leading=True,nTrackChunks=_NTRACKCHUNKS,
                      tdisrupt=4.5/bovy_conversion.time_in_Gyr(220.,8.))
        sdft= streamdf(_SIGV/220.,progenitor=Orbit(obs),pot=lp,aA=aAI,
                       leading=False,nTrackChunks=_NTRACKCHUNKS,
                       tdisrupt=4.5/bovy_conversion.time_in_Gyr(220.,8.))
    #Plot
    bovy_plot.bovy_print()
    bovy_plot.bovy_plot(data[:,1],data[:,2],'k,',
                        xlabel=r'$X\,(\mathrm{kpc})$',
                        ylabel=r'$Z\,(\mathrm{kpc})$',
                        xrange=[0.,16.],
                        yrange=[-0.5,11.])
    if includeorbit:
        bovy_plot.bovy_plot(pox,poz,'o',color='0.5',mec='none',overplot=True,ms=8)
        bovy_plot.bovy_plot(pvec[0,:],pvec[1,:],'k--',overplot=True,lw=1.)
    if includetrack:
        d1= 'x'
        d2= 'z'
        sdf.plotTrack(d1=d1,d2=d2,interp=True,color='k',spread=0,
                      overplot=True,lw=1.,scaleToPhysical=True)
        sdft.plotTrack(d1=d1,d2=d2,interp=True,color='k',spread=0,
                       overplot=True,lw=1.,scaleToPhysical=True)
        #Also create inset
        pyplot.plot([12.,12.],[0.5,7.5],'k-')
        pyplot.plot([14.5,14.5],[0.5,7.5],'k-')
        pyplot.plot([12.,14.5],[0.5,0.5],'k-')
        pyplot.plot([12.,14.5],[7.5,7.5],'k-')
        pyplot.plot([12.,8.8],[7.5,7.69],'k:')
        pyplot.plot([12.,8.8],[0.5,-0.21],'k:')
        insetAxes= pyplot.axes([0.15,0.12,0.4,0.55])
        pyplot.sca(insetAxes)
        bovy_plot.bovy_plot(data[:,1],data[:,2],'k,',
                            overplot=True)
        bovy_plot.bovy_plot(pvec[0,:],pvec[1,:],'k--',overplot=True,lw=1.)
        sdf.plotTrack(d1=d1,d2=d2,interp=True,color='k',spread=0,
                      overplot=True,lw=1.,scaleToPhysical=True)
        nullfmt   = NullFormatter()         # no labels
        insetAxes.xaxis.set_major_formatter(nullfmt)
        insetAxes.yaxis.set_major_formatter(nullfmt)
        insetAxes.set_xlim(12.,14.5)
        insetAxes.set_ylim(.5,7.5)
        pyplot.tick_params(\
            axis='both',          # changes apply to the x-axis
            which='both',      # both major and minor ticks are affected
            bottom='off',      # ticks along the bottom edge are off
            top='off',         # ticks along the top edge are off
            left='off',      # ticks along the bottom edge are off
            right='off')         # ticks along the top edge are off
    bovy_plot.bovy_end_print(plotfilename)

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
def plot_stream_lb(plotfilename):
    #Read stream
    data= numpy.loadtxt(os.path.join(_STREAMSNAPDIR,'gd1_evol_hitres_01312.dat'),
                        delimiter=',')
    aadata= numpy.loadtxt(os.path.join(_STREAMSNAPAADIR,
                                       'gd1_evol_hitres_aa_01312.dat'),
                          delimiter=',')
    thetar= aadata[:,6]
    thetar= (numpy.pi+(thetar-numpy.median(thetar))) % (2.*numpy.pi)
    sindx= numpy.fabs(thetar-numpy.pi) > (1.5*numpy.median(numpy.fabs(thetar-numpy.median(thetar)))) #stars in the stream
    #Transform to (l,b)
    XYZ= bovy_coords.galcenrect_to_XYZ(data[:,1],data[:,3],data[:,2],Xsun=8.)
    lbd= bovy_coords.XYZ_to_lbd(XYZ[0],XYZ[1],XYZ[2],degree=True)
    vXYZ= bovy_coords.galcenrect_to_vxvyvz(data[:,4],data[:,6],data[:,5],
                                           vsun=[0.,30.24*8.,0.])
    vlbd= bovy_coords.vxvyvz_to_vrpmllpmbb(vXYZ[0],vXYZ[1],vXYZ[2],
                                           lbd[:,0],lbd[:,1],lbd[:,2],
                                           degree=True)
    includeorbit= True
    if includeorbit:
        npts= 201
        pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9)
        pts= numpy.linspace(0.,4.,npts)
        #Calculate progenitor orbit around this point
        pox= numpy.median(data[:,1])
        poy= numpy.median(data[:,3])
        poz= numpy.median(data[:,2])
        povx= numpy.median(data[:,4])
        povy= numpy.median(data[:,6])
        povz= numpy.median(data[:,5])
        pR,pphi,pZ= bovy_coords.rect_to_cyl(pox,poy,poz)
        pvR,pvT,pvZ= bovy_coords.rect_to_cyl_vec(povx,povy,povz,pR,
                                                 pphi,pZ,cyl=True)
        ppo= Orbit([pR/8.,pvR/220.,pvT/220.,pZ/8.,pvZ/220.,pphi])
        pno= Orbit([pR/8.,-pvR/220.,-pvT/220.,pZ/8.,-pvZ/220.,pphi])
        ppo.integrate(pts,pot)
        pno.integrate(pts,pot)
        pvec= numpy.zeros((6,npts*2-1))
        pvec[0,:npts-1]= pno.x(pts)[::-1][:-1]
        pvec[1,:npts-1]= pno.z(pts)[::-1][:-1]
        pvec[2,:npts-1]= pno.y(pts)[::-1][:-1]
        pvec[0,npts-1:]= ppo.x(pts)
        pvec[1,npts-1:]= ppo.z(pts)
        pvec[2,npts-1:]= ppo.y(pts)
        pvec[3,:npts-1]= -pno.vx(pts)[::-1][:-1]
        pvec[4,:npts-1]= -pno.vz(pts)[::-1][:-1]
        pvec[5,:npts-1]= -pno.vy(pts)[::-1][:-1]
        pvec[3,npts-1:]= ppo.vx(pts)
        pvec[4,npts-1:]= ppo.vz(pts)
        pvec[5,npts-1:]= ppo.vy(pts)
        pvec[:3,:]*= 8.
        pvec[3:,:]*= 220.
        pXYZ= bovy_coords.galcenrect_to_XYZ(pvec[0,:],pvec[2,:],pvec[1,:],
                                            Xsun=8.)
        plbd= bovy_coords.XYZ_to_lbd(pXYZ[0],pXYZ[1],pXYZ[2],degree=True)
        pvXYZ= bovy_coords.galcenrect_to_vxvyvz(pvec[3,:],pvec[5,:],pvec[4,:],
                                                vsun=[0.,30.24*8.,0.])
        pvlbd= bovy_coords.vxvyvz_to_vrpmllpmbb(pvXYZ[0],pvXYZ[1],pvXYZ[2],
                                                plbd[:,0],plbd[:,1],plbd[:,2],
                                                degree=True)
    includetrack= True
    if includetrack:
        #Setup stream model
        lp= potential.LogarithmicHaloPotential(q=0.9,normalize=1.)
        aAI= actionAngleIsochroneApprox(b=0.8,pot=lp)
        obs= numpy.array([1.56148083,0.35081535,-1.15481504,
                          0.88719443,-0.47713334,0.12019596])
        sdf= streamdf(_SIGV/220.,progenitor=Orbit(obs),pot=lp,aA=aAI,
                      leading=True,nTrackChunks=_NTRACKCHUNKS,
                      vsun=[0.,30.24*8.,0.],
                      tdisrupt=4.5/bovy_conversion.time_in_Gyr(220.,8.),
                      multi=_NTRACKCHUNKS)
        sdft= streamdf(_SIGV/220.,progenitor=Orbit(obs),pot=lp,aA=aAI,
                       leading=False,nTrackChunks=_NTRACKCHUNKS,
                       vsun=[0.,30.24*8.,0.],
                       tdisrupt=4.5/bovy_conversion.time_in_Gyr(220.,8.),
                       multi=_NTRACKCHUNKS)
    #Plot
    bovy_plot.bovy_print(fig_width=8.25,fig_height=3.5)
    if 'ld' in plotfilename:
        lbindx= 2
        ylabel=r'$\mathrm{Distance}\,(\mathrm{kpc})$'
        yrange=[0.,30.]
    elif 'lvlos' in plotfilename:
        lbindx= 0
        ylabel=r'$V_\mathrm{los}\,(\mathrm{km\,s}^{-1})$'
        yrange=[-500.,500.]
    elif 'lpmll' in plotfilename:
        lbindx= 1
        ylabel=r'$\mu_{l}\cos b\,(\mathrm{mas\,yr}^{-1})$'
        yrange=[-2.,13.5]
    elif 'lpmbb' in plotfilename:
        lbindx= 2
        ylabel=r'$\mu_{b}\,(\mathrm{mas\,yr}^{-1})$'
        yrange=[-8.,7.]
    else:
        lbindx= 1 
        yrange=[-10.,60.]
        ylabel=r'$\mathrm{Galactic\ latitude}\,(\mathrm{deg})$'
    if 'vlos' in plotfilename or 'pm' in plotfilename:
#.........这里部分代码省略.........

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import y [as 别名]
def plot_stream_xz(plotfilename):
    #Read stream
    data= numpy.loadtxt(os.path.join(_STREAMSNAPDIR,'gd1-hisigv_evol_00041.dat'),
                        delimiter=',')
    includeorbit= True
    if includeorbit:
        npts= 201
        pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9)
        pts= numpy.linspace(0.,17.,npts)
        #Calculate progenitor orbit around this point
        pox= numpy.median(data[:,1])
        poy= numpy.median(data[:,3])
        poz= numpy.median(data[:,2])
        povx= numpy.median(data[:,4])
        povy= numpy.median(data[:,6])
        povz= numpy.median(data[:,5])
        pR,pphi,pZ= bovy_coords.rect_to_cyl(pox,poy,poz)
        pvR,pvT,pvZ= bovy_coords.rect_to_cyl_vec(povx,povy,povz,pR,
                                                 pphi,pZ,cyl=True)
        ppo= Orbit([pR/8.,pvR/220.,pvT/220.,pZ/8.,pvZ/220.,pphi])
        pno= Orbit([pR/8.,-pvR/220.,-pvT/220.,pZ/8.,-pvZ/220.,pphi])
        ppo.integrate(pts,pot)
        pno.integrate(pts,pot)
        pvec= numpy.zeros((3,npts*2-1))
        pvec[0,:npts-1]= pno.x(pts)[::-1][:-1]
        pvec[1,:npts-1]= pno.z(pts)[::-1][:-1]
        pvec[2,:npts-1]= pno.y(pts)[::-1][:-1]
        pvec[0,npts-1:]= ppo.x(pts)
        pvec[1,npts-1:]= ppo.z(pts)
        pvec[2,npts-1:]= ppo.y(pts)
        pvec*= 8.
    includetrack= True
    if includetrack:
        #Setup stream model
        lp= potential.LogarithmicHaloPotential(q=0.9,normalize=1.)
        aAI= actionAngleIsochroneApprox(b=0.8,pot=lp)
        obs= numpy.array([1.56148083,0.35081535,-1.15481504,
                          0.88719443,-0.47713334,0.12019596])
        sdf= streamdf(_SIGV/220.,progenitor=Orbit(obs),pot=lp,aA=aAI,
                      leading=True,nTrackChunks=_NTRACKCHUNKS,
                      tdisrupt=4.5/bovy_conversion.time_in_Gyr(220.,8.),
                      deltaAngleTrack=13.5,multi=_NTRACKCHUNKS)
        sdft= streamdf(_SIGV/220.,progenitor=Orbit(obs),pot=lp,aA=aAI,
                       leading=False,nTrackChunks=_NTRACKCHUNKS,
                       tdisrupt=4.5/bovy_conversion.time_in_Gyr(220.,8.),
                       deltaAngleTrack=13.5,multi=_NTRACKCHUNKS)
    #Plot
    bovy_plot.bovy_print()
    bovy_plot.bovy_plot(data[:,1],data[:,2],'k,',
                        xlabel=r'$X\,(\mathrm{kpc})$',
                        ylabel=r'$Z\,(\mathrm{kpc})$',
                        xrange=[-30.,30.],
                        yrange=[-20.,20])
    if includeorbit:
        bovy_plot.bovy_plot(pox,poz,'o',color='0.5',mec='none',overplot=True,ms=8)
        bovy_plot.bovy_plot(pvec[0,:],pvec[1,:],'k--',overplot=True,lw=1.)
    if includetrack:
        d1= 'x'
        d2= 'z'
        sdf.plotTrack(d1=d1,d2=d2,interp=True,color='k',spread=0,
                      overplot=True,lw=1.,scaleToPhysical=True)
        sdft.plotTrack(d1=d1,d2=d2,interp=True,color='k',spread=0,
                       overplot=True,lw=1.,scaleToPhysical=True)
    bovy_plot.bovy_text(r'$M^p = 2\times 10^7\,M_\odot$'+'\n'+
                        r'$\sigma_v^p = 14\,\mathrm{km\,s}^{-1}$',
                        top_left=True,size=16.)
    bovy_plot.bovy_end_print(plotfilename)