Python Orbit.z方法代碼示例

# 需要導入模塊: from galpy.orbit import Orbit [as 別名]
# 或者: from galpy.orbit.Orbit import z [as 別名]
def test_surfacesection():
    #Preliminary code
    import numpy
    from galpy.potential import MWPotential2014
    from galpy.potential import evaluatePotentials as evalPot
    from galpy.orbit import Orbit
    E, Lz= -1.25, 0.6
    o1= Orbit([0.8,0.,Lz/0.8,0.,numpy.sqrt(2.*(E-evalPot(0.8,0.,MWPotential2014)-(Lz/0.8)**2./2.)),0.])
    ts= numpy.linspace(0.,100.,2001)
    o1.integrate(ts,MWPotential2014)
    o2= 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.])
    o2.integrate(ts,MWPotential2014)
    def surface_section(Rs,zs,vRs):
        # Find points where the orbit crosses z from - to +
        shiftzs= numpy.roll(zs,-1)
        indx= (zs[:-1] < 0.)*(shiftzs[:-1] > 0.)
        return (Rs[:-1][indx],vRs[:-1][indx])
    # Calculate and plot the surface of section
    ts= numpy.linspace(0.,1000.,20001) # long integration
    o1.integrate(ts,MWPotential2014)
    o2.integrate(ts,MWPotential2014)
    sect1Rs,sect1vRs=surface_section(o1.R(ts),o1.z(ts),o1.vR(ts))
    sect2Rs,sect2vRs=surface_section(o2.R(ts),o2.z(ts),o2.vR(ts))
    from matplotlib.pyplot import plot, xlim, ylim
    plot(sect1Rs,sect1vRs,'bo',mec='none')
    xlim(0.3,1.); ylim(-0.69,0.69)
    plot(sect2Rs,sect2vRs,'yo',mec='none')
    return None

# 需要導入模塊: from galpy.orbit import Orbit [as 別名]
# 或者: from galpy.orbit.Orbit import z [as 別名]
def test_estimateDelta():

    #_____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')


    #____estimate Focal length Delta_____
    #for each time step individually:
    deltas_estimate = numpy.zeros(len(ts))
    for ii in range(len(ts)):
        deltas_estimate[ii] = estimateDeltaStaeckel(pot,o.R(ts[ii]),o.z(ts[ii]))

    assert numpy.all(numpy.fabs(deltas_estimate - Delta) < 10.**-8), \
            'Focal length Delta estimated along the orbit is not constant.'

    #for all time steps together:
    delta_estimate = estimateDeltaStaeckel(pot,o.R(ts),o.z(ts))
    
    assert numpy.fabs(delta_estimate - Delta) < 10.**-8, \
            'Focal length Delta estimated from the orbit is not the same as the input focal length.'

    return None

# 需要導入模塊: from galpy.orbit import Orbit [as 別名]
# 或者: from galpy.orbit.Orbit import z [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 z [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 z [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 z [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 z [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 z [as 別名]
def calc_delta_MWPotential2014(savefilename,plotfilename):
    Lmin, Lmax= 0.01, 10.
    if not os.path.exists(savefilename):
        #Setup grid
        nL, nE= 101,101
        Ls= 10.**numpy.linspace(numpy.log10(Lmin),numpy.log10(Lmax),nL)
        #Integration times
        ts= numpy.linspace(0.,20.,1001)
        deltas= numpy.empty((nL,nE))
        Einf= evalPot(10.**12.,0.,MWPotential2014)
        print Einf
        for ii in range(nL):
            #Calculate Ec
            Rc= rl(MWPotential2014,Ls[ii])
            print ii, "Rc = ", Rc*8.
            Ec= evalPot(Rc,0.,MWPotential2014)+0.5*Ls[ii]**2./Rc**2.
            Es= Ec+(Einf-Ec)*10.**numpy.linspace(-2.,0.,nE)
            for jj in range(nE):
                #Setup an orbit with this energy and angular momentum
                Er= 2.*(Es[jj]-Ec) #Random energy times 2 = vR^2 + vz^2
                vR= numpy.sqrt(4./5.*Er)
                vz= numpy.sqrt(1./5.*Er)
                o= Orbit([Rc,vR,Ls[ii]/Rc,0.,vz,0.])
                o.integrate(ts,MWPotential2014,method='symplec4_c')
                deltas[ii,jj]= estimateDeltaStaeckel(o.R(ts),o.z(ts),
                                                     pot=MWPotential2014)
        #Save
        save_pickles(savefilename,deltas)
    else:
        savefile= open(savefilename,'rb')
        deltas= pickle.load(savefile)
        savefile.close()
    #Plot
    print numpy.nanmax(deltas)
    bovy_plot.bovy_print()
    bovy_plot.bovy_dens2d(deltas.T,origin='lower',cmap='jet',
                          xrange=[numpy.log10(Lmin),numpy.log10(Lmax)],
                          yrange=[-2,0.],
                          xlabel=r'$\log_{10} L$',
                          ylabel=r'$\log_{10}\left(\frac{E-E_c(L)}{E(\infty)-E_c(L)}\right)$',
                          colorbar=True,shrink=0.78,
                          zlabel=r'$\mathrm{Approximate\ focal\ length}$',
#                          interpolation='nearest',
                          vmin=0.,vmax=0.6)
    bovy_plot.bovy_end_print(plotfilename)
    return None

# 需要導入模塊: from galpy.orbit import Orbit [as 別名]
# 或者: from galpy.orbit.Orbit import z [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 z [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 z [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 z [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 z [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 z [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 z [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$'),
#.........這裏部分代碼省略.........