Python Orbit.integrate方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [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 integrate [as 别名]
def test_correctInitialSolarMotion():
    age = np.pi # Half a galactic revolution
    
    # For reference: the expected solar motion as provided by Schoenrich
    REFERENCE_SCHOENRICH_SOLAR_MOTION = np.array([-11.1, -12.24, -7.25])
    ntimesteps = 10

    # vxvv is in cylindrical coordinates here. [R,vR,vT(,z,vz,phi)]
    lsr_orbit = Orbit(vxvv=[1,0,1,0,0,0],vo=220,ro=8,solarmotion='schoenrich')
    lsr_orbit.integrate(np.linspace(0.,age,ntimesteps), MWPotential2014)

    U = lsr_orbit.U(0)
    V = lsr_orbit.V(0)
    W = lsr_orbit.W(0)

    results = np.array([U, V, W]).reshape(3)
    for i, vel in enumerate('UVW'):
        print("For velocity {}:".format(vel))
        print("    expected: {}".format(REFERENCE_SCHOENRICH_SOLAR_MOTION[i]))
        print("    received: {}".format(results[i]))

    assert (np.allclose(REFERENCE_SCHOENRICH_SOLAR_MOTION, results)),\
        '!!! Using galpy version {} Need galpy version 1.1 !!!'.format(
            galpy.__version__
        )

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [as 别名]
def traceback2(params,times):
    """Trace forward a cluster. First column of returned array is the
    position of the cluster at a given age.

    Parameters
    ----------
    times: float array
        Times to trace forward, in Myr. Note that positive numbers are going
        forward in time.
    params: float array
        [RA,DE,Plx,PM(RA),PM(DE),RV]
        RA = Right Ascension (Deg)
        DE = Declination (Deg)
        Plx = Paralax (Mas)
        PM(RA) = Proper motion (Right Ascension) (mas/yr)
        PM(DE) = Proper motion (Declination) (mas/yr)
        RV = Radial Velocity (km/s)
    age: Age of cluster, in Myr
    
    """
    #FIXME: This is very out of date and should be deleted!!!
    #Times in Myr
    ts = -(times/1e3)/bovy_conversion.time_in_Gyr(220.,8.)
    nts = len(times)

    #Positions and velocities in the co-rotating solar reference frame.
    xyzuvw = np.zeros( (1,nts,6) )

    #Trace forward the local standard of rest.
    lsr_orbit= Orbit(vxvv=[1.,0,1,0,0.,0],vo=220,ro=8)
    lsr_orbit.integrate(ts,MWPotential2014)#,method='odeint')

    xyzuvw = integrate_xyzuvw(params,ts,lsr_orbit,MWPotential2014)
    return xyzuvw

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [as 别名]
def test_orbitIntegrationC():

    #_____initialize some KKSPot_____
    Delta = 1.0
    pot = KuzminKutuzovStaeckelPotential(ac=20.,Delta=Delta,normalize=True)

    #_____initialize an orbit (twice)_____
    vxvv = [1.,0.1,1.1,0.,0.1]
    o_P= Orbit(vxvv=vxvv)
    o_C= Orbit(vxvv=vxvv)

    #_____integrate the orbit with python and C_____
    ts= numpy.linspace(0,100,101)
    o_P.integrate(ts,pot,method='leapfrog')  #python
    o_C.integrate(ts,pot,method='leapfrog_c')#C

    for ii in range(5):
        exp3= -1.7
        if   ii == 0: Python, CC, string, exp1, exp2 = o_P.R(ts) , o_C.R(ts) , 'R' , -5., -10.
        elif ii == 1: Python, CC, string, exp1, exp2 = o_P.z(ts) , o_C.z(ts) , 'z' , -3.25, -4.
        elif ii == 2: Python, CC, string, exp1, exp2 = o_P.vR(ts), o_C.vR(ts), 'vR', -3., -10.
        elif ii == 3: Python, CC, string, exp1, exp2, exp3 = o_P.vz(ts), o_C.vz(ts), 'vz', -3., -4., -1.3
        elif ii == 4: Python, CC, string, exp1, exp2 = o_P.vT(ts), o_C.vT(ts), 'vT', -5., -10.

        rel_diff = numpy.fabs((Python-CC)/CC) < 10.**exp1
        abs_diff = (numpy.fabs(Python-CC) < 10.**exp2) * (numpy.fabs(Python) < 10.**exp3)
        assert numpy.all(rel_diff+abs_diff), \
            'Orbit integration for '+string+' coordinate different in ' + \
            'C and Python implementation.'

    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [as 别名]
def test_estimateDelta():

    #_____initialize some KKSPot_____
    Delta = 1.0
    pot = KuzminKutuzovStaeckelPotential(ac=20.,Delta=Delta,normalize=True)

    #_____initialize an orbit (twice)_____
    vxvv = [1.,0.1,1.1,0.01,0.1]
    o= Orbit(vxvv=vxvv)

    #_____integrate the orbit with C_____
    ts= numpy.linspace(0,101,100)
    o.integrate(ts,pot,method='leapfrog_c')


    #____estimate Focal length Delta_____
    #for each time step individually:
    deltas_estimate = numpy.zeros(len(ts))
    for ii in range(len(ts)):
        deltas_estimate[ii] = estimateDeltaStaeckel(pot,o.R(ts[ii]),o.z(ts[ii]))

    assert numpy.all(numpy.fabs(deltas_estimate - Delta) < 10.**-8), \
            'Focal length Delta estimated along the orbit is not constant.'

    #for all time steps together:
    delta_estimate = estimateDeltaStaeckel(pot,o.R(ts),o.z(ts))
    
    assert numpy.fabs(delta_estimate - Delta) < 10.**-8, \
            'Focal length Delta estimated from the orbit is not the same as the input focal length.'

    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [as 别名]
def test_actionConservation():

    #_____initialize some KKSPot_____
    Delta = 1.0
    pot = KuzminKutuzovStaeckelPotential(ac=20.,Delta=Delta,normalize=True)

    #_____initialize an orbit (twice)_____
    vxvv = [1.,0.1,1.1,0.01,0.1]
    o= Orbit(vxvv=vxvv)

    #_____integrate the orbit with C_____
    ts= numpy.linspace(0,101,100)
    o.integrate(ts,pot,method='leapfrog_c')

    #_____Setup ActionAngle object and calculate actions (Staeckel approximation)_____
    aAS = actionAngleStaeckel(pot=pot,delta=Delta,c=True)
    jrs,lzs,jzs = aAS(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts))

    assert numpy.all(numpy.fabs(jrs - jrs[0]) < 10.**-8.), \
        'Radial action is not conserved along orbit.'

    assert numpy.all(numpy.fabs(lzs - lzs[0]) < 10.**-8.), \
        'Angular momentum is not conserved along orbit.'

    assert numpy.all(numpy.fabs(jzs - jzs[0]) < 10.**-8.), \
        'Vertical action is not conserved along orbit.'

    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [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 integrate [as 别名]
def evolveorbit(icon, ti, tau, pot):
    global x
    global y
    o = Orbit(vxvv=icon) # [R,vR,vT,z,vz,phi]
    tf = ti+tau
    ts = np.linspace(ti,tf,100)
    o.integrate(ts, pot, method = 'leapfrog')
    x.append(o.x(ts[0]))
    y.append(o.y(ts[0]))
    return [o.R(tf),o.vR(tf),o.vT(tf),o.z(tf),o.vz(tf),o.phi(tf)]

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [as 别名]
def plot_aagrid(plotfilename1,plotfilename2):
    #Setup orbit
    E, Lz= -1.25, 0.6
    o= 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.])
    delta= 0.434
    #Integrate the orbit, setup Staeckel already to calculate the period
    aAS= actionAngleStaeckel(pot=MWPotential2014,delta=delta,c=True)
    orbt= 2.*numpy.pi/aAS.actionsFreqs(o)[4]  
    norb= 5.
    nt= 501
    ts= numpy.linspace(0.,norb*orbt,nt)
    o.integrate(ts,MWPotential2014,method='symplec4_c')
    #First do adiabatic
    aAA= actionAngleAdiabatic(pot=MWPotential2014,gamma=1.,c=True)
    aAAG= actionAngleAdiabaticGrid(pot=MWPotential2014,gamma=1.,c=True,
                                   nR=31,nEz=31,nEr=51,nLz=51)
    jfa= aAA(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
    jfag= aAAG(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
    #First do adiabatic
    #aAS already setup
    aASG= actionAngleStaeckelGrid(pot=MWPotential2014,delta=delta,c=True,
                                  nE=51,npsi=51,nLz=51)
    jfs= aAS(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
    jfsg= aASG(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
    bovy_plot.bovy_print()
    line1= bovy_plot.bovy_plot(jfa[0],jfa[2],'r.',
                               xrange=[0.045,0.055],
                               yrange=[0.0075,0.011],
                               xlabel=r'$J_R$',ylabel=r'$J_z$',zorder=2)
    line2= bovy_plot.bovy_plot(jfag[0],jfag[2],'rx',overplot=True,zorder=1)
    bovy_plot.bovy_plot(jfs[0],jfs[2],'k,',overplot=True)
    pyplot.legend((line1[0],line2[0]),
                  (r'$\mathrm{\texttt{actionAngleAdiabatic}}$',
                   r'$\mathrm{\texttt{actionAngleAdiabaticGrid}}$',),
                  loc='upper right',#bbox_to_anchor=(.91,.375),
                  numpoints=1,
                  prop={'size':14},
                  frameon=False)
    bovy_plot.bovy_end_print(plotfilename1)
    #Zoom of Staeckel
    line1= bovy_plot.bovy_plot(jfs[0],jfs[2],'k.',
                               xrange=[0.05025,0.05145],
                               yrange=[0.0086,0.00933],
                               xlabel=r'$J_R$',ylabel=r'$J_z$')
    line2= bovy_plot.bovy_plot(jfsg[0],jfsg[2],'kx',overplot=True)
    pyplot.legend((line1[0],line2[0]),
                  (r'$\mathrm{\texttt{actionAngleStaeckel}}$',
                   r'$\mathrm{\texttt{actionAngleStaeckelGrid}}$',),
                  loc='upper right',#bbox_to_anchor=(.91,.375),
                  numpoints=1,
                  prop={'size':14},
                  frameon=False)
    bovy_plot.bovy_end_print(plotfilename2)
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [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 integrate [as 别名]
def integrate_xyzuvw(params,ts,lsr_orbit,MWPotential2014):
    """Convenience function"""
    vxvv = params.copy()
    vxvv[2]=1.0/params[2]
    o = Orbit(vxvv=vxvv, radec=True, solarmotion='schoenrich')
    o.integrate(ts,MWPotential2014,method='odeint')
    xyzuvw = np.zeros( (len(ts),6) )
    xyzuvw[:,0] = 1e3*(o.x(ts)-lsr_orbit.x(ts))
    xyzuvw[:,1] = 1e3*(o.y(ts)-lsr_orbit.y(ts))
    xyzuvw[:,2] = 1e3*(o.z(ts))
    xyzuvw[:,3] = o.U(ts) - lsr_orbit.U(ts)
    xyzuvw[:,4] = o.V(ts) - lsr_orbit.V(ts)
    xyzuvw[:,5] = o.W(ts)
    return xyzuvw

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [as 别名]
def test_surfacesection():
    #Preliminary code
    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)
    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)
    def surface_section(Rs,zs,vRs):
        # Find points where the orbit crosses z from - to +
        shiftzs= numpy.roll(zs,-1)
        indx= (zs[:-1] < 0.)*(shiftzs[:-1] > 0.)
        return (Rs[:-1][indx],vRs[:-1][indx])
    # Calculate and plot the surface of section
    ts= numpy.linspace(0.,1000.,20001) # long integration
    o1.integrate(ts,MWPotential2014)
    o2.integrate(ts,MWPotential2014)
    sect1Rs,sect1vRs=surface_section(o1.R(ts),o1.z(ts),o1.vR(ts))
    sect2Rs,sect2vRs=surface_section(o2.R(ts),o2.z(ts),o2.vR(ts))
    from matplotlib.pyplot import plot, xlim, ylim
    plot(sect1Rs,sect1vRs,'bo',mec='none')
    xlim(0.3,1.); ylim(-0.69,0.69)
    plot(sect2Rs,sect2vRs,'yo',mec='none')
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [as 别名]
def get_lsr_orbit(times, Potential=MWPotential2014):
    """Integrate the local standard of rest backwards in time, in order to provide
    a reference point to our XYZUVW coordinates.
    
    Parameters
    ----------
    times: numpy float array
        Times at which we want the orbit of the local standard of rest.
        
    Potential: galpy potential object.
    """
    ts = -(times/1e3)/bovy_conversion.time_in_Gyr(220.,8.)
    lsr_orbit= Orbit(vxvv=[1.,0,1,0,0.,0],vo=220,ro=8, solarmotion='schoenrich')
    lsr_orbit.integrate(ts,Potential)#,method='odeint')
    return lsr_orbit

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import integrate [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 integrate [as 别名]
def integrate_xyzuvw(params,times,lsr_orbit=None,Potential=MWPotential2014):
    """Convenience function. Integrates the motion of a star backwards in time.
    
    Parameters
    ----------
    params: numpy array 
        Kinematic parameters (RAdeg,DEdeg,Plx,pmRA,pmDE,RV)
        
    ts: times for back propagation, in Myr.
        
    lsr_orbit:
        WARNING: messy...
        lsr_orbit= Orbit(vxvv=[1.,0,1,0,0.,0],vo=220,ro=8, solarmotion='schoenrich')
        lsr_orbit.integrate(ts,MWPotential2014,method='odeint')
    MWPotential2014:
        WARNING: messy...
        from galpy.potential import MWPotential2014        
    """
    #We allow MWPotential and lsr_orbit to be passed for speed, but can compute/import 
    #now.
    if not lsr_orbit:
        lsr_orbit = get_lsr_orbit(times)
        
    #Convert times to galpy units:
    ts = -(times/1e3)/bovy_conversion.time_in_Gyr(220.,8.)
        
    params = np.array(params)
    vxvv = params.copy()
    #We'd prefer to pass parallax in mas, but distance in kpc is accepted. This
    #reciporical could be done elsewhere...
    vxvv[2]=1.0/params[2]   
    o = Orbit(vxvv=vxvv, radec=True, solarmotion='schoenrich')
    
    o.integrate(ts,Potential)#,method='odeint')
    xyzuvw = np.zeros( (len(ts),6) )
    xyzuvw[:,0] = 1e3*(o.x(ts)-lsr_orbit.x(ts))
    xyzuvw[:,1] = 1e3*(o.y(ts)-lsr_orbit.y(ts))
    #The lsr_orbit is zero in the z direction by definition - the local standard of 
    #rest is at the midplane of the Galaxy
    xyzuvw[:,2] = 1e3*(o.z(ts))  
    #UVW is relative to the sun. We *could* have used vx, vy and vz. Would these
    #have been relative to the LSR?
    xyzuvw[:,3] = o.U(ts) - lsr_orbit.U(ts)
    xyzuvw[:,4] = o.V(ts) - lsr_orbit.V(ts)
    xyzuvw[:,5] = o.W(ts) - lsr_orbit.W(ts) #NB This line changed !!!
    return xyzuvw