<jupyter_start><jupyter_text>Simplified, idealized simulations to start: all at the same time, no forward evolution To perform some simple simulations to start we will do the following:        * We'll assume that the stream is frozen at the time of impact (i.e., no new stars are added after the impacts, which would dilute the signal)    * We still use a single mass for each mass-bin (and a single scale radius)    * We include hits with impact parameters between 0 and 2 x the scale radius (Carlberg 2012)    * We draw subhalo velocities from a Gaussian with an isotropic dispersion of 150 km/s    * All impacts happen at the same time    * We only look at the stream at the time of impact    We will load the galpy model for the gap in Sanders, Bovy, & Erkal (2015) as the stream model, hitting the stream at 10 Gyr<jupyter_code>lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
xv_prog_init= numpy.array([30.,0.,0.,0., 105.74895, 105.74895])
def convert_to_cylindrical(xv):
    R,phi,Z= bovy_coords.rect_to_cyl(xv[:,0],xv[:,1],xv[:,2])
    vR,vT,vZ= bovy_coords.rect_to_cyl_vec(xv[:,3],xv[:,4],xv[:,5],R,phi,Z,cyl=True)
    out= numpy.empty_like(xv)
    # Preferred galpy arrangement of cylindrical coordinates
    out[:,0]= R
    out[:,1]= vR
    out[:,2]= vT
    out[:,3]= Z
    out[:,4]= vZ
    out[:,5]= phi
    return out
sigv= 0.365*(10./2.)**(1./3.)
RvR_prog_init= convert_to_cylindrical(xv_prog_init[:,numpy.newaxis].T)[0,:]
prog_init= Orbit([RvR_prog_init[0]/R0,RvR_prog_init[1]/V0,RvR_prog_init[2]/V0,
                  RvR_prog_init[3]/R0,RvR_prog_init[4]/V0,RvR_prog_init[5]])
times= numpy.linspace(0.,10.88/bovy_conversion.time_in_Gyr(V0,R0),10001)
prog_init.integrate(times,lp)
xv_unp_peri_prog= [prog_init.x(times[-1])[0],prog_init.y(times[-1])[0],prog_init.z(times[-1]),
              prog_init.vx(times[-1])[0],prog_init.vy(times[-1])[0],prog_init.vz(times[-1])]
RvR_unp_peri_prog= convert_to_cylindrical(numpy.array(xv_unp_peri_prog)[:,numpy.newaxis].T)[0,:]
prog_unp_peri= Orbit([RvR_unp_peri_prog[0],RvR_unp_peri_prog[1],RvR_unp_peri_prog[2],
                      RvR_unp_peri_prog[3]+0.0,RvR_unp_peri_prog[4],RvR_unp_peri_prog[5]])
aAI= actionAngleIsochroneApprox(pot=lp,b=0.8)
sdfl_per_1e8_peri= galpy.df.streamgapdf(sigv/V0,progenitor=prog_unp_peri,pot=lp,aA=aAI,
                                        leading=False,nTrackChunks=26,nTrackIterations=1,
                                        sigMeanOffset=4.5,
                                        tdisrupt=10.88/bovy_conversion.time_in_Gyr(V0,R0),
                                        Vnorm=V0,Rnorm=R0,
                                        impactb=0.,
                                        subhalovel=numpy.array([6.82200571,132.7700529,149.4174464])/V0, # from Denis' email
                                        timpact=0.88/bovy_conversion.time_in_Gyr(V0,R0),
                                        impact_angle=-2.34,
                                        GM=10.**-2./bovy_conversion.mass_in_1010msol(V0,R0),
                                        rs=0.625/R0)<jupyter_output>galpyWarning: Using C implementation to integrate orbits<jupyter_text>Functions that give the scale radius and number of hits as a function of mass<jupyter_code>def rs(M):
    """Give scale radius for a given mass"""
    return (M/10**7.)**0.4*0.25/R0
def N(M):
    return 10.*(10.**7./M)<jupyter_output><empty_output><jupyter_text>The following function computes the velocity kicks due to a random set of impacts of a given mass, happening at a given rate<jupyter_code>thetapars= sdfl_per_1e8_peri._kick_interpolatedThetasTrack
def randomdvs(GM,sdf,rate=None):
    """Construct a random set of impacts; inputs:
    
    M: mass of subhalo (/Msolar)
    sdf: a streamgapdf object that we can use to compute kicks
    rate: expected number of impacts
    """
    if rate is None: rate= N(GM)
    # First we sample angles at which the impacts occur, starting at the progenitor
    angles= numpy.cumsum(numpy.random.exponential(scale=1./rate,size=int(numpy.ceil(rate))))
    angles= angles[angles < 1.]*sdf._deltaAngleTrackImpact*sdf._gap_sigMeanSign
    nimpacts= len(angles)
    # Now we sample for each impact a set of impact parameters and halo velocities and compute dv
    dvs= numpy.zeros((nimpacts,len(thetapars),3)) 
    dvdotvs= numpy.zeros((nimpacts,len(thetapars)))
    for ii, angle in enumerate(angles):
        # Sample b from a linearly increasing probability
        b= numpy.random.uniform()**0.5*2.*rs(GM)
        if numpy.random.uniform() < 0.5: b*= -1.
        v= numpy.random.normal(scale=150./V0,size=3)
        # Get the stream at the closest-approach point
        kick_ObsTrackXY_closest= numpy.array([\
                sdf._kick_interpTrackX(angle),
                sdf._kick_interpTrackY(angle),
                sdf._kick_interpTrackZ(angle),
                sdf._kick_interpTrackvX(angle),
                sdf._kick_interpTrackvY(angle),
                sdf._kick_interpTrackvZ(angle)])
        dvs[ii]= impulse_deltav_plummer_curvedstream(
            sdf._kick_interpolatedObsTrackXY[:,3:],
            sdf._kick_interpolatedObsTrackXY[:,:3],
            b,v,
            kick_ObsTrackXY_closest[:3],
            kick_ObsTrackXY_closest[3:],
            GM/bovy_conversion.mass_in_msol(V0,R0),rs(GM))
        dvdotvs[ii]= numpy.dot(dvs[ii],kick_ObsTrackXY_closest[3:])
    return (dvs,dvdotvs)<jupyter_output><empty_output><jupyter_text>Single mass ranges Let's look at delta vs<jupyter_code>figsize(16,4)
subplot(1,4,1)
dvs, dvdotvs= randomdvs(10**8,sdfl_per_1e8_peri)
plot(thetapars,dvdotvs[0,:])
plot(thetapars,numpy.sum(dvdotvs,axis=0),zorder=0,color='k')
xlabel(r'$\theta_\parallel$')
ylabel(r'$\delta v\cdot v$')
ylim(-0.005,0.005)
bovy_plot.bovy_text(r'$M = 10^8\,M_\odot$',top_left=True)
subplot(1,4,2)
dvs, dvdotvs= randomdvs(10**7,sdfl_per_1e8_peri)
plot(thetapars,dvdotvs[0,:])
plot(thetapars,dvdotvs[3,:])
plot(thetapars,dvdotvs[5,:])
plot(thetapars,numpy.sum(dvdotvs,axis=0),zorder=0,color='k')
xlabel(r'$\theta_\parallel$')
ylim(-0.005,0.005)
bovy_plot.bovy_text(r'$M = 10^7\,M_\odot$',top_left=True)
subplot(1,4,3)
dvs, dvdotvs= randomdvs(10**6,sdfl_per_1e8_peri)
plot(thetapars,dvdotvs[0,:])
plot(thetapars,dvdotvs[30,:])
plot(thetapars,dvdotvs[50,:])
plot(thetapars,numpy.sum(dvdotvs,axis=0),zorder=0,color='k')
xlabel(r'$\theta_\parallel$')
ylim(-0.005,0.005)
bovy_plot.bovy_text(r'$M = 10^6\,M_\odot$',top_left=True)
subplot(1,4,4)
dvs, dvdotvs= randomdvs(10**5,sdfl_per_1e8_peri)
plot(thetapars,dvdotvs[0,:])
plot(thetapars,dvdotvs[300,:])
plot(thetapars,dvdotvs[500,:])
plot(thetapars,numpy.sum(dvdotvs,axis=0),zorder=0,color='k')
xlabel(r'$\theta_\parallel$')
ylim(-0.005,0.005)
bovy_plot.bovy_text(r'$M = 10^5\,M_\odot$',top_left=True)
tight_layout()<jupyter_output><empty_output><jupyter_text>The following figure shows simulations where all mass ranges hit at the same time:<jupyter_code>figsize(16,4)
subplot(1,4,1)
dvs, dvdotvs8= randomdvs(10**8,sdfl_per_1e8_peri)
dvs, dvdotvs7= randomdvs(10**7,sdfl_per_1e8_peri)
dvs, dvdotvs6= randomdvs(10**6,sdfl_per_1e8_peri)
dvs, dvdotvs5= randomdvs(10**5,sdfl_per_1e8_peri)
plot(thetapars,numpy.sum(dvdotvs8,axis=0)
     +numpy.sum(dvdotvs7,axis=0)
     +numpy.sum(dvdotvs6,axis=0)
     +numpy.sum(dvdotvs5,axis=0),
     zorder=0,color='k')
xlabel(r'$\theta_\parallel$')
ylabel(r'$\delta v\cdot v$')
ylim(-0.005,0.005)
subplot(1,4,2)
dvs, dvdotvs8= randomdvs(10**8,sdfl_per_1e8_peri)
dvs, dvdotvs7= randomdvs(10**7,sdfl_per_1e8_peri)
dvs, dvdotvs6= randomdvs(10**6,sdfl_per_1e8_peri)
dvs, dvdotvs5= randomdvs(10**5,sdfl_per_1e8_peri)
plot(thetapars,numpy.sum(dvdotvs8,axis=0)
     +numpy.sum(dvdotvs7,axis=0)
     +numpy.sum(dvdotvs6,axis=0)
     +numpy.sum(dvdotvs5,axis=0),
     zorder=0,color='k')
xlabel(r'$\theta_\parallel$')
ylabel(r'$\delta v\cdot v$')
ylim(-0.005,0.005)
subplot(1,4,3)
dvs, dvdotvs8= randomdvs(10**8,sdfl_per_1e8_peri)
dvs, dvdotvs7= randomdvs(10**7,sdfl_per_1e8_peri)
dvs, dvdotvs6= randomdvs(10**6,sdfl_per_1e8_peri)
dvs, dvdotvs5= randomdvs(10**5,sdfl_per_1e8_peri)
plot(thetapars,numpy.sum(dvdotvs8,axis=0)
     +numpy.sum(dvdotvs7,axis=0)
     +numpy.sum(dvdotvs6,axis=0)
     +numpy.sum(dvdotvs5,axis=0),
     zorder=0,color='k')
xlabel(r'$\theta_\parallel$')
ylabel(r'$\delta v\cdot v$')
ylim(-0.005,0.005)
subplot(1,4,4)
dvs, dvdotvs8= randomdvs(10**8,sdfl_per_1e8_peri)
dvs, dvdotvs7= randomdvs(10**7,sdfl_per_1e8_peri)
dvs, dvdotvs6= randomdvs(10**6,sdfl_per_1e8_peri)
dvs, dvdotvs5= randomdvs(10**5,sdfl_per_1e8_peri)
plot(thetapars,numpy.sum(dvdotvs8,axis=0)
     +numpy.sum(dvdotvs7,axis=0)
     +numpy.sum(dvdotvs6,axis=0)
     +numpy.sum(dvdotvs5,axis=0),
     zorder=0,color='k')
xlabel(r'$\theta_\parallel$')
ylabel(r'$\delta v\cdot v$')
ylim(-0.005,0.005)
tight_layout()<jupyter_output><empty_output><jupyter_text>Now we look at the (median) power spectra of these deviations. The black line on the left is when hitting the stream with all mass ranges:<jupyter_code>figsize(16,8)
ntrials= 1001
# Just to get the sampling frequencies
px= signal.welch(numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[0]
ppy= numpy.zeros((len(px),ntrials))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**5.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
subplot(1,2,1)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)))
subplot(1,2,2)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)/numpy.amax(numpy.median(ppy,axis=1))))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**6.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
subplot(1,2,1)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)))
subplot(1,2,2)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)/numpy.amax(numpy.median(ppy,axis=1))))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
subplot(1,2,1)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)))
subplot(1,2,2)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)/numpy.amax(numpy.median(ppy,axis=1))))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**7.85,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
subplot(1,2,1)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**5.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**6.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**8.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)),'k')
ylim(10**-5.,10.**-3.)
subplot(1,2,2)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)/numpy.amax(numpy.median(ppy,axis=1))))
ylim(10**-2.,1.)<jupyter_output><empty_output><jupyter_text>More realistic subhalo numbers (assume N~1/sqrt(M) to take into account shrinking sphere of influence)<jupyter_code>def N(M):
    return 10.*(10.**7./M)**0.5
figsize(16,8)
ntrials= 1001
# Just to get the sampling frequencies
px= signal.welch(numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[0]
ppy= numpy.zeros((len(px),ntrials))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**5.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
subplot(1,2,1)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)))
subplot(1,2,2)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)/numpy.amax(numpy.median(ppy,axis=1))))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**6.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
subplot(1,2,1)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)))
subplot(1,2,2)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)/numpy.amax(numpy.median(ppy,axis=1))))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
subplot(1,2,1)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)))
subplot(1,2,2)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)/numpy.amax(numpy.median(ppy,axis=1))))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**7.85,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
subplot(1,2,1)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)))
for ii in range(ppy.shape[1]):
    ppy[:,ii]= signal.welch(numpy.sum(randomdvs(10**5.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**6.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**8.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)),'k')
ylim(10**-6.,10.**-3.)
subplot(1,2,2)
loglog(px,numpy.sqrt(numpy.median(ppy,axis=1)/numpy.amax(numpy.median(ppy,axis=1))))
ylim(10**-3.,1.)
# Go back to original definition
def N(M):
    return 10.*(10.**7./M)<jupyter_output><empty_output><jupyter_text>Back to N~M, but showing individual streams<jupyter_code>figsize(16,8)
ntrials= 2
plotkwargs= {'alpha':.9,'lw':1.}
# Just to get the sampling frequencies
px= signal.welch(numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[0]
for ii in range(ntrials):
    ppy= signal.welch(numpy.sum(randomdvs(10**5.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
    subplot(1,2,1)
    loglog(px,numpy.sqrt(ppy),color=sns.color_palette()[0],**plotkwargs)
    subplot(1,2,2)
    loglog(px,numpy.sqrt(ppy/numpy.amax(ppy)),color=sns.color_palette()[0],**plotkwargs)
for ii in range(ntrials):
    ppy= signal.welch(numpy.sum(randomdvs(10**6.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
    subplot(1,2,1)
    loglog(px,numpy.sqrt(ppy),color=sns.color_palette()[1],**plotkwargs)
    subplot(1,2,2)
    loglog(px,numpy.sqrt(ppy/numpy.amax(ppy)),color=sns.color_palette()[1],**plotkwargs)
for ii in range(ntrials):
    ppy= signal.welch(numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
    subplot(1,2,1)
    loglog(px,numpy.sqrt(ppy),color=sns.color_palette()[2],**plotkwargs)
    subplot(1,2,2)
    loglog(px,numpy.sqrt(ppy/numpy.amax(ppy)),color=sns.color_palette()[2],**plotkwargs)
for ii in range(ntrials):
    ppy= signal.welch(numpy.sum(randomdvs(10**5.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**6.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**7.,sdfl_per_1e8_peri)[1],axis=0)
                            +numpy.sum(randomdvs(10**8.,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
    subplot(1,2,1)
    loglog(px,numpy.sqrt(ppy),color='k',**plotkwargs)
for ii in range(ntrials):
    ppy= signal.welch(numpy.sum(randomdvs(10**7.85,sdfl_per_1e8_peri)[1],axis=0),
                                  fs=1./(thetapars[1]-thetapars[0]),scaling='spectrum')[1]
    subplot(1,2,1)
    loglog(px,numpy.sqrt(ppy),color=sns.color_palette()[3],**plotkwargs)
    ylim(10**-6.,10.**-3.)
    subplot(1,2,2)
    loglog(px,numpy.sqrt(ppy/numpy.amax(ppy)),color=sns.color_palette()[3],**plotkwargs)
    ylim(10**-3.,1.)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1742304416595693, 'fraction_numerical': 0.05730624163523543, 'mean_word_length': 4.080037082818294, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '314784', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7387, 'n_tokens_neox': 6864, 'n_words': 818, 'path': 'dev/bell.ipynb', 'repo_name': 'jobovy/streamgap-pepper', 'size': '16438'}