Python Orbit.e方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import e [as 别名]
def test_orbmethods():
    from galpy.orbit import Orbit
    from galpy.potential import MWPotential2014
    o= Orbit([0.8,0.3,0.75,0.,0.2,0.]) # setup R,vR,vT,z,vz,phi
    times= numpy.linspace(0.,10.,1001) # Output times
    o.integrate(times,MWPotential2014) # Integrate
    o.E() # Energy
    assert numpy.fabs(o.E()+1.2547650648697966) < 10.**-5., 'Orbit method does not work as expected'
    o.L() # Angular momentum
    assert numpy.all(numpy.fabs(o.L()-numpy.array([[ 0.  , -0.16,  0.6 ]])) < 10.**-5.), 'Orbit method does not work as expected'
    o.Jacobi(OmegaP=0.65) #Jacobi integral E-OmegaP Lz
    assert numpy.fabs(o.Jacobi(OmegaP=0.65)-numpy.array([-1.64476506])) < 10.**-5., 'Orbit method does not work as expected'
    o.ER(times[-1]), o.Ez(times[-1]) # Rad. and vert. E at end
    assert numpy.fabs(o.ER(times[-1])+1.27601734263047) < 10.**-5., 'Orbit method does not work as expected'
    assert numpy.fabs(o.Ez(times[-1])-0.021252201847851909) < 10.**-5.,  'Orbit method does not work as expected'
    o.rperi(), o.rap(), o.zmax() # Peri-/apocenter r, max. |z|
    assert numpy.fabs(o.rperi()-0.44231993168097) < 10.**-5., 'Orbit method does not work as expected'
    assert numpy.fabs(o.rap()-0.87769030382105) < 10.**-5., 'Orbit method does not work as expected'
    assert numpy.fabs(o.zmax()-0.077452357289016) < 10.**-5., 'Orbit method does not work as expected'
    o.e() # eccentricity (rap-rperi)/(rap+rperi)
    assert numpy.fabs(o.e()-0.32982348199330563) < 10.**-5., 'Orbit method does not work as expected'
    o.R(2.,ro=8.) # Cylindrical radius at time 2. in kpc
    assert numpy.fabs(o.R(2.,ro=8.)-3.5470772876920007) < 10.**-3., 'Orbit method does not work as expected'
    o.vR(5.,vo=220.) # Cyl. rad. velocity at time 5. in km/s
    assert numpy.fabs(o.vR(5.,vo=220.)-45.202530965094553) < 10.**-3., 'Orbit method does not work as expected'
    o.ra(1.), o.dec(1.) # RA and Dec at t=1. (default settings)
    # 5/12/2016: test weakened, because improved galcen<->heliocen 
    #            transformation has changed these, but still close
    assert numpy.fabs(o.ra(1.)-numpy.array([ 288.19277])) < 10.**-1., 'Orbit method does not work as expected'
    assert numpy.fabs(o.dec(1.)-numpy.array([ 18.98069155])) < 10.**-1., 'Orbit method does not work as expected'
    o.jr(type='adiabatic'), o.jz() # R/z actions (ad. approx.)
    assert numpy.fabs(o.jr(type='adiabatic')-0.05285302231137586) < 10.**-3., 'Orbit method does not work as expected'
    assert numpy.fabs(o.jz()-0.006637988850751242) < 10.**-3., 'Orbit method does not work as expected'
    # Rad. period w/ Staeckel approximation w/ focal length 0.5,
    o.Tr(type='staeckel',delta=0.5,ro=8.,vo=220.) # in Gyr  
    assert numpy.fabs(o.Tr(type='staeckel',delta=0.5,ro=8.,vo=220.)-0.1039467864018446) < 10.**-3., 'Orbit method does not work as expected'
    o.plot(d1='R',d2='z') # Plot the orbit in (R,z)
    o.plot3d() # Plot the orbit in 3D, w/ default [x,y,z]
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import e [as 别名]
def illustrate_adiabatic_invariance(plotfilename1,plotfilename2):
    # Initialize two different IsochronePotentials
    ip1= IsochronePotential(normalize=1.,b=1.)
    ip2= IsochronePotential(normalize=0.5,b=1.)
    # Use TimeInterpPotential to interpolate smoothly between the two
    tip= TimeInterpPotential(ip1,ip2,dt=100.,tform=50.)
    # Integrate the orbit, in three parts
    # 1) Orbit in the first isochrone potential
    o1= Orbit([1.,0.1,1.1,0.0,0.1,0.])
    ts= numpy.linspace(0.,50.,1001)
    o1.integrate(ts,tip)
    bovy_plot.bovy_print()
    o1.plot(d1='x',d2='y',xrange=[-1.6,1.6],yrange=[-1.6,1.6],color='b',
            gcf=True)
    # 2) Orbit in the transition
    o2= o1(ts[-1]) # Last time step = initial time step of the next integration
    ts2= numpy.linspace(50.,150.,1001)
    o2.integrate(ts2,tip)
    o2.plot(d1='x',d2='y',overplot=True,color='g')
    # 3) Orbit in the second isochrone potential
    o3= o2(ts2[-1])
    ts3= numpy.linspace(150.,200.,1001)
    o3.integrate(ts3,ip2)
    o3.plot(d1='x',d2='y',overplot=True,color='r')
    bovy_plot.bovy_end_print(plotfilename1)
    # Also plot the R,z projection
    bovy_plot.bovy_print(fig_height=2.3333)
    o1.plot(d1='R',d2='z',xrange=[0.9,1.65],yrange=[-.175,.175],color='b',
            gcf=True)
    o2.plot(d1='R',d2='z',overplot=True,color='g')
    o3.plot(d1='R',d2='z',overplot=True,color='r')
    bovy_plot.bovy_end_print(plotfilename2)   
    # Now we calculate the energy, eccentricity, mean radius, and maximum height
    print o1.E(pot=ip1), o1.e(), 0.5*(o1.rperi()+o1.rap()), o1.zmax()
    print o3.E(pot=ip2), o3.e(), 0.5*(o3.rperi()+o3.rap()), o3.zmax()
   # The orbit has clearly moved to larger radii, the actions are however conserved
    aAI1= actionAngleIsochrone(ip=ip1)
    aAI2= actionAngleIsochrone(ip=ip2)
    print aAI1(o1)
    print aAI2(o3)
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import e [as 别名]
def calc_eccentricity(args, options):
    table = os.path.join(args[0],'table2.dat')
    readme = os.path.join(args[0],'ReadMe')
    dierickx = ascii.read(table, readme=readme)
    vxvv = np.dstack([dierickx['RAdeg'], dierickx['DEdeg'], dierickx['Dist']/1e3, dierickx['pmRA'], dierickx['pmDE'], dierickx['HRV']])[0]
    ro, vo, zo = 8., 220., 0.025
    ra, dec= vxvv[:,0], vxvv[:,1]
    lb= bovy_coords.radec_to_lb(ra,dec,degree=True)
    pmra, pmdec= vxvv[:,3], vxvv[:,4]
    pmllpmbb= bovy_coords.pmrapmdec_to_pmllpmbb(pmra,pmdec,ra,dec,degree=True)
    d, vlos= vxvv[:,2], vxvv[:,5]
    rectgal= bovy_coords.sphergal_to_rectgal(lb[:,0],lb[:,1],d,vlos,pmllpmbb[:,0], pmllpmbb[:,1],degree=True)
    vsolar= np.array([-10.1,4.0,6.7])
    vsun= np.array([0.,1.,0.,])+vsolar/vo
    X = rectgal[:,0]/ro
    Y = rectgal[:,1]/ro
    Z = rectgal[:,2]/ro
    vx = rectgal[:,3]/vo
    vy = rectgal[:,4]/vo
    vz = rectgal[:,5]/vo
    vsun= np.array([0.,1.,0.,])+vsolar/vo
    Rphiz= bovy_coords.XYZ_to_galcencyl(X,Y,Z,Zsun=zo/ro)
    vRvTvz= bovy_coords.vxvyvz_to_galcencyl(vx,vy,vz,Rphiz[:,0],Rphiz[:,1],Rphiz[:,2],vsun=vsun,Xsun=1.,Zsun=zo/ro,galcen=True)
    #do the integration and individual analytic estimate for each object
    ts= np.linspace(0.,20.,10000)
    lp= LogarithmicHaloPotential(normalize=1.)
    e_ana = numpy.zeros(len(vxvv))
    e_int = numpy.zeros(len(vxvv))
    print('Performing orbit integration and analytic parameter estimates for Dierickx et al. sample...')
    for i in tqdm(range(len(vxvv))):
        try:
            orbit = Orbit(vxvv[i], radec=True, vo=220., ro=8.)
            e_ana[i] = orbit.e(analytic=True, pot=lp, c=True)
        except UnboundError:
            e_ana[i] = np.nan
        orbit.integrate(ts, lp)
        e_int[i] = orbit.e(analytic=False)
    fig = plt.figure()
    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
    plt.scatter(e_int, e_ana,  s=1, color='Black', lw=0.)
    plt.xlabel(r'$\mathrm{galpy\ integrated}\ e$')
    plt.ylabel(r'$\mathrm{galpy\ analytic}\ e$')
    plt.xlim(0.,1.)
    plt.ylim(0.,1.)
    fig.tight_layout()
    plt.savefig(os.path.join(args[0],'dierickx-integratedeanalytice.png'), format='png', dpi=200)
    fig = plt.figure()
    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
    plt.hist(e_int, bins=30)
    plt.xlim(0.,1.)
    plt.xlabel(r'$\mathrm{galpy}\ e$')
    fig.tight_layout()
    plt.savefig(os.path.join(args[0], 'dierickx-integratedehist.png'), format='png', dpi=200)
    fig = plt.figure()
    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
    plt.scatter(dierickx['e'], e_int,  s=1, color='Black', lw=0.)
    plt.xlabel(r'$\mathrm{Dierickx\ et\ al.}\ e$')
    plt.ylabel(r'$\mathrm{galpy\ integrated}\ e$')
    plt.xlim(0.,1.)
    plt.ylim(0.,1.)
    fig.tight_layout()
    plt.savefig(os.path.join(args[0],'dierickx-integratedee.png'), format='png', dpi=200)
    fig = plt.figure()
    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
    plt.scatter(dierickx['e'], e_ana,  s=1, color='Black', lw=0.)
    plt.xlabel(r'$\mathrm{Dierickx\ et\ al.}\ e$')
    plt.ylabel(r'$\mathrm{galpy\ estimated}\ e$')
    plt.xlim(0.,1.)
    plt.ylim(0.,1.)
    fig.tight_layout()
    plt.savefig(os.path.join(args[0],'dierickx-analyticee.png'), format='png', dpi=200)
    arr = numpy.recarray(len(e_ana), dtype=[('analytic_e', float), ('integrated_e', float)])
    arr['analytic_e'] = e_ana
    arr['integrated_e'] = e_int
    with open(os.path.join(args[0],'eccentricities.dat'), 'w') as file:
        pickle.dump(arr, file)
        file.close()

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import e [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 e [as 别名]
def calc_es():
    savefilename= 'myes.sav'
    if os.path.exists(savefilename):
        savefile= open(savefilename,'rb')
        mye= pickle.load(savefile)
        e= pickle.load(savefile)
        savefile.close()
    else:
       #Read data
        dialect= csv.excel
        dialect.skipinitialspace=True
        reader= csv.reader(open('../data/Dierickx-etal-tab2.txt','r'),delimiter=' ',dialect=dialect)
        vxvs= []
        es= []
        vphis= []
        vxs= []
        vys= []
        vzs= []
        ls= []
        for row in reader:
            thisra= float(row[3])
            thisdec= float(row[4])
            thisd= float(row[17])/1000.
            thispmra= float(row[13])
            thispmdec= float(row[15])
            thisvlos= float(row[11])
            thise= float(row[26])
            vxvs.append([thisra,thisdec,thisd,thispmra,thispmdec,thisvlos])
            es.append(thise)
            vphis.append(float(row[25]))
            vxs.append(float(row[19]))
            vys.append(float(row[21]))
            vzs.append(float(row[23]))
            ls.append(float(row[5]))
        vxvv= nu.array(vxvs)
        e= nu.array(es)
        vphi= nu.array(vphis)
        vx= nu.array(vxs)
        vy= nu.array(vys)
        vz= nu.array(vzs)
        l= nu.array(ls)

        #Define potential
        lp= LogarithmicHaloPotential(normalize=1.)
        mp= MiyamotoNagaiPotential(a=0.5,b=0.0375,amp=1.,normalize=.6)
        np= NFWPotential(a=4.5,normalize=.35)
        hp= HernquistPotential(a=0.6/8,normalize=0.05)
        ts= nu.linspace(0.,20.,10000)
        
        mye= nu.zeros(len(e))
        for ii in range(len(e)):
           #Integrate the orbit
            o= Orbit(vxvv[ii,:],radec=True,vo=220.,ro=8.)
            o.integrate(ts,lp)
            mye[ii]= o.e()
            

        #Save
        savefile= open(savefilename,'wb')
        pickle.dump(mye,savefile)
        pickle.dump(e,savefile)
        savefile.close()

    #plot
    plot.bovy_print()
    plot.bovy_plot(nu.array([0.,1.]),nu.array([0.,1.]),'k-',
                   xlabel=r'$\mathrm{Dierickx\ et\ al.}\ e$',
                   ylabel=r'$\mathrm{galpy}\ e$')
    plot.bovy_plot(e,mye,'k,',overplot=True)
    plot.bovy_end_print('myee.png')

    plot.bovy_print()
    plot.bovy_hist(e,bins=30,xlabel=r'$\mathrm{Dierickx\ et\ al.}\ e$')
    plot.bovy_end_print('ehist.png')

    plot.bovy_print()
    plot.bovy_hist(mye,bins=30,xlabel=r'$\mathrm{galpy}\ e$')
    plot.bovy_end_print('myehist.png')

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import e [as 别名]
def calcj(rotcurve):
    if rotcurve == 'flat':
        savefilename= 'myjs.sav'
    elif rotcurve == 'power':
        savefilename= 'myjs_power.sav'
    if os.path.exists(savefilename):
        savefile= open(savefilename,'rb')
        myjr= pickle.load(savefile)
        myjp= pickle.load(savefile)
        mywr= pickle.load(savefile)
        mywp= pickle.load(savefile)
        mye= pickle.load(savefile)
        myzmax= pickle.load(savefile)
        e= pickle.load(savefile)
        zmax= pickle.load(savefile)
        savefile.close()
    else:
        dialect= csv.excel
        dialect.skipinitialspace=True
        reader= csv.reader(open('../data/gcs.tsv','r'),delimiter='|',dialect=dialect)
        vxvs= []
        es= []
        zmaxs= []
        for row in reader:
            if row[0][0] == '#':
                continue
            thisra= row[0]
            thisdec= row[1]
            thisd= read_float(row[2])/1000.
            if thisd > 0.2: continue
            thisu= read_float(row[3])
            thisv= read_float(row[4])
            thisw= read_float(row[5])
            thise= read_float(row[6])
            thiszmax= read_float(row[7])
            if thisd == -9999 or thisu == -9999 or thisv == -9999 or thisw == -9999:
                continue
            vxvs.append([hms_to_rad(thisra),dms_to_rad(thisdec),
                         thisd,thisu,thisv,thisw])
            es.append(thise)
            zmaxs.append(thiszmax)
        vxvv= nu.array(vxvs)
        e= nu.array(es)
        zmax= nu.array(zmaxs)

        #Define potential
        lp= LogarithmicHaloPotential(normalize=1.)
        pp= PowerSphericalPotential(normalize=1.,alpha=-2.)
        mp= MiyamotoNagaiPotential(a=0.5,b=0.0375,amp=1.,normalize=.6)
        np= NFWPotential(a=4.5,normalize=.35)
        hp= HernquistPotential(a=0.6/8,normalize=0.05)
        ts= nu.linspace(0.,20.,10000)

        myjr= nu.zeros(len(e))
        myjp= nu.zeros(len(e))
        mywr= nu.zeros(len(e))
        mywp= nu.zeros(len(e))
        mye= nu.zeros(len(e))
        myzmax= nu.zeros(len(e))
        for ii in range(len(e)):
           #Integrate the orbit
            o= Orbit(vxvv[ii,:],radec=True,uvw=True,vo=220.,ro=8.)
            #o.integrate(ts,[mp,np,hp])
            if rotcurve == 'flat':
                o.integrate(ts,lp)
                mye[ii]= o.e()
                myzmax[ii]= o.zmax()*8.
                print e[ii], mye[ii], zmax[ii], myzmax[ii]
                o= o.toPlanar()
                myjr[ii]= o.jr(lp)[0]
            else:
                o= o.toPlanar()
                myjr[ii]= o.jr(pp)[0]
            myjp[ii]= o.jp()[0]
            mywr[ii]= o.wr()[0]
            mywp[ii]= o.wp()

        #Save
        savefile= open(savefilename,'wb')
        pickle.dump(myjr,savefile)
        pickle.dump(myjp,savefile)
        pickle.dump(mywr,savefile)
        pickle.dump(mywp,savefile)
        pickle.dump(mye,savefile)
        pickle.dump(myzmax,savefile)
        pickle.dump(e,savefile)
        pickle.dump(zmax,savefile)
        savefile.close()

    #plot
    if rotcurve == 'flat':
        plot.bovy_print()
        plot.bovy_plot(nu.array([0.,1.]),nu.array([0.,1.]),'k-',
                       xlabel=r'$\mathrm{Holmberg\ et\ al.}\ e$',
                       ylabel=r'$\mathrm{galpy}\ e$')
        plot.bovy_plot(e,mye,'k,',overplot=True)
        plot.bovy_end_print('myee.png')
        
        plot.bovy_print()
        plot.bovy_plot(nu.array([0.,2.5]),
#.........这里部分代码省略.........

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import e [as 别名]
def calcj(rotcurve):
    if rotcurve == "flat":
        savefilename = "myjs.sav"
    elif rotcurve == "power":
        savefilename = "myjs_power.sav"
    if os.path.exists(savefilename):
        savefile = open(savefilename, "rb")
        myjr = pickle.load(savefile)
        myjp = pickle.load(savefile)
        mye = pickle.load(savefile)
        myzmax = pickle.load(savefile)
        e = pickle.load(savefile)
        zmax = pickle.load(savefile)
        savefile.close()
    else:
        dialect = csv.excel
        dialect.skipinitialspace = True
        reader = csv.reader(open("../data/gcs.tsv", "r"), delimiter="|", dialect=dialect)
        vxvs = []
        es = []
        zmaxs = []
        for row in reader:
            if row[0][0] == "#":
                continue
            thisra = row[0]
            thisdec = row[1]
            thisd = read_float(row[2]) / 1000.0
            if thisd > 0.2:
                continue
            thisu = read_float(row[3])
            thisv = read_float(row[4])
            thisw = read_float(row[5])
            thise = read_float(row[6])
            thiszmax = read_float(row[7])
            if thisd == -9999 or thisu == -9999 or thisv == -9999 or thisw == -9999:
                continue
            vxvs.append([hms_to_rad(thisra), dms_to_rad(thisdec), thisd, thisu, thisv, thisw])
            es.append(thise)
            zmaxs.append(thiszmax)
        vxvv = nu.array(vxvs)
        e = nu.array(es)
        zmax = nu.array(zmaxs)

        # Define potential
        lp = LogarithmicHaloPotential(normalize=1.0)
        pp = PowerSphericalPotential(normalize=1.0, alpha=-2.0)
        ts = nu.linspace(0.0, 100.0, 10000)

        myjr = nu.zeros(len(e))
        myjp = nu.zeros(len(e))
        mye = nu.zeros(len(e))
        myzmax = nu.zeros(len(e))
        for ii in range(len(e)):
            # Integrate the orbit
            o = Orbit(vxvv[ii, :], radec=True, uvw=True, vo=220.0, ro=8.0)
            if rotcurve == "flat":
                o.integrate(ts, lp)
                mye[ii] = o.e()
                myzmax[ii] = o.zmax() * 8.0
                print e[ii], mye[ii], zmax[ii], myzmax[ii]
                myjr[ii] = o.jr(lp)
            else:
                myjr[ii] = o.jr(pp)
            myjp[ii] = o.jp()

        # Save
        savefile = open(savefilename, "wb")
        pickle.dump(myjr, savefile)
        pickle.dump(myjp, savefile)
        pickle.dump(mye, savefile)
        pickle.dump(myzmax, savefile)
        pickle.dump(e, savefile)
        pickle.dump(zmax, savefile)
        savefile.close()

    # plot
    if rotcurve == "flat":
        plot.bovy_print()
        plot.bovy_plot(
            nu.array([0.0, 1.0]),
            nu.array([0.0, 1.0]),
            "k-",
            xlabel=r"$\mathrm{Holmberg\ et\ al.}\ e$",
            ylabel=r"$\mathrm{galpy}\ e$",
        )
        plot.bovy_plot(e, mye, "k,", overplot=True)
        plot.bovy_end_print("myee.png")

        plot.bovy_print()
        plot.bovy_plot(
            nu.array([0.0, 2.5]),
            nu.array([0.0, 2.5]),
            "k-",
            xlabel=r"$\mathrm{Holmberg\ et\ al.}\ z_{\mathrm{max}}$",
            ylabel=r"$\mathrm{galpy}\ z_{\mathrm{max}}$",
        )
        plot.bovy_plot(zmax, myzmax, "k,", overplot=True)
        plot.bovy_end_print("myzmaxzmax.png")

    plot.bovy_print()
#.........这里部分代码省略.........