Skip to content

Commit a52d0a1

Browse files
Update icbuilder.py
1 parent f392cb0 commit a52d0a1

File tree

1 file changed

+51
-0
lines changed

1 file changed

+51
-0
lines changed

tools/icbuilder.py

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -377,3 +377,54 @@ def spiral_ic(N, A=1.0, B=1.0, arms=2, mass=1.0, noise=0.02):
377377
particles.append((x, y, z, vx, vy, vz, mass / N, 0))
378378

379379
return particles
380+
381+
def spiral_dm_ic(N, A=1.0, B=1.0, arms=2,
382+
mass=1.0, dm_fraction=0.85,
383+
noise=0.02, dm_scale=1.5):
384+
"""
385+
Logarithmic spiral ICs with baryons (type=0) + DM (type=1).
386+
387+
r(phi) = A * log( B * tan(phi / (2*arms)) )
388+
389+
- N: number of baryon particles (DM will also be N)
390+
- A, B: spiral parameters
391+
- arms: number of spiral arms
392+
- mass: total mass (baryons + DM)
393+
- dm_fraction: fraction of mass in dark matter
394+
- noise: positional jitter
395+
- dm_scale: DM radius multiplier (DM more extended)
396+
"""
397+
398+
import math, random
399+
particles = []
400+
401+
M_b = mass * (1 - dm_fraction)
402+
M_dm = mass * dm_fraction
403+
404+
# -----------------------
405+
# BARYONS (type = 0)
406+
# -----------------------
407+
for _ in range(N):
408+
phi = random.uniform(0.1, math.pi * 2 * arms - 0.1)
409+
r = A * math.log(B * math.tan(phi / (2 * arms)))
410+
411+
x = r * math.cos(phi) + noise * (random.random() - 0.5)
412+
y = r * math.sin(phi) + noise * (random.random() - 0.5)
413+
z = noise * (random.random() - 0.5)
414+
415+
particles.append((x, y, z, 0, 0, 0, M_b / N, 0))
416+
417+
# -----------------------
418+
# DARK MATTER (type = 1)
419+
# -----------------------
420+
for _ in range(N):
421+
phi = random.uniform(0.1, math.pi * 2 * arms - 0.1)
422+
r = dm_scale * A * math.log(B * math.tan(phi / (2 * arms)))
423+
424+
x = r * math.cos(phi) + noise * (random.random() - 0.5)
425+
y = r * math.sin(phi) + noise * (random.random() - 0.5)
426+
z = noise * (random.random() - 0.5)
427+
428+
particles.append((x, y, z, 0, 0, 0, M_dm / N, 1))
429+
430+
return particles

0 commit comments

Comments
 (0)