Python Orbit.rperi方法代码示例

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