Python Orbit.plot方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import plot [as 别名]
def test_orbitint():
    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)
    o1.plot(xrange=[0.3,1.],yrange=[-0.2,0.2],color='k')
    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)
    o2.plot(xrange=[0.3,1.],yrange=[-0.2,0.2],color='k')
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import plot [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 plot [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 plot [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 plot [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 plot [as 别名]
def demo_lsr_and_sun_cal():
    """
    Litte demo showing how one would calculate the orbit of the LSR and sun
    :return:
    """
    perimeter = 2 * np.pi * 8 * u.kpc
    velocity  = 220 * u.km/ u.s
    # for reference, LSR (at 8 kpc, with V = 220 km/s) should take this long
    # to complete one orbit
    orbit_time = (perimeter / velocity).to("Gyr")

    max_age = 100 * orbit_time / bovy_conversion.time_in_Gyr(220., 8.) # Gyr
    ntimes = 10000
    ts = np.linspace(0,max_age,ntimes)

    # INITIALISING SUN COORDINATES AND ORBIT
    #deg, deg, kpc,  mas/yr, mas/yr, km/s
    ra, dec, dist, mu_ra,  mu_dec, vlos = 0., 0., 0., 0., 0., 0.
    solar_coords = [ra, dec, dist, mu_ra, mu_dec, vlos]
    sun = Orbit(vxvv=solar_coords, radec=True, solarmotion='schoenrich') # should just be the sun's orbit
    sun.integrate(ts,mp,method='odeint')

    # get the orbit [R, vR, vT, z, vz, phi] (pos scaled by ro, vel scaled by vo)
    sun_data = sun.getOrbit()

    # plots the sun's motion with respect to Galactic Centre
    sunR = 8 * sun_data[:,0]
    sunphi = sun_data[:,5]
    sunX = sunR * np.cos(sunphi)
    sunY = sunR * np.sin(sunphi)
    sunZ = 8 * sun_data[:,3]
    plt.clf()
    plt.plot(sunX, sunY)
    plt.savefig('temp_plots/sunXY.png')

    plt.clf()
    plt.plot(sunX, sunZ)
    plt.savefig('temp_plots/sunXZ.png')

    # plot the XY of the sun's motion using galpy's plot function (w.r.t GC)
    plt.clf()
    sun.plot(d1='x', d2='y')
    plt.savefig('temp_plots/galpy_sunXY.png')

    sun.plot(d1='x', d2='z')
    plt.savefig('temp_plots/galpy_sunXZ.png')

    plt.clf()
    sun.plot(d1='R', d2='z')
    plt.savefig('temp_plots/galpy_sunRZ.png')

    #                                                        kpc, km/s
    # INITIALISING THE LSR (at XYZUVW (w.r.t sun) of [0,0,-0.025,0,220,0]
    R, vR, vT, z, vz, phi = 1., 0., 1., 0., 0., 0. # <--- Galpy units
    LSR_coords = [R, vR, vT, z, vz, phi]
    lsr = Orbit(vxvv=LSR_coords, solarmotion='schoenrich', vo=220, ro=8)
    lsr.integrate(ts, mp, method='odeint')

    # plots a perfect circle
    plt.clf()
    lsr.plot(d1='x', d2='y')
    plt.savefig('temp_plots/galpy_lsrXY.png')

    plt.clf()
    lsr.plot(d1='x', d2='z')
    plt.savefig('temp_plots/galpy_lsrXZ.png')

    # Manually reconstructing orbit
    lsr_data = lsr.getOrbit()
    lsrR = 8 * lsr_data[:,0]
    lsrphi = lsr_data[:,5]

    lsrX = lsrR * np.cos(lsrphi)
    lsrY = lsrR * np.sin(lsrphi)
    lsrZ = 8 * lsr_data[:,3]

    plt.clf()
    plt.plot(lsrX, lsrY)
    plt.savefig('temp_plots/lsrXY.png')
    plt.clf()
    plt.plot(lsrX, lsrZ)
    plt.savefig('temp_plots/lsrXZ.png')

    # plotting both sun and lsr
    plt.clf()
    plt.plot(lsrX, lsrY)
    plt.plot(sunX, sunY)
    plt.savefig('temp_plots/combXY.png')
    plt.clf()
    plt.plot(lsrX, lsrZ)
    plt.plot(sunX, sunZ)
    plt.savefig('temp_plots/combXZ.png')

    # Finding sun's path w.r.t the LSR in non-corotating frame
    relsunX = sunX - lsrX
    relsunY = sunY - lsrY
    relsunZ = sunZ - lsrZ

    plt.clf()
    plt.plot(relsunX, relsunY)
#.........这里部分代码省略.........

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import plot [as 别名]
#! /usr/bin/env python

from galpy.potential import MWPotential2014 as mp
from galpy.orbit import Orbit

import matplotlib.pyplot as plt
import numpy as np
import pdb

o = Orbit(vxvv=[1., 0.1, 1.1, 0., 0.1])
ts = np.linspace(0,100,10000)
o.integrate(ts, mp, method='odeint')

pdb.set_trace()
plot = o.plot()
o.plotE(normed=True)