本文整理汇总了Python中galpy.orbit.Orbit.wr方法的典型用法代码示例。如果您正苦于以下问题：Python Orbit.wr方法的具体用法？Python Orbit.wr怎么用？Python Orbit.wr使用的例子？那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类galpy.orbit.Orbit的用法示例。
在下文中一共展示了Orbit.wr方法的1个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: calcj
# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import wr [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]),
#.........这里部分代码省略.........