@@ -330,4 +330,78 @@ def random_solar_system(
330330
331331 particles .append ((mx , my , mz , mvx , mvy , mvz , mm ))
332332
333+ return particles
334+
335+ def big_bang_ic (
336+ N ,
337+ radius = 1.0 ,
338+ mass = 1.0 ,
339+ hubble_k = 1.0 ,
340+ center_bias = 1.0 ,
341+ perturb_amp = 0.05 ,
342+ vel_jitter = 0.02
343+ ):
344+ """
345+ Big-Bang-like initial conditions:
346+ - Expanding sphere
347+ - Small density perturbations
348+ - Radial Hubble flow
349+ - Isotropic velocity jitter
350+ """
351+
352+ particles = []
353+ m = mass / N
354+
355+ for _ in range (N ):
356+ # --- 1. Sample radius with optional central bias ---
357+ u = random .random ()
358+ r = radius * (u ** (1.0 / (3.0 * center_bias )))
359+
360+ # Random isotropic direction
361+ z = 2 * random .random () - 1
362+ t = 2 * math .pi * random .random ()
363+ s = math .sqrt (1 - z * z )
364+
365+ ux = s * math .cos (t )
366+ uy = s * math .sin (t )
367+ uz = z
368+
369+ # Position
370+ x = r * ux
371+ y = r * uy
372+ z = r * uz
373+
374+ # --- 2. Add density perturbations ---
375+ # Small random displacement to seed structure formation
376+ dx = perturb_amp * radius * (random .random () - 0.5 )
377+ dy = perturb_amp * radius * (random .random () - 0.5 )
378+ dz = perturb_amp * radius * (random .random () - 0.5 )
379+
380+ x += dx
381+ y += dy
382+ z += dz
383+
384+ # --- 3. Hubble-like expansion ---
385+ v = hubble_k * r
386+
387+ vx = v * ux
388+ vy = v * uy
389+ vz = v * uz
390+
391+ # --- 4. Isotropic velocity jitter ---
392+ # Proper spherical noise
393+ jz = 2 * random .random () - 1
394+ jt = 2 * math .pi * random .random ()
395+ js = math .sqrt (1 - jz * jz )
396+
397+ jx = js * math .cos (jt )
398+ jy = js * math .sin (jt )
399+ jz = jz
400+
401+ vx += vel_jitter * v * jx
402+ vy += vel_jitter * v * jy
403+ vz += vel_jitter * v * jz
404+
405+ particles .append ((x , y , z , vx , vy , vz , m ))
406+
333407 return particles
0 commit comments