Skip to content

Commit dfb551a

Browse files
Added more features to the IC-Builder
1 parent c55ead7 commit dfb551a

1 file changed

Lines changed: 70 additions & 0 deletions

File tree

tools/icbuilder.py

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,4 +63,74 @@ def disk(N, radius=10.0, mass=1.0, thickness=0.1):
6363

6464
particles.append((x, y, z, vx, vy, vz, mass / N))
6565

66+
return particles
67+
68+
def hernquist(N, scale=1.0, mass=1.0):
69+
"""
70+
Generates a Hernquist sphere with scale radius 'scale' and total mass 'mass'.
71+
Velocities are set to zero (cold start).
72+
"""
73+
particles = []
74+
M = mass
75+
a = scale
76+
77+
for _ in range(N):
78+
# Sample radius from Hernquist distribution
79+
u = random.random()
80+
r = a * math.sqrt(u) / (1 - math.sqrt(u))
81+
82+
# Random isotropic angles
83+
theta = math.acos(2*random.random() - 1)
84+
phi = 2 * math.pi * random.random()
85+
86+
# Convert to Cartesian
87+
x = r * math.sin(theta) * math.cos(phi)
88+
y = r * math.sin(theta) * math.sin(phi)
89+
z = r * math.cos(theta)
90+
91+
particles.append((x, y, z, 0, 0, 0, M / N))
92+
93+
return particles
94+
95+
def rotate_z(x, y, z, angle):
96+
ca = math.cos(angle)
97+
sa = math.sin(angle)
98+
return (ca*x - sa*y, sa*x + ca*y, z)
99+
100+
def two_galaxies(
101+
N1, R1, M1, scale1,
102+
N2, R2, M2, scale2,
103+
angle_deg=0.0,
104+
v_rel=1.0
105+
):
106+
"""
107+
Builds a two-galaxy merger:
108+
- Galaxy 1: Hernquist sphere at (-R1, 0, 0)
109+
- Galaxy 2: Hernquist sphere at (+R2, 0, 0), rotated by angle_deg
110+
- v_rel: relative velocity along y-axis
111+
"""
112+
113+
g1 = hernquist(N1, scale1, M1)
114+
g2 = hernquist(N2, scale2, M2)
115+
116+
angle = math.radians(angle_deg)
117+
particles = []
118+
119+
# Galaxy 1: shift left, give velocity +v_rel/2
120+
for (x, y, z, vx, vy, vz, m) in g1:
121+
particles.append((
122+
x - R1, y, z,
123+
vx, vy + v_rel/2, vz,
124+
m
125+
))
126+
127+
# Galaxy 2: rotate, shift right, give velocity -v_rel/2
128+
for (x, y, z, vx, vy, vz, m) in g2:
129+
xr, yr, zr = rotate_z(x, y, z, angle)
130+
particles.append((
131+
xr + R2, yr, zr,
132+
vx, vy - v_rel/2, vz,
133+
m
134+
))
135+
66136
return particles

0 commit comments

Comments
 (0)