Python Orbit.vz方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vz [as 别名]
def test_orbitIntegrationC():

    #_____initialize some KKSPot_____
    Delta = 1.0
    pot = KuzminKutuzovStaeckelPotential(ac=20.,Delta=Delta,normalize=True)

    #_____initialize an orbit (twice)_____
    vxvv = [1.,0.1,1.1,0.,0.1]
    o_P= Orbit(vxvv=vxvv)
    o_C= Orbit(vxvv=vxvv)

    #_____integrate the orbit with python and C_____
    ts= numpy.linspace(0,100,101)
    o_P.integrate(ts,pot,method='leapfrog')  #python
    o_C.integrate(ts,pot,method='leapfrog_c')#C

    for ii in range(5):
        exp3= -1.7
        if   ii == 0: Python, CC, string, exp1, exp2 = o_P.R(ts) , o_C.R(ts) , 'R' , -5., -10.
        elif ii == 1: Python, CC, string, exp1, exp2 = o_P.z(ts) , o_C.z(ts) , 'z' , -3.25, -4.
        elif ii == 2: Python, CC, string, exp1, exp2 = o_P.vR(ts), o_C.vR(ts), 'vR', -3., -10.
        elif ii == 3: Python, CC, string, exp1, exp2, exp3 = o_P.vz(ts), o_C.vz(ts), 'vz', -3., -4., -1.3
        elif ii == 4: Python, CC, string, exp1, exp2 = o_P.vT(ts), o_C.vT(ts), 'vT', -5., -10.

        rel_diff = numpy.fabs((Python-CC)/CC) < 10.**exp1
        abs_diff = (numpy.fabs(Python-CC) < 10.**exp2) * (numpy.fabs(Python) < 10.**exp3)
        assert numpy.all(rel_diff+abs_diff), \
            'Orbit integration for '+string+' coordinate different in ' + \
            'C and Python implementation.'

    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vz [as 别名]
def plot_aagrid(plotfilename1,plotfilename2):
    #Setup orbit
    E, Lz= -1.25, 0.6
    o= Orbit([0.8,0.3,Lz/0.8,0.,numpy.sqrt(2.*(E-evalPot(0.8,0.,MWPotential2014)-(Lz/0.8)**2./2.-0.3**2./2.)),0.])
    delta= 0.434
    #Integrate the orbit, setup Staeckel already to calculate the period
    aAS= actionAngleStaeckel(pot=MWPotential2014,delta=delta,c=True)
    orbt= 2.*numpy.pi/aAS.actionsFreqs(o)[4]  
    norb= 5.
    nt= 501
    ts= numpy.linspace(0.,norb*orbt,nt)
    o.integrate(ts,MWPotential2014,method='symplec4_c')
    #First do adiabatic
    aAA= actionAngleAdiabatic(pot=MWPotential2014,gamma=1.,c=True)
    aAAG= actionAngleAdiabaticGrid(pot=MWPotential2014,gamma=1.,c=True,
                                   nR=31,nEz=31,nEr=51,nLz=51)
    jfa= aAA(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
    jfag= aAAG(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
    #First do adiabatic
    #aAS already setup
    aASG= actionAngleStaeckelGrid(pot=MWPotential2014,delta=delta,c=True,
                                  nE=51,npsi=51,nLz=51)
    jfs= aAS(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
    jfsg= aASG(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
    bovy_plot.bovy_print()
    line1= bovy_plot.bovy_plot(jfa[0],jfa[2],'r.',
                               xrange=[0.045,0.055],
                               yrange=[0.0075,0.011],
                               xlabel=r'$J_R$',ylabel=r'$J_z$',zorder=2)
    line2= bovy_plot.bovy_plot(jfag[0],jfag[2],'rx',overplot=True,zorder=1)
    bovy_plot.bovy_plot(jfs[0],jfs[2],'k,',overplot=True)
    pyplot.legend((line1[0],line2[0]),
                  (r'$\mathrm{\texttt{actionAngleAdiabatic}}$',
                   r'$\mathrm{\texttt{actionAngleAdiabaticGrid}}$',),
                  loc='upper right',#bbox_to_anchor=(.91,.375),
                  numpoints=1,
                  prop={'size':14},
                  frameon=False)
    bovy_plot.bovy_end_print(plotfilename1)
    #Zoom of Staeckel
    line1= bovy_plot.bovy_plot(jfs[0],jfs[2],'k.',
                               xrange=[0.05025,0.05145],
                               yrange=[0.0086,0.00933],
                               xlabel=r'$J_R$',ylabel=r'$J_z$')
    line2= bovy_plot.bovy_plot(jfsg[0],jfsg[2],'kx',overplot=True)
    pyplot.legend((line1[0],line2[0]),
                  (r'$\mathrm{\texttt{actionAngleStaeckel}}$',
                   r'$\mathrm{\texttt{actionAngleStaeckelGrid}}$',),
                  loc='upper right',#bbox_to_anchor=(.91,.375),
                  numpoints=1,
                  prop={'size':14},
                  frameon=False)
    bovy_plot.bovy_end_print(plotfilename2)
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vz [as 别名]
def test_actionConservation():

    #_____initialize some KKSPot_____
    Delta = 1.0
    pot = KuzminKutuzovStaeckelPotential(ac=20.,Delta=Delta,normalize=True)

    #_____initialize an orbit (twice)_____
    vxvv = [1.,0.1,1.1,0.01,0.1]
    o= Orbit(vxvv=vxvv)

    #_____integrate the orbit with C_____
    ts= numpy.linspace(0,101,100)
    o.integrate(ts,pot,method='leapfrog_c')

    #_____Setup ActionAngle object and calculate actions (Staeckel approximation)_____
    aAS = actionAngleStaeckel(pot=pot,delta=Delta,c=True)
    jrs,lzs,jzs = aAS(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts))

    assert numpy.all(numpy.fabs(jrs - jrs[0]) < 10.**-8.), \
        'Radial action is not conserved along orbit.'

    assert numpy.all(numpy.fabs(lzs - lzs[0]) < 10.**-8.), \
        'Angular momentum is not conserved along orbit.'

    assert numpy.all(numpy.fabs(jzs - jzs[0]) < 10.**-8.), \
        'Vertical action is not conserved along orbit.'

    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vz [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 vz [as 别名]
def test_actionAngleTorus_orbit():
    from galpy.actionAngle import actionAngleTorus
    from galpy.potential import MWPotential2014
    from galpy.orbit import Orbit
    # Set up instance
    aAT= actionAngleTorus(pot=MWPotential2014,tol=10.**-5.)
    jr,jphi,jz= 0.05,1.1,0.025
    # First calculate frequencies and the initial RvR
    RvRom= aAT.xvFreqs(jr,jphi,jz,
                       numpy.array([0.]),
                       numpy.array([1.]),
                       numpy.array([2.]))
    om= RvRom[1:]
    # Angles along an orbit
    ts= numpy.linspace(0.,100.,1001)
    angler= ts*om[0]
    anglephi= 1.+ts*om[1]
    anglez= 2.+ts*om[2]
    # Calculate the orbit using actionAngleTorus
    RvR= aAT(jr,jphi,jz,angler,anglephi,anglez).T
    # Calculate the orbit using orbit integration
    orb= Orbit([RvRom[0][:,0],RvRom[0][:,1],RvRom[0][:,2],
                RvRom[0][:,3],RvRom[0][:,4],RvRom[0][:,5]])
    orb.integrate(ts,MWPotential2014)
    # Compare
    tol= -3.
    assert numpy.all(numpy.fabs(orb.R(ts)-RvR[0]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in R'
    assert numpy.all(numpy.fabs(orb.vR(ts)-RvR[1]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in vR'
    assert numpy.all(numpy.fabs(orb.vT(ts)-RvR[2]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in vT'
    assert numpy.all(numpy.fabs(orb.z(ts)-RvR[3]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in z'
    assert numpy.all(numpy.fabs(orb.vz(ts)-RvR[4]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in vz'
    assert numpy.all(numpy.fabs(orb.phi(ts)-RvR[5]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in phi'
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vz [as 别名]
def calc_progenitor_actions(savefilename):
    # Setup potential
    lp = potential.LogarithmicHaloPotential(normalize=1.0, q=0.9)
    # Setup orbit
    x, z, y, vx, vz, vy = -11.63337239, -10.631736273934635, -20.76235661, -128.8281653, 79.172383882274971, 42.88727925
    R, phi, z = bovy_coords.rect_to_cyl(x, y, z)
    vR, vT, vz = bovy_coords.rect_to_cyl_vec(vx, vy, vz, R, phi, z, cyl=True)
    R /= 8.0
    z /= 8.0
    vR /= 220.0
    vT /= 220.0
    vz /= 220.0
    o = Orbit([R, vR, vT, z, vz, phi])
    ts = numpy.linspace(0.0, 5.125 * 220.0 / 8.0, 1313)  # times of the snapshots
    o.integrate(ts, lp, method="dopr54_c")
    if "aas" in savefilename:
        aA = actionAngleStaeckel(pot=lp, delta=1.20, c=True)
    else:
        aA = actionAngleIsochroneApprox(pot=lp, b=0.8)
    # Now calculate actions, frequencies, and angles for all positions
    Rs = o.R(ts)
    vRs = o.vR(ts)
    vTs = o.vT(ts)
    zs = o.z(ts)
    vzs = o.vz(ts)
    phis = o.phi(ts)
    csvfile = open(savefilename, "wb")
    writer = csv.writer(csvfile, delimiter=",")
    for ii in range(len(ts)):
        acfs = aA.actionsFreqsAngles(Rs[ii], vRs[ii], vTs[ii], zs[ii], vzs[ii], phis[ii])
        writer.writerow(
            [acfs[0][0], acfs[1][0], acfs[2][0], acfs[3][0], acfs[4][0], acfs[5][0], acfs[6][0], acfs[7][0], acfs[8][0]]
        )
        csvfile.flush()
    csvfile.close()
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vz [as 别名]
def calcOrbits(parser):
    options,args= parser.parse_args()
    #Read data
    XYZ,vxvyvz,cov_vxvyvz,rawdata= readData(metal='allall',
                                            sample=options.sample,
                                            loggmin=4.2,
                                            snmin=15.,
                                            select=options.select)
    #Define potential
    if options.logp:
        pot= LogarithmicHaloPotential(normalize=1.)
    else:
        pot= MWPotential
    ts= numpy.linspace(0.,_MAXT,10000) #times to integrate
    if os.path.exists(args[0]):#Load savefile
        savefile= open(args[0],'rb')
        orbits= pickle.load(savefile)
        _ORBITSLOADED= True
        try:
            samples= pickle.load(savefile)
        except EORError:
            _SAMPLESLOADED= False
        else:
            _SAMPLESLOADED= True
        finally:
            savefile.close()
    else:
        _ORBITSLOADED= False
    if not _ORBITSLOADED:
        #First calculate orbits
        es, rmeans, rperis, raps, zmaxs = [], [], [], [], []
        densrmeans, vzrmeans= [], []
        for ii in range(len(rawdata)):
            sys.stdout.write('\r'+"Working on object %i/%i" % (ii,len(rawdata)))
            sys.stdout.flush()
            #Integrate the orbit
            data= rawdata[ii]
            o= Orbit([data.ra,data.dec,data.dist,data.pmra,data.pmdec,data.vr],
                     radec=True,vo=220.,ro=8.,zo=_ZSUN)
            o.integrate(ts,pot)
            es.append(o.e())
            rperis.append(o.rperi())
            raps.append(o.rap())
            zmaxs.append(o.zmax())
            rmeans.append(0.5*(o.rperi()+o.rap()))
            Rs= o.R(ts)
            vz2= o.vz(ts)**2.
            dens= evaluateDensities(Rs,0.*Rs,pot)
            densrmeans.append(numpy.sum(dens*Rs)/numpy.sum(dens))
            vzrmeans.append(numpy.sum(vz2*Rs)/numpy.sum(vz2))
#            print " ", rmeans[-1], densrmeans[-1], vzrmeans[-1]
        sys.stdout.write('\r'+_ERASESTR+'\r')
        sys.stdout.flush()
        es= numpy.array(es)
        rmeans= numpy.array(rmeans)
        rperis= numpy.array(rperis)
        raps= numpy.array(raps)
        zmaxs= numpy.array(zmaxs)
        orbits= _append_field_recarray(rawdata,'e',es)
        orbits= _append_field_recarray(orbits,'rmean',rmeans)
        orbits= _append_field_recarray(orbits,'rperi',rperis)
        orbits= _append_field_recarray(orbits,'rap',raps)
        orbits= _append_field_recarray(orbits,'zmax',zmaxs)
        orbits= _append_field_recarray(orbits,'densrmean',densrmeans)
        orbits= _append_field_recarray(orbits,'vzrmean',vzrmeans)
        #Pickle
        savefile= open(args[0],'wb')
        pickle.dump(orbits,savefile)
        savefile.close()
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vz [as 别名]
def plot_aaspher_conservation(plotfilename1,plotfilename2):
    #Setup orbit
    E, Lz= -1.25, 0.6
    o= Orbit([0.8,0.3,Lz/0.8,0.,numpy.sqrt(2.*(E-evalPot(0.8,0.,MWPotential2014)-(Lz/0.8)**2./2.-0.3**2./2.)),0.])
    #Integrate the orbit to estimate an equivalent b
    nt= 1001
    ts= numpy.linspace(0.,20.,nt)
    o.integrate(ts,MWPotential2014,method='symplec4_c')
    b= estimateBIsochrone(o.R(ts),o.z(ts),pot=MWPotential2014)
    print b
    b= 0.3
    #Now integrate the orbit in the isochronePotential
    ip= IsochronePotential(normalize=1.,b=b)
    aAI= actionAngleIsochrone(ip=ip)
    orbt= 2.*numpy.pi/aAI.actionsFreqs(o)[4]    
    norb= 200.
    nt= 20001
    ts= numpy.linspace(0.,norb*orbt,nt)
    o.integrate(ts,ip,method='symplec4_c')
    #Calculate actions, frequencies, and angles
    jfa= aAI.actionsFreqsAngles(o.R(ts),o.vR(ts),o.vT(ts),
                                o.z(ts),o.vz(ts),o.phi(ts))
    dJs= numpy.fabs((jfa[0]-numpy.mean(jfa[0]))/numpy.mean(jfa[0]))
    dOrs= numpy.fabs((jfa[3]-numpy.mean(jfa[3]))/numpy.mean(jfa[3]))
    dOzs= numpy.fabs((jfa[5]-numpy.mean(jfa[5]))/numpy.mean(jfa[5]))
    print "frequencies", numpy.mean(dOrs), numpy.mean(dOzs)
    ar= dePeriod(numpy.reshape(jfa[6],(1,len(ts)))).flatten()
    az= dePeriod(numpy.reshape(jfa[8],(1,len(ts)))).flatten()
    danglers= numpy.fabs(ar-numpy.mean(jfa[3])*ts-jfa[6][0])/2./numpy.pi
    danglezs= numpy.fabs(az-numpy.mean(jfa[5])*ts-jfa[8][0])/2./numpy.pi
    #Break up
    breakt= 50.
    pts= parse_break(ts,ts < breakt)
    pdJs= parse_break(dJs,ts < breakt)
    pdanglers= parse_break(danglers,ts < breakt)
    pdanglezs= parse_break(danglezs,ts < breakt)
    #dAngles
    bovy_plot.bovy_print()
    pyplot.subplot(2,1,1)
    bovy_plot.bovy_plot(pts/orbt,
                        pdJs,
                        color='k',loglog=True,gcf=True,
                        xrange=[0.5,norb],
                        yrange=[10.**-12.,1.])
    bovy_plot.bovy_text(r'$\texttt{actionAngleIsochrone}$',
                        top_left=True,size=14.)
    ax= pyplot.gca()
    ax.yaxis.set_ticks([10.**-12.,10.**-8.,10.**-4.,1.])
    nullfmt   = NullFormatter()         # no labels
    ax.xaxis.set_major_formatter(nullfmt)
    #Same for actionAngleSpherical
    aAS= actionAngleSpherical(pot=ip)
    tts= ts[::1]
    jfa= aAS.actionsFreqsAngles(o.R(tts),o.vR(tts),o.vT(tts),
                                o.z(tts),o.vz(tts),o.phi(tts),
                                fixed_quad=True)
    #dJr
    dJs= numpy.fabs((jfa[0]-numpy.mean(jfa[0]))/numpy.mean(jfa[0]))
    dOrs= numpy.fabs((jfa[3]-numpy.mean(jfa[3]))/numpy.mean(jfa[3]))
    dOzs= numpy.fabs((jfa[5]-numpy.mean(jfa[5]))/numpy.mean(jfa[5]))
    print "frequencies", numpy.mean(dOrs), numpy.mean(dOzs)
    #dAngles
    ar= dePeriod(numpy.reshape(jfa[6],(1,len(tts)))).flatten()
    az= dePeriod(numpy.reshape(jfa[8],(1,len(tts)))).flatten()
    danglers= numpy.fabs(ar-numpy.mean(jfa[3])*tts-jfa[6][0])/2./numpy.pi
    danglezs= numpy.fabs(az-numpy.mean(jfa[5])*tts-jfa[8][0])/2./numpy.pi
    print numpy.mean(danglers)
    print numpy.mean(danglezs)
    ptts= parse_break(tts,tts < breakt)
    pdJs= parse_break(dJs,tts < breakt)
    pyplot.subplot(2,1,2)
    bovy_plot.bovy_plot(ptts/orbt,
                        pdJs,
                        color='k',loglog=True,gcf=True,
                        xrange=[0.5,norb],
                        yrange=[10.**-12.,1.],
                        xlabel=r'$\mathrm{Number\ of\ orbital\ periods}$')
    bovy_plot.bovy_text(r'$\texttt{actionAngleSpherical}$',
                        top_left=True,size=14.)
    bovy_plot.bovy_text(0.175,10.**2.,r'$\left|\Delta J_R / J_R\right|$',
                        fontsize=16.,
                        rotation='vertical')
    ax= pyplot.gca()
    ax.xaxis.set_major_formatter(ticker.FormatStrFormatter(r'$%0.f$'))
    ax.yaxis.set_ticks([10.**-12.,10.**-8.,10.**-4.,1.])
    bovy_plot.bovy_end_print(plotfilename1)    
    #Now plot the deviations in the angles
    bovy_plot.bovy_print()
    pyplot.subplot(2,1,1)
    liner= bovy_plot.bovy_plot(pts/orbt,
                               pdanglers,
                               color='k',ls='-',loglog=True,gcf=True,
                               xrange=[0.5,norb],
                               yrange=[10.**-12.,1.])
    linez= bovy_plot.bovy_plot(pts/orbt,
                               pdanglezs,
                               color='k',ls='--',overplot=True)
    legend1= pyplot.legend((liner[0],linez[0]),
                           (r'$\theta_R$',
                            r'$\theta_z$'),
#.........这里部分代码省略.........

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

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

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