Python Orbit.vT方法代码示例

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vT [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 vT [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 vT [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 vT [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 vT [as 别名]
def test_actionAngleTorus_orbit():
    from galpy.actionAngle import actionAngleTorus
    from galpy.potential import MWPotential2014
    from galpy.orbit import Orbit
    # Set up instance
    aAT= actionAngleTorus(pot=MWPotential2014,tol=10.**-5.)
    jr,jphi,jz= 0.05,1.1,0.025
    # First calculate frequencies and the initial RvR
    RvRom= aAT.xvFreqs(jr,jphi,jz,
                       numpy.array([0.]),
                       numpy.array([1.]),
                       numpy.array([2.]))
    om= RvRom[1:]
    # Angles along an orbit
    ts= numpy.linspace(0.,100.,1001)
    angler= ts*om[0]
    anglephi= 1.+ts*om[1]
    anglez= 2.+ts*om[2]
    # Calculate the orbit using actionAngleTorus
    RvR= aAT(jr,jphi,jz,angler,anglephi,anglez).T
    # Calculate the orbit using orbit integration
    orb= Orbit([RvRom[0][:,0],RvRom[0][:,1],RvRom[0][:,2],
                RvRom[0][:,3],RvRom[0][:,4],RvRom[0][:,5]])
    orb.integrate(ts,MWPotential2014)
    # Compare
    tol= -3.
    assert numpy.all(numpy.fabs(orb.R(ts)-RvR[0]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in R'
    assert numpy.all(numpy.fabs(orb.vR(ts)-RvR[1]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in vR'
    assert numpy.all(numpy.fabs(orb.vT(ts)-RvR[2]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in vT'
    assert numpy.all(numpy.fabs(orb.z(ts)-RvR[3]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in z'
    assert numpy.all(numpy.fabs(orb.vz(ts)-RvR[4]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in vz'
    assert numpy.all(numpy.fabs(orb.phi(ts)-RvR[5]) < 10.**tol), \
        'Integrated orbit does not agree with torus orbit in phi'
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vT [as 别名]
def calc_progenitor_actions(savefilename):
    # Setup potential
    lp = potential.LogarithmicHaloPotential(normalize=1.0, q=0.9)
    # Setup orbit
    x, z, y, vx, vz, vy = -11.63337239, -10.631736273934635, -20.76235661, -128.8281653, 79.172383882274971, 42.88727925
    R, phi, z = bovy_coords.rect_to_cyl(x, y, z)
    vR, vT, vz = bovy_coords.rect_to_cyl_vec(vx, vy, vz, R, phi, z, cyl=True)
    R /= 8.0
    z /= 8.0
    vR /= 220.0
    vT /= 220.0
    vz /= 220.0
    o = Orbit([R, vR, vT, z, vz, phi])
    ts = numpy.linspace(0.0, 5.125 * 220.0 / 8.0, 1313)  # times of the snapshots
    o.integrate(ts, lp, method="dopr54_c")
    if "aas" in savefilename:
        aA = actionAngleStaeckel(pot=lp, delta=1.20, c=True)
    else:
        aA = actionAngleIsochroneApprox(pot=lp, b=0.8)
    # Now calculate actions, frequencies, and angles for all positions
    Rs = o.R(ts)
    vRs = o.vR(ts)
    vTs = o.vT(ts)
    zs = o.z(ts)
    vzs = o.vz(ts)
    phis = o.phi(ts)
    csvfile = open(savefilename, "wb")
    writer = csv.writer(csvfile, delimiter=",")
    for ii in range(len(ts)):
        acfs = aA.actionsFreqsAngles(Rs[ii], vRs[ii], vTs[ii], zs[ii], vzs[ii], phis[ii])
        writer.writerow(
            [acfs[0][0], acfs[1][0], acfs[2][0], acfs[3][0], acfs[4][0], acfs[5][0], acfs[6][0], acfs[7][0], acfs[8][0]]
        )
        csvfile.flush()
    csvfile.close()
    return None

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vT [as 别名]
def plot_aaspher_conservation(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.])
    #Integrate the orbit to estimate an equivalent b
    nt= 1001
    ts= numpy.linspace(0.,20.,nt)
    o.integrate(ts,MWPotential2014,method='symplec4_c')
    b= estimateBIsochrone(o.R(ts),o.z(ts),pot=MWPotential2014)
    print b
    b= 0.3
    #Now integrate the orbit in the isochronePotential
    ip= IsochronePotential(normalize=1.,b=b)
    aAI= actionAngleIsochrone(ip=ip)
    orbt= 2.*numpy.pi/aAI.actionsFreqs(o)[4]    
    norb= 200.
    nt= 20001
    ts= numpy.linspace(0.,norb*orbt,nt)
    o.integrate(ts,ip,method='symplec4_c')
    #Calculate actions, frequencies, and angles
    jfa= aAI.actionsFreqsAngles(o.R(ts),o.vR(ts),o.vT(ts),
                                o.z(ts),o.vz(ts),o.phi(ts))
    dJs= numpy.fabs((jfa[0]-numpy.mean(jfa[0]))/numpy.mean(jfa[0]))
    dOrs= numpy.fabs((jfa[3]-numpy.mean(jfa[3]))/numpy.mean(jfa[3]))
    dOzs= numpy.fabs((jfa[5]-numpy.mean(jfa[5]))/numpy.mean(jfa[5]))
    print "frequencies", numpy.mean(dOrs), numpy.mean(dOzs)
    ar= dePeriod(numpy.reshape(jfa[6],(1,len(ts)))).flatten()
    az= dePeriod(numpy.reshape(jfa[8],(1,len(ts)))).flatten()
    danglers= numpy.fabs(ar-numpy.mean(jfa[3])*ts-jfa[6][0])/2./numpy.pi
    danglezs= numpy.fabs(az-numpy.mean(jfa[5])*ts-jfa[8][0])/2./numpy.pi
    #Break up
    breakt= 50.
    pts= parse_break(ts,ts < breakt)
    pdJs= parse_break(dJs,ts < breakt)
    pdanglers= parse_break(danglers,ts < breakt)
    pdanglezs= parse_break(danglezs,ts < breakt)
    #dAngles
    bovy_plot.bovy_print()
    pyplot.subplot(2,1,1)
    bovy_plot.bovy_plot(pts/orbt,
                        pdJs,
                        color='k',loglog=True,gcf=True,
                        xrange=[0.5,norb],
                        yrange=[10.**-12.,1.])
    bovy_plot.bovy_text(r'$\texttt{actionAngleIsochrone}$',
                        top_left=True,size=14.)
    ax= pyplot.gca()
    ax.yaxis.set_ticks([10.**-12.,10.**-8.,10.**-4.,1.])
    nullfmt   = NullFormatter()         # no labels
    ax.xaxis.set_major_formatter(nullfmt)
    #Same for actionAngleSpherical
    aAS= actionAngleSpherical(pot=ip)
    tts= ts[::1]
    jfa= aAS.actionsFreqsAngles(o.R(tts),o.vR(tts),o.vT(tts),
                                o.z(tts),o.vz(tts),o.phi(tts),
                                fixed_quad=True)
    #dJr
    dJs= numpy.fabs((jfa[0]-numpy.mean(jfa[0]))/numpy.mean(jfa[0]))
    dOrs= numpy.fabs((jfa[3]-numpy.mean(jfa[3]))/numpy.mean(jfa[3]))
    dOzs= numpy.fabs((jfa[5]-numpy.mean(jfa[5]))/numpy.mean(jfa[5]))
    print "frequencies", numpy.mean(dOrs), numpy.mean(dOzs)
    #dAngles
    ar= dePeriod(numpy.reshape(jfa[6],(1,len(tts)))).flatten()
    az= dePeriod(numpy.reshape(jfa[8],(1,len(tts)))).flatten()
    danglers= numpy.fabs(ar-numpy.mean(jfa[3])*tts-jfa[6][0])/2./numpy.pi
    danglezs= numpy.fabs(az-numpy.mean(jfa[5])*tts-jfa[8][0])/2./numpy.pi
    print numpy.mean(danglers)
    print numpy.mean(danglezs)
    ptts= parse_break(tts,tts < breakt)
    pdJs= parse_break(dJs,tts < breakt)
    pyplot.subplot(2,1,2)
    bovy_plot.bovy_plot(ptts/orbt,
                        pdJs,
                        color='k',loglog=True,gcf=True,
                        xrange=[0.5,norb],
                        yrange=[10.**-12.,1.],
                        xlabel=r'$\mathrm{Number\ of\ orbital\ periods}$')
    bovy_plot.bovy_text(r'$\texttt{actionAngleSpherical}$',
                        top_left=True,size=14.)
    bovy_plot.bovy_text(0.175,10.**2.,r'$\left|\Delta J_R / J_R\right|$',
                        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.])
    bovy_plot.bovy_end_print(plotfilename1)    
    #Now plot the deviations in the angles
    bovy_plot.bovy_print()
    pyplot.subplot(2,1,1)
    liner= bovy_plot.bovy_plot(pts/orbt,
                               pdanglers,
                               color='k',ls='-',loglog=True,gcf=True,
                               xrange=[0.5,norb],
                               yrange=[10.**-12.,1.])
    linez= bovy_plot.bovy_plot(pts/orbt,
                               pdanglezs,
                               color='k',ls='--',overplot=True)
    legend1= pyplot.legend((liner[0],linez[0]),
                           (r'$\theta_R$',
                            r'$\theta_z$'),
#.........这里部分代码省略.........

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vT [as 别名]
def impulse_deltav_general_fullplummerintegration(v,x,b,w,x0,v0,galpot,GM,rs,
                                                  tmaxfac=10.,N=1000,
                                                integrate_method='symplec4_c'):
    """
    NAME:

       impulse_deltav_general_fullplummerintegration

    PURPOSE:

       calculate the delta velocity to due an encounter with a moving Plummer sphere and galactic potential relative to just in galactic potential

    INPUT:

       v - velocity of the stream (nstar,3)

       x - position along the stream (nstar,3)

       b - impact parameter

       w - velocity of the subhalo (3)

       x0 - position of closest approach (3)

       v0 - velocity of stream at closest approach (3)

       galpot - Galaxy Potential object

       GM - mass of Plummer

       rs - scale of Plummer

       tmaxfac(10) - multiple of rs/|w-v0| to use for time integration interval

       N(1000) - number of forward integration points

       integrate_method('symplec4_c') - orbit integrator to use (see Orbit.integrate)

    OUTPUT:

       deltav (nstar,3)

    HISTORY:

       2015-08-18 - SANDERS

    """
    if len(v.shape) == 1: v= numpy.reshape(v,(1,3))
    if len(x.shape) == 1: x= numpy.reshape(x,(1,3))

    nstar,ndim=numpy.shape(v)
    b0 = numpy.cross(w,v0)
    b0 *= b/numpy.sqrt(numpy.sum(b0**2))
    X = x0-b0

    # Setup Plummer orbit
    R, phi, z= bovy_coords.rect_to_cyl(X[0],X[1],X[2])
    vR, vp, vz= bovy_coords.rect_to_cyl_vec(w[0],w[1],w[2],R,phi,z,cyl=True)
    tmax = tmaxfac*rs/numpy.sqrt(numpy.sum((w-v0)**2))
    times = numpy.linspace(0.,tmax,N)
    dtimes = numpy.linspace(-tmax,tmax,2*N)
    o = Orbit(vxvv=[R,-vR,-vp,z,-vz,phi])
    o.integrate(times,galpot,method=integrate_method)
    oplum = Orbit(vxvv=[o.R(times[-1]),-o.vR(times[-1]),-o.vT(times[-1]),o.z(times[-1]),-o.vz(times[-1]),o.phi(times[-1])])
    oplum.integrate(dtimes,galpot,method=integrate_method)
    plumpot = MovingObjectPotential(orbit=oplum, GM=GM, softening_model='plummer', softening_length=rs)

    # Now integrate each particle backwards in galaxy potential, forwards in combined potential and backwards again in galaxy and take diff

    deltav = numpy.zeros((nstar,3))
    R, phi, z= bovy_coords.rect_to_cyl(x[:,0],x[:,1],x[:,2])
    vR, vp, vz= bovy_coords.rect_to_cyl_vec(v[:,0],v[:,1],v[:,2],
                                            R,phi,z,cyl=True)
    for i in range(nstar):
      ostar= Orbit(vxvv=[R[i],-vR[i],-vp[i],z[i],-vz[i],phi[i]])
      ostar.integrate(times,galpot,method=integrate_method)
      oboth = Orbit(vxvv=[ostar.R(times[-1]),-ostar.vR(times[-1]),-ostar.vT(times[-1]),ostar.z(times[-1]),-ostar.vz(times[-1]),ostar.phi(times[-1])])
      oboth.integrate(dtimes,[galpot,plumpot],method=integrate_method)
      ogalpot = Orbit(vxvv=[oboth.R(times[-1]),-oboth.vR(times[-1]),-oboth.vT(times[-1]),oboth.z(times[-1]),-oboth.vz(times[-1]),oboth.phi(times[-1])])
      ogalpot.integrate(times,galpot,method=integrate_method)
      deltav[i][0] = -ogalpot.vx(times[-1]) - v[i][0]
      deltav[i][1] = -ogalpot.vy(times[-1]) - v[i][1]
      deltav[i][2] = -ogalpot.vz(times[-1]) - v[i][2]
    return deltav

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vT [as 别名]
ps= potential.PowerSphericalPotentialwCutoff(alpha=1.8,rc=1.9/8.,normalize=0.05)
mn= potential.MiyamotoNagaiPotential(a=3./8.,b=0.28/8.,normalize=.6)
MWPotential= [ps,mn,scf]

##Creating the orbit
o = Orbit([229.018,-0.124,23.2,-2.296,-2.257,-58.7],radec=True,ro=8.,vo=220.,solarmotion=[-11.1,24.,7.25])
o.turn_physical_off()

##Integrate backwards in time
o = o.flip()
ts= nu.arange(0,(TIME + step_size).value,step_size.value)*units.Myr

o.integrate(ts, MWPotential, method='dopr54_c')

##Integrating Forward in time
newOrbit = Orbit([o.R(TIME), -o.vR(TIME), -o.vT(TIME), o.z(TIME), -o.vz(TIME), o.phi(TIME)],ro=8.,vo=220.)
newOrbit.turn_physical_off()
newOrbit.integrate(ts, MWPotential, method='dopr54_c')


def randomVelocity(std=.001):
    if type(std).__name__ == "Quantity":
        return nu.random.normal(scale=std.value)*std.unit
    return nu.random.normal(scale=std)

time1 = nu.arange(0, TIME.value, dt.value)*units.Myr
orbits_pos = nu.empty((len(time1) + 1,9,len(ts)), dtype=units.quantity.Quantity)
orbits_pos[0, :, :] = ts, newOrbit.x(ts), newOrbit.y(ts), newOrbit.z(ts), newOrbit.vx(ts), newOrbit.vy(ts), newOrbit.vz(ts), newOrbit.ra(ts), newOrbit.dec(ts)
orbits_pos[:,:,:] = orbits_pos[0,:,:]
i = 0
std = 0.004

# 需要导入模块: from galpy.orbit import Orbit [as 别名]
# 或者: from galpy.orbit.Orbit import vT [as 别名]
                dec,
                distance,
                pm_ra,
                pm_dec,
                velocity,
                orbit.x() * r0,
                orbit.y() * r0,
                orbit.z() * r0,
                orbit.U()[0],
                orbit.V()[0],
                orbit.W()[0],
                orbit.vx()[0] * v0_mc,
                orbit.vy()[0] * v0_mc,
                orbit.vz() * v0_mc,
                orbit.vR() * v0_mc,
                orbit.vT() * v0_mc,
                orbit.vphi()[0] * v0_mc,
                orbit.L()[0][2] * v0 * r0, # Lz, 18
                orbit.E(pot=MWPotential) * pow(v0, 2)/2, # Energy
                pow(np.sum(pow(orbit.L()[0][:2], 2)), 0.5) * v0 * r0 # Lperp
            ]
            
        velocities[i, :] = velocity_field


    monte_carlo_data[j, :, :] = velocities

    x_pos, y_pos = velocities[:, 6], velocities[:, 7]

    ax2.scatter(x_pos, y_pos, marker='o', s=3, edgecolor='none', facecolor=star['color'], alpha=0.3)
    invisible.append(ax2.scatter(0, 0, marker='o', s=30, edgecolor='none', facecolor=star['color'], label=star['designation']))