Python Orbit.vy方法代码示例

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

       impulse_deltav_general_fullplummerintegration

    PURPOSE:

       calculate the delta velocity to due an encounter with a moving Plummer sphere and galactic potential relative to just in galactic potential

    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)

       galpot - Galaxy Potential object

       GM - mass of Plummer

       rs - scale of Plummer

       tmaxfac(10) - multiple of rs/|w-v0| to use for time integration interval

       N(1000) - number of forward integration points

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

    OUTPUT:

       deltav (nstar,3)

    HISTORY:

       2015-08-18 - 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))
    X = x0-b0

    # Setup Plummer orbit
    R, phi, z= bovy_coords.rect_to_cyl(X[0],X[1],X[2])
    vR, vp, vz= bovy_coords.rect_to_cyl_vec(w[0],w[1],w[2],R,phi,z,cyl=True)
    tmax = tmaxfac*rs/numpy.sqrt(numpy.sum((w-v0)**2))
    times = numpy.linspace(0.,tmax,N)
    dtimes = numpy.linspace(-tmax,tmax,2*N)
    o = Orbit(vxvv=[R,-vR,-vp,z,-vz,phi])
    o.integrate(times,galpot,method=integrate_method)
    oplum = Orbit(vxvv=[o.R(times[-1]),-o.vR(times[-1]),-o.vT(times[-1]),o.z(times[-1]),-o.vz(times[-1]),o.phi(times[-1])])
    oplum.integrate(dtimes,galpot,method=integrate_method)
    plumpot = MovingObjectPotential(orbit=oplum, GM=GM, softening_model='plummer', softening_length=rs)

    # Now integrate each particle backwards in galaxy potential, forwards in combined potential and backwards again in galaxy and take diff

    deltav = numpy.zeros((nstar,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):
      ostar= Orbit(vxvv=[R[i],-vR[i],-vp[i],z[i],-vz[i],phi[i]])
      ostar.integrate(times,galpot,method=integrate_method)
      oboth = Orbit(vxvv=[ostar.R(times[-1]),-ostar.vR(times[-1]),-ostar.vT(times[-1]),ostar.z(times[-1]),-ostar.vz(times[-1]),ostar.phi(times[-1])])
      oboth.integrate(dtimes,[galpot,plumpot],method=integrate_method)
      ogalpot = Orbit(vxvv=[oboth.R(times[-1]),-oboth.vR(times[-1]),-oboth.vT(times[-1]),oboth.z(times[-1]),-oboth.vz(times[-1]),oboth.phi(times[-1])])
      ogalpot.integrate(times,galpot,method=integrate_method)
      deltav[i][0] = -ogalpot.vx(times[-1]) - v[i][0]
      deltav[i][1] = -ogalpot.vy(times[-1]) - v[i][1]
      deltav[i][2] = -ogalpot.vz(times[-1]) - v[i][2]
    return deltav

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vy [as 别名]
        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
            ]
            
        velocities[i, :] = velocity_field


    monte_carlo_data[j, :, :] = velocities

    x_pos, y_pos = velocities[:, 6], velocities[:, 7]

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vy [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 vy [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:
#.........这里部分代码省略.........