本文整理汇总了Python中galpy.orbit.Orbit.integrate_dxdv方法的典型用法代码示例。如果您正苦于以下问题：Python Orbit.integrate_dxdv方法的具体用法？Python Orbit.integrate_dxdv怎么用？Python Orbit.integrate_dxdv使用的例子？那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类galpy.orbit.Orbit的用法示例。
在下文中一共展示了Orbit.integrate_dxdv方法的2个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: plot_liouville
# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate_dxdv [as 别名]
def plot_liouville(plotfilename):
E, Lz= -1.25, 0.6
aAS= actionAngleStaeckel(pot=MWPotential2014,c=True,delta=0.5)
#planarOrbit w/ the same energy
o= Orbit([0.8,numpy.sqrt(2.*(E-evalPot(0.8,0.,MWPotential2014)-(Lz/0.8)**2./2.)),Lz/0.8,0.])
#For the orbital period
fo= Orbit([0.8,numpy.sqrt(2.*(E-evalPot(0.8,0.,MWPotential2014)-(Lz/0.8)**2./2.)),Lz/0.8,0.,0.001,0.])
orbt= 2.*numpy.pi/aAS.actionsFreqs(fo)[4]
norb= 200.
nt= 20001
ts= numpy.linspace(0.,norb*orbt,nt)
start= time.time()
integrator= 'dopr54_c'
o.integrate_dxdv([1.,0.,0.,0.],ts,MWPotential2014,
method=integrator,
rectIn=True,rectOut=True)
dx= o.getOrbit_dxdv()[:,:]
o.integrate_dxdv([0.,1.,0.,0.],ts,MWPotential2014,
method=integrator,
rectIn=True,rectOut=True)
dy= o.getOrbit_dxdv()[:,:]
o.integrate_dxdv([0.,0.,1.,0.],ts,MWPotential2014,
method=integrator,
rectIn=True,rectOut=True)
dvx= o.getOrbit_dxdv()[:,:]
o.integrate_dxdv([0.,0.,0.,1.],ts,MWPotential2014,
method=integrator,
rectIn=True,rectOut=True)
dvy= o.getOrbit_dxdv()[:,:]
print integrator, time.time()-start
#Calculate Jacobian
jacs= numpy.array([numpy.linalg.det(numpy.array([dx[ii],dy[ii],
dvx[ii],dvy[ii]])) for ii in range(nt)])
breakt= 8.
pts= list(ts[ts < breakt])
pts.extend(list((ts[ts >= breakt])[::10]))
pts= numpy.array(pts)
pjacs= list(jacs[ts < breakt])
pjacs.extend(list((jacs[ts >= breakt])[::10]))
pjacs= numpy.array(pjacs)
print integrator, numpy.mean(jacs)-1.
bovy_plot.bovy_print(fig_width=3.25,fig_height=4.5)
pyplot.subplot(4,1,4)
bovy_plot.bovy_plot(pts/orbt,numpy.fabs(pjacs-1.),color='k',
loglog=True,gcf=True,
xrange=[0.5,norb],
yrange=[10.**-12.,10.**0.],
xlabel=r'$\mathrm{Number\ of\ orbital\ periods}$')
bovy_plot.bovy_text(r'$\texttt{dopr54\_c}$',
top_left=True,size=14.)
bovy_plot.bovy_text(0.1,10.**44.5,
r'$|\mathrm{Determinant\ of\ volume\ transformation}-1|$',
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.])
nullfmt = NullFormatter() # no labels
other_integrators= ['odeint',
'rk4_c','rk6_c']
for ii,integrator in enumerate(other_integrators):
start= time.time()
o.integrate_dxdv([1.,0.,0.,0.],ts,MWPotential2014,
method=integrator,
rectIn=True,rectOut=True)
dx= o.getOrbit_dxdv()[:,:]
o.integrate_dxdv([0.,1.,0.,0.],ts,MWPotential2014,
method=integrator,
rectIn=True,rectOut=True)
dy= o.getOrbit_dxdv()[:,:]
o.integrate_dxdv([0.,0.,1.,0.],ts,MWPotential2014,
method=integrator,
rectIn=True,rectOut=True)
dvx= o.getOrbit_dxdv()[:,:]
o.integrate_dxdv([0.,0.,0.,1.],ts,MWPotential2014,
method=integrator,
rectIn=True,rectOut=True)
dvy= o.getOrbit_dxdv()[:,:]
print integrator, time.time()-start
jacs= numpy.array([numpy.linalg.det(numpy.array([dx[jj],dy[jj],
dvx[jj],dvy[jj]])) for jj in range(nt)])
pts= list(ts[ts < breakt])
pts.extend(list((ts[ts >= breakt])[::10]))
pts= numpy.array(pts)
pjacs= list(jacs[ts < breakt])
pjacs.extend(list((jacs[ts >= breakt])[::10]))
pjacs= numpy.array(pjacs)
print integrator, numpy.mean(jacs)-1.
pyplot.subplot(4,1,ii+1)
if integrator == 'odeint':
yrange=[10.**-8.,10.**4.]
yticks= [10.**-8.,10.**-4.,1.,10.**4.]
else:
yrange=[10.**-12.,10.**0.]
yticks= [10.**-12.,10.**-8.,10.**-4.,1.]
bovy_plot.bovy_plot(pts/orbt,numpy.fabs(pjacs-1.),color='k',
loglog=True,gcf=True,
xrange=[0.5,norb],
yrange=yrange)
thisax= pyplot.gca()
#.........这里部分代码省略.........
示例2: MiyamotoNagaiPotential
# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate_dxdv [as 别名]
import sys
import numpy
from galpy.orbit import Orbit
from galpy.potential import MiyamotoNagaiPotential
if __name__ == '__main__':
# python orbitint4sigint.py symplec4_c full
mp= MiyamotoNagaiPotential(normalize=1.,a=0.5,b=0.05)
ts= numpy.linspace(0.,10000000.,1000001)
if sys.argv[2] == 'full' or sys.argv[2] == 'planar':
if sys.argv[2] == 'full':
o= Orbit([1.,0.1,1.1,0.1,0.1,0.])
elif sys.argv[2] == 'planar':
o= Orbit([1.,0.1,1.1,0.1])
o.integrate(ts,mp,method=sys.argv[1])
elif sys.argv[2] == 'planardxdv':
o= Orbit([1.,0.1,1.1,0.1])
o.integrate_dxdv([0.1,0.1,0.1,0.1],ts,mp,method=sys.argv[1])
sys.exit(0)