Python Orbit.zmax方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import zmax [as 别名]
def test_adinvariance():
    from galpy.potential import IsochronePotential
    from galpy.orbit import Orbit
    from galpy.actionAngle import actionAngleIsochrone
    # Initialize two different IsochronePotentials
    ip1= IsochronePotential(normalize=1.,b=1.)
    ip2= IsochronePotential(normalize=0.5,b=1.)
    # Use TimeInterpPotential to interpolate smoothly
    tip= TimeInterpPotential(ip1,ip2,dt=100.,tform=50.)
    # Integrate: 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)
    o1.plot(d1='x',d2='y',xrange=[-1.6,1.6],yrange=[-1.6,1.6],
            color='b')
    # 2) Orbit in the transition
    o2= o1(ts[-1]) # Last time step => initial time step
    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,tip)
    o3.plot(d1='x',d2='y',overplot=True,color='r')
    # Now we calculate energy, maximum height, and mean radius
    print(o1.E(pot=ip1), (o1.rperi()+o1.rap())/2, o1.zmax())
    assert numpy.fabs(o1.E(pot=ip1)+2.79921356237) < 10.**-4., 'Energy in the adiabatic invariance test is different'
    assert numpy.fabs((o1.rperi()+o1.rap())/2-1.07854158141) < 10.**-4., 'mean radius in the adiabatic invariance test is different'
    assert numpy.fabs(o1.zmax()-0.106331362938) < 10.**-4., 'zmax in the adiabatic invariance test is different'
    print(o3.E(pot=ip2), (o3.rperi()+o3.rap())/2, o3.zmax())
    assert numpy.fabs(o3.E(pot=ip2)+1.19677002624) < 10.**-4., 'Energy in the adiabatic invariance test is different'
    assert numpy.fabs((o3.rperi()+o3.rap())/2-1.39962036137) < 10.**-4., 'mean radius in the adiabatic invariance test is different'
    assert numpy.fabs(o3.zmax()-0.138364269321) < 10.**-4., 'zmax in the adiabatic invariance test is different'
    # The orbit has clearly moved to larger radii,
    # the actions are however conserved from beginning to end
    aAI1= actionAngleIsochrone(ip=ip1); print(aAI1(o1))
    js= aAI1(o1)
    assert numpy.fabs(js[0]-numpy.array([ 0.00773779])) < 10.**-4., 'action in the adiabatic invariance test is different'
    assert numpy.fabs(js[1]-numpy.array([ 1.1])) < 10.**-4., 'action in the adiabatic invariance test is different'
    assert numpy.fabs(js[2]-numpy.array([ 0.0045361])) < 10.**-4., 'action in the adiabatic invariance test is different'
    aAI2= actionAngleIsochrone(ip=ip2); print(aAI2(o3))  
    js= aAI2(o3)
    assert numpy.fabs(js[0]-numpy.array([ 0.00773812])) < 10.**-4., 'action in the adiabatic invariance test is different'
    assert numpy.fabs(js[1]-numpy.array([ 1.1])) < 10.**-4., 'action in the adiabatic invariance test is different'
    assert numpy.fabs(js[2]-numpy.array([ 0.0045361])) < 10.**-4., 'action in the adiabatic invariance test is different'
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import zmax [as 别名]
def plot_twoorbits(plotfilename1,plotfilename2):
    #First 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)
    print "First orbit: E, L = %f,%f" % (o1.E(),o1.L()[:,2])
    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)
    print "Second orbit: E, L = %f,%f" % (o2.E(),o2.L()[:,2])
    print "First orbit: zmax = %f" % (o1.zmax())
    print "Second orbit: zmax = %f" % (o2.zmax())
    o1.plot(xrange=[0.3,1.],yrange=[-0.2,0.2],color='k')
    bovy_plot.bovy_end_print(plotfilename1)
    o2.plot(xrange=[0.3,1.],yrange=[-0.2,0.2],color='k')
    bovy_plot.bovy_end_print(plotfilename2)
    return (o1,o2)

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