@@ -244,4 +244,90 @@ def four_body_cross(mass=1.0):
244244 (x2 , y2 , 0.0 , vx2 , vy2 , 0.0 , m ),
245245 (x3 , y3 , 0.0 , vx3 , vy3 , 0.0 , m ),
246246 (x4 , y4 , 0.0 , vx4 , vy4 , 0.0 , m ),
247- ]
247+ ]
248+
249+ def random_solar_system (
250+ seed_value = None ,
251+ star_mass_range = (0.5 , 2.0 ),
252+ planet_count_range = (1 , 12 ),
253+ planet_mass_range = (1e-6 , 1e-3 ),
254+ semi_major_axis_range = (0.3 , 40.0 ),
255+ eccentricity_range = (0.0 , 0.3 ),
256+ inclination_range = (0.0 , 5.0 ), # degrees
257+ moon_chance = 0.25 ,
258+ max_moons = 3
259+ ):
260+ """
261+ Generates a random solar system using Keplerian orbits.
262+ Returns a list of (x, y, z, vx, vy, vz, m) particles.
263+ """
264+
265+ if seed_value is not None :
266+ random .seed (seed_value )
267+
268+ particles = []
269+
270+ # --- 1. STAR ---
271+ star_mass = random .uniform (* star_mass_range )
272+ particles .append ((0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , star_mass ))
273+
274+ # --- 2. PLANETS ---
275+ num_planets = random .randint (* planet_count_range )
276+
277+ for _ in range (num_planets ):
278+ m = random .uniform (* planet_mass_range )
279+
280+ # Orbital parameters
281+ a = random .uniform (* semi_major_axis_range ) # semi-major axis
282+ e = random .uniform (* eccentricity_range )
283+ inc = math .radians (random .uniform (* inclination_range ))
284+ phase = random .uniform (0 , 2 * math .pi )
285+
286+ # Distance at given phase
287+ r = a * (1 - e ** 2 ) / (1 + e * math .cos (phase ))
288+
289+ # Position in orbital plane
290+ x = r * math .cos (phase )
291+ y = r * math .sin (phase )
292+ z = 0.0
293+
294+ # Rotate by inclination
295+ z = x * math .sin (inc )
296+ x = x * math .cos (inc )
297+
298+ # Orbital velocity magnitude (vis-viva)
299+ v = math .sqrt (star_mass * (2 / r - 1 / a ))
300+
301+ # Velocity direction (perpendicular to radius vector)
302+ vx = - v * math .sin (phase )
303+ vy = v * math .cos (phase )
304+ vz = 0.0
305+
306+ # Rotate velocity by inclination
307+ vz = vx * math .sin (inc )
308+ vx = vx * math .cos (inc )
309+
310+ particles .append ((x , y , z , vx , vy , vz , m ))
311+
312+ # --- 3. MOONS ---
313+ if random .random () < moon_chance :
314+ moon_count = random .randint (1 , max_moons )
315+ for _ in range (moon_count ):
316+ mm = m * random .uniform (0.001 , 0.05 )
317+ da = random .uniform (0.01 , 0.2 ) # moon orbit radius
318+ phase_m = random .uniform (0 , 2 * math .pi )
319+
320+ # Moon position relative to planet
321+ mx = x + da * math .cos (phase_m )
322+ my = y + da * math .sin (phase_m )
323+ mz = z
324+
325+ # Moon velocity relative to planet
326+ vm = math .sqrt (m / da )
327+ mvx = vx - vm * math .sin (phase_m )
328+ mvy = vy + vm * math .cos (phase_m )
329+ mvz = vz
330+
331+ particles .append ((mx , my , mz , mvx , mvy , mvz , mm ))
332+
333+ return particles
0 commit comments