Python Orbit.getOrbit方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import getOrbit [as 别名]
def observed_to_xyzuvw_orbit(obs, ts, lsr_orbit=None):
    """
    Convert six-parameter astrometric solution to XYZUVW orbit.

    Parameters
    ----------
    obs :   [RA (deg), DEC (deg), pi (mas),
             mu_ra (mas/yr), mu_dec (mas/yr), vlos (km/s)]
        Current kinematics
    ts : [ntimes] array
        times (in Myr) to traceback to

    lsr_orbit : Orbit
        the orbit of the local standard of rest for comparison, if None can
        calculate on the fly

    XYZUVW : [ntimes, 6] array
        The space position and velocities of the star in a co-rotating frame
        centred on the LSR
    """
    #convert times from Myr into bovy_time
    bovy_ts = ts / (bovy_conversion.time_in_Gyr(220.,8.)*1000) # bovy-time/Myr
    logging.info("max time in Myr: {}".format(np.max(ts)))
    logging.info("max time in Bovy time: {}".format(np.max(bovy_ts)))

    o = Orbit(vxvv=obs, radec=True, solarmotion='schoenrich')
    o.integrate(bovy_ts,mp,method='odeint')
    data = o.getOrbit()
    XYZUVW = galpy_coords_to_xyzuvw(data, bovy_ts)
    return XYZUVW

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import getOrbit [as 别名]
def trace_cartesian_orbit(xyzuvw_start, times=None, single_age=True,
                          potential=MWPotential2014, ro=8., vo=220.):
    """
    Given a star's XYZUVW relative to the LSR (at any time), project its
    orbit forward (or backward) to each of the times listed in *times*

    Positive times --> traceforward
    Negative times --> traceback

    TODO: Primary source of inefficiencies, 1366.2 (s)

    Parameters
    ----------
    xyzuvw : [pc,pc,pc,km/s,km/s,km/s]
    times : (float) or ([ntimes] float array)
        Myr - time of 0.0 must be present in the array. Times need not be
        spread linearly.
    single_age: (bool) {True}
        Set this flag if only providing a single age to trace to

    Returns
    -------
    xyzuvw_tf : [ntimes, 6] array
        [pc, pc, pc, km/s, km/s, km/s] - the traced orbit with positions
        and velocities

    Notes
    -----
    Profiling comments have been left in for future reference, but note
    that the profiling was done with previous versions of coordinate
    functions - ones that utilised astropy.units (and thus symbolic algebra)
    """
    if single_age:
        # replace 0 with some tiny number
        if times == 0.:
            times = 1e-15
        times = np.array([0., times])
    else:
        times = np.array(times)

    xyzuvw_start = np.copy(xyzuvw_start).astype(np.float)

    bovy_times = convert_myr2bovytime(times)

    # since the LSR is constant in chron coordinates, the starting point
    # is always treated as time 0
    galpy_coords = convert_cart2galpycoords(xyzuvw_start, ts=0.,
                                            ro=ro, vo=vo)
    o = Orbit(vxvv=galpy_coords, ro=ro, vo=vo)
    o.integrate(bovy_times, potential, method='odeint')

    xyzuvw = convert_galpycoords2cart(o.getOrbit(), bovy_times,
                                      ro=ro, vo=vo)
    if single_age:
        return xyzuvw[-1]
    return xyzuvw

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import getOrbit [as 别名]
def trace_orbit_xyzuvw(xyzuvw_then, times):
    """
    Given a star's XYZUVW relative to the LSR (at any time), project its
    orbit forward to each of the times listed in *times*

    Parameters
    ----------
    xyzuvw : [pc,pc,pc,km/s,km/s,km/s]
    times : [ntimes] array
        Myr

    Returns
    -------
    xyzuvw_tf : [ntimes, 6] array
        the traced-forward positions and velocities
    """
    XYZUVW_sun_now = np.array([0.,0.,0.025,11.1,12.24,7.25])

    # convert positions to kpc
    xyzuvw_then[:3] *= 1e-3

    # convert times from Myr to bovy units
    bovy_times = times*1e-3 / bovy_conversion.time_in_Gyr(220., 8.)
    logging.debug("Tracing up to {} Myr".format(times[-1]))
    logging.debug("Tracing up to {} Bovy yrs".format(bovy_times[-1]))

    # Galpy coordinates are from the vantage point of the sun.
    # So even though the sun wasn't where it is,  we still "observe" the star
    # from a solar height and motion
    XYZUVW_gp_tf = xyzuvw_then - XYZUVW_sun_now
    logging.debug("Galpy vector: {}".format(XYZUVW_gp_tf))

    l,b,dist = lbdist_from_XYZ(XYZUVW_gp_tf)
    vxvv = [l,b,dist,XYZUVW_gp_tf[3],XYZUVW_gp_tf[4],XYZUVW_gp_tf[5]]
    logging.debug("vxvv: {}".format(vxvv))
    otf = Orbit(vxvv=vxvv, lb=True, uvw=True, solarmotion='schoenrich')

    otf.integrate(bovy_times,mp,method='odeint')
    data_tf = otf.getOrbit()
    XYZUVW_tf = galpy_coords_to_xyzuvw(data_tf, bovy_times)
    logging.debug("Started orbit at {}".format(XYZUVW_tf[0]))
    logging.debug("Finished orbit at {}".format(XYZUVW_tf[-1]))
    return XYZUVW_tf

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import getOrbit [as 别名]
# V = lsr_orbit.V(0)
# W = lsr_orbit.W(0)
#
# results = np.array([U, V, W]).reshape(3)
# print('-- Inspecting the initialised velocities --')
# for i, vel in enumerate('UVW'):
#     print("For velocity {}:".format(vel))
#     print("    expected: {}".format(REFERENCE_SCHOENRICH_SOLAR_MOTION[i]))
#     print("    received: {}".format(results[i]))

# print("Since we initialised at the LSR we expect z to stay constant")
# print("z: {}".format([lsr_orbit.z(t) for t in np.linspace(0,age,10)]))
# print("Since z is constant, we expect W to be constant")
# print("z: {}".format([lsr_orbit.W(t) for t in np.linspace(0,age,10)]))

orbit = lsr_orbit.getOrbit()
back_orbit = back_lsr_orbit.getOrbit()

dim_labels = ['R','vR','vT','z','vz','phi']
print('____ PRINTING ORBIT, DIMENSION BY DIMENSION ____')
print('_____ forward _____')
for i in range(6):
    print('----- {}:{:5} -----'.format(i, dim_labels[i]))
    print(orbit[:,i])
print('_____ backward _____')
for i in range(6):
    print('----- {}:{:5} -----'.format(i, dim_labels[i]))
    print(back_orbit[:,i])

import sys
sys.path.insert(0, '..')

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import getOrbit [as 别名]
def test_careful_traceback_and_forward():
    """Step by step, project orbit forward, then backward"""
    bovy_times = np.array([0., np.pi/3.])
    chron_times = torb.convert_bovytime2myr(bovy_times)

    init_pos_chron = np.array([
        4000, 8000.*np.sqrt(3)/2, 0,
        np.sin(np.pi/3)*220.,
        -np.cos(np.pi/3)*220.,
        0
    ])
    init_pos_galpy = torb.convert_cart2galpycoords(init_pos_chron, ts=0.)

    assert np.allclose(np.array([1.,0,1,0,0,np.pi/3.]),
                       init_pos_galpy)

    o = Orbit(vxvv=init_pos_galpy, ro=8., vo=220.)
    o.integrate(bovy_times, MWPotential2014, method='odeint')

    orbit_galpy = o.getOrbit()
    assert np.allclose(init_pos_galpy, orbit_galpy[0])
    assert np.allclose(init_pos_galpy
                       + np.array([0.,0.,0.,0.,0.,bovy_times[-1]]),
                       orbit_galpy[-1])

    orbit_chron = torb.convert_galpycoords2cart(orbit_galpy,
                                                ts=bovy_times)
    assert np.allclose(init_pos_chron, orbit_chron[0])
    assert np.allclose(init_pos_chron,
                       orbit_chron[-1])

    # Setup for backwards time integration
    # Currently at time of PI/3
    back_init_pos_chron = orbit_chron[-1]
    back_init_pos_galpy = torb.convert_cart2galpycoords(
        back_init_pos_chron,
        bovy_times=bovy_times[-1],
    )

    assert np.allclose(back_init_pos_galpy,
                       torb.convert_cart2galpycoords(
                           back_init_pos_chron,
                           bovy_times=bovy_times[-1]
                       ))

    back_o = Orbit(vxvv=back_init_pos_galpy, ro=8., vo=220.)
    back_o.integrate(-1*bovy_times, MWPotential2014, method='odeint')

    back_orbit_galpy = back_o.getOrbit()

    assert np.allclose(back_init_pos_galpy, back_orbit_galpy[0])
    assert np.allclose(back_init_pos_galpy
                       - np.array([0.,0.,0.,0.,0.,bovy_times[-1]]),
                       back_orbit_galpy[-1])
    assert np.allclose(init_pos_galpy, back_orbit_galpy[-1])

    back_orbit_chron = torb.convert_galpycoords2cart(
        back_orbit_galpy,
        ts=bovy_times[::-1],
    )

    assert np.allclose(init_pos_chron, back_orbit_chron[-1])

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import getOrbit [as 别名]
 def _parse_args(self,freqsAngles=True,_firstFlip=False,*args):
     """Helper function to parse the arguments to the __call__ and actionsFreqsAngles functions"""
     from galpy.orbit import Orbit
     RasOrbit= False
     integrated= True #whether the orbit was already integrated when given
     if len(args) == 5 or len(args) == 3:
         raise IOError("Must specify phi for actionAngleIsochroneApprox")
     if len(args) == 6 or len(args) == 4:
         if len(args) == 6:
             R,vR,vT, z, vz, phi= args
         else:
             R,vR,vT, phi= args
             z, vz= 0., 0.
         if isinstance(R,float):
             o= Orbit([R,vR,vT,z,vz,phi])
             o.integrate(self._tsJ,pot=self._pot,method=self._integrate_method)
             this_orbit= o.getOrbit()
             R= nu.reshape(this_orbit[:,0],(1,self._ntintJ))
             vR= nu.reshape(this_orbit[:,1],(1,self._ntintJ))
             vT= nu.reshape(this_orbit[:,2],(1,self._ntintJ))
             z= nu.reshape(this_orbit[:,3],(1,self._ntintJ))
             vz= nu.reshape(this_orbit[:,4],(1,self._ntintJ))           
             phi= nu.reshape(this_orbit[:,5],(1,self._ntintJ))           
             integrated= False
         if len(R.shape) == 1: #not integrated yet
             os= [Orbit([R[ii],vR[ii],vT[ii],z[ii],vz[ii],phi[ii]]) for ii in range(R.shape[0])]
             RasOrbit= True
             integrated= False
     if isinstance(args[0],Orbit) \
             or (isinstance(args[0],list) and isinstance(args[0][0],Orbit)) \
             or RasOrbit:
         if RasOrbit:
             pass
         elif not isinstance(args[0],list):
             os= [args[0]]
             if len(os[0]._orb.vxvv) == 3 or len(os[0]._orb.vxvv) == 5:
                 raise IOError("Must specify phi for actionAngleIsochroneApprox")
         else:
             os= args[0]
             if len(os[0]._orb.vxvv) == 3 or len(os[0]._orb.vxvv) == 5:
                 raise IOError("Must specify phi for actionAngleIsochroneApprox")
         if not hasattr(os[0],'orbit'): #not integrated yet
             if _firstFlip:
                 for o in os:
                     o._orb.vxvv[1]= -o._orb.vxvv[1]
                     o._orb.vxvv[2]= -o._orb.vxvv[2]
                     o._orb.vxvv[4]= -o._orb.vxvv[4]
             [o.integrate(self._tsJ,pot=self._pot,
                          method=self._integrate_method) for o in os]
             if _firstFlip:
                 for o in os:
                     o._orb.vxvv[1]= -o._orb.vxvv[1]
                     o._orb.vxvv[2]= -o._orb.vxvv[2]
                     o._orb.vxvv[4]= -o._orb.vxvv[4]
                     o._orb.orbit[:,1]= -o._orb.orbit[:,1]
                     o._orb.orbit[:,2]= -o._orb.orbit[:,2]
                     o._orb.orbit[:,4]= -o._orb.orbit[:,4]
             integrated= False
         ntJ= os[0].getOrbit().shape[0]
         no= len(os)
         R= nu.empty((no,ntJ))
         vR= nu.empty((no,ntJ))
         vT= nu.empty((no,ntJ))
         z= nu.zeros((no,ntJ))+10.**-7. #To avoid numpy warnings for
         vz= nu.zeros((no,ntJ))+10.**-7. #planarOrbits
         phi= nu.empty((no,ntJ))
         for ii in range(len(os)):
             this_orbit= os[ii].getOrbit()
             R[ii,:]= this_orbit[:,0]
             vR[ii,:]= this_orbit[:,1]
             vT[ii,:]= this_orbit[:,2]
             if this_orbit.shape[1] == 6:
                 z[ii,:]= this_orbit[:,3]
                 vz[ii,:]= this_orbit[:,4]
                 phi[ii,:]= this_orbit[:,5]
             else:
                 phi[ii,:]= this_orbit[:,3]
     if freqsAngles and not integrated: #also integrate backwards in time, such that the requested point is not at the edge
         no= R.shape[0]
         nt= R.shape[1]
         oR= nu.empty((no,2*nt-1))
         ovR= nu.empty((no,2*nt-1))
         ovT= nu.empty((no,2*nt-1))
         oz= nu.zeros((no,2*nt-1))+10.**-7. #To avoid numpy warnings for
         ovz= nu.zeros((no,2*nt-1))+10.**-7. #planarOrbits
         ophi= nu.empty((no,2*nt-1))
         if _firstFlip:
             oR[:,:nt]= R[:,::-1]
             ovR[:,:nt]= vR[:,::-1]
             ovT[:,:nt]= vT[:,::-1]
             oz[:,:nt]= z[:,::-1]
             ovz[:,:nt]= vz[:,::-1]
             ophi[:,:nt]= phi[:,::-1]
         else:
             oR[:,nt-1:]= R
             ovR[:,nt-1:]= vR
             ovT[:,nt-1:]= vT
             oz[:,nt-1:]= z
             ovz[:,nt-1:]= vz
             ophi[:,nt-1:]= phi
#.........这里部分代码省略.........

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import getOrbit [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 getOrbit [as 别名]
def lsr_nonsense():
    """
    Mucking around, showing how to get the LSR orbit
    """

    RO = 8.
    VO = 220.
    BOVY_TIME_CONVERSION = bovy_conversion.time_in_Gyr(VO, RO) * 1000 # Myr/bovy_time

    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("Myr")

    max_age = orbit_time.value / BOVY_TIME_CONVERSION
    ntimes = 100
    ts = np.linspace(0, max_age, ntimes)

    # demo a star (with vT=220, vR=0, vZ=0, z=0, phi=0.1 pi) staying
    # fixed in our coordinate frame
    R, vR, vT, z, vz, phi = 1., 0., 1., 0., 0., 0.
    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')

    lsr_data = lsr.getOrbit()
    lsrR = RO * lsr_data[:,0]
    lsrphi = lsr_data[:,5]

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

    R, vR, vT, z, vz, phi = 1., 0., 1., 0., 0., 0.25*np.pi
    rot_lsr_coords = [R, vR, vT, z, vz, phi]
    rot_lsr = Orbit(vxvv=rot_lsr_coords, solarmotion='schoenrich', vo=220, ro=8)
    rot_lsr.integrate(ts, mp, method='odeint')

    rot_lsr_data = rot_lsr.getOrbit()

    # putting into corotating cartesian system centred on LSR
    XYZUVW_rot = galpy_coords_to_xyzuvw(rot_lsr_data, ts)
    plt.clf()
    plt.plot(XYZUVW_rot[:,0], XYZUVW_rot[:,1])
    plt.savefig("temp_plots/rotXY.png")


    orbit_time = (perimeter / velocity).to("Myr")
    ts = np.linspace(0., 10*orbit_time.value, 1000) / BOVY_TIME_CONVERSION
    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()
    XYZUVW_sun = galpy_coords_to_xyzuvw(sun_data, ts)
    plt.clf()
    plt.plot(XYZUVW_sun[:,0], XYZUVW_sun[:,1])
    plt.savefig("temp_plots/sunXY.png")
    plt.clf()
    plt.plot(XYZUVW_sun[:,0], XYZUVW_sun[:,2])
    plt.savefig("temp_plots/sunXZ.png")