Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 66 additions & 19 deletions garnets/garnets.py
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,10 @@ def coalesce_planetesimals(disk, planets, canidate, do_moons):
and canidate.mass > .0001 * earth_mass \
and planet.mass_of_moons < planet.mass * .05 \
and planet.mass > canidate.mass:
# TODO: Remove planet.mass > canidate.mass distinction, just switch the canidate and planet!
# TODO: Refactor moon capture logic. If `canidate.mass > planet.mass`,
# swap `candidate` and `planet` to ensure the more massive body is treated
# as the primary object capturing the smaller one as a moon.
# This avoids the `planet.mass > canidate.mass` check.
planet.add_moon(
convert_planetesimal_to_protomoon(
canidate, planet))
Expand All @@ -219,9 +222,11 @@ def coalesce_planetesimals(disk, planets, canidate, do_moons):
finished = True
break
else:
# TODO: Reasons.
# TODO: Enhance logging for failed moon capture. Specify which condition(s)
# (e.g., candidate mass too high/low, planet's existing moon mass too high,
# candidate more massive than planet) prevented capture.
logging.info(
"Did not capture potential moon at %s. Collision imminent." % quantity_repr(
"Did not capture potential moon at %s. Conditions not met or collision imminent." % quantity_repr(
planet.orbit.a, au)
)

Expand All @@ -234,7 +239,11 @@ def coalesce_planetesimals(disk, planets, canidate, do_moons):
planet.gas_mass = planet.gas_mass + canidate.gas_mass # + new_gas
finished = True

# Accrete MORE DUST! TODO: Refactor to this, mark the planet as needing to re-accrete.
# TODO: Consider refactoring the post-collision accretion step.
# Instead of immediate re-accretion via `disk.accrete_dust(planet)`,
# explore marking the planet as 'dirty' or 'needs_re_accretion' to
# handle this in a subsequent simulation phase. This could allow for
# more complex interactions or batch processing of accretion events.
disk.accrete_dust(planet)

logging.info("Conglomerate is now %s at %s." % (
Expand All @@ -243,9 +252,11 @@ def coalesce_planetesimals(disk, planets, canidate, do_moons):
))

if not finished:
# TODO: Extra info.
logging.info("New Protoplanet at %s." %
quantity_repr(canidate.orbit.a, au))
# TODO: Enhance logging for new protoplanet formation. Include
# additional information such as its mass and potentially key
# compositional details (e.g., dust/gas ratio if available).
logging.info("New Protoplanet formed at %s with mass %s." %
(quantity_repr(canidate.orbit.a, au), mass_repr(canidate.mass)))
planets.append(convert_planetesimal_to_protoplanet(canidate))


Expand All @@ -259,7 +270,11 @@ def calculate_gases(star, planet, planet_id):

for i in range(len(GASES)):

# TODO(woursler): WTH is this?
# TODO: Clarify the derivation and meaning of the variable `yp`
# and the formula used for its calculation. This formula appears to adjust
# boiling points based on pressure. Document the source of the constants
# (e.g., -5050.5, 373., 0.001) and explain its physical basis
# (e.g., Clausius-Clapeyron approximation for water, adapted for other gases).
yp = GASES[i].boil / \
(373. * ((log((pressure) + 0.001) / -5050.5) + (1.0 / 373.)))

Expand All @@ -274,7 +289,11 @@ def calculate_gases(star, planet, planet_id):
fract = 1.0
pres2 = 1.0

# TODO(woursler): This needs to go in the chemtable somehow.
# TODO: Refactor gas reactivity logic. Move the hardcoded rules
# for specific gases (Ar, He, O, O2, CO2) and their dependence on
# stellar age and pressure into `chemtable.py`, possibly by adding
# reactivity parameters or functions to the `Gas` class or a
# related data structure. This would centralize chemical properties.

if GASES[i].symbol == "Ar":
react = .15 * star.age / 4e9 / year
Expand Down Expand Up @@ -383,7 +402,10 @@ def roche_limit(planet, moon):


def hill_sphere(planet, star):
# TODO(woursler): Confirm this
# TODO: The formula used is a common approximation for the Hill sphere
# radius ( R_H = a * (m / (3*M))^(1/3) ). Confirm if this level of precision
# is adequate for the simulation or if more exact formulations are needed for
# specific scenarios (e.g., highly eccentric orbits, larger moon/planet mass ratios).
return planet.orbit.a * pow(
(planet.mass / (3.0 * star.mass)), (1.0 / 3.0))

Expand Down Expand Up @@ -431,12 +453,19 @@ def generate_planet(protoplanet,
planet.radius = volume_radius(planet.mass, planet.density)

planet.surf_accel = acceleration(planet.mass, planet.radius)
planet.surf_grav = planet.surf_accel # TODO(woursler): Remove ambiguity.
# TODO: Clarify the distinction between `surf_accel` (calculated
# from M, R) and `surf_grav`. If they are intended to be identical
# (gravitational acceleration at the surface), consider using only one
# variable (e.g., `surface_gravity`). If `surf_grav` might later include
# other effects (e.g., centrifugal due to rotation), document this.
planet.surf_grav = planet.surf_accel

planet.molec_weight = min_molec_weight(planet)

if ((planet.mass > 1 * earth_mass)
# TODO(woursler): Understand where this comes from.
# TODO: Document the source or rationale for the 0.05 (5%)
# gas mass fraction threshold used in classifying gas giant types.
# Cite relevant planetary science literature or model assumptions.
and ((planet.gas_mass / planet.mass) > 0.05)
and (planet.molec_weight <= lookup_gas('He').weight)):

Expand Down Expand Up @@ -555,6 +584,11 @@ def generate_planet(protoplanet,
),
atm,
)'''
# TODO: Review the surface pressure calculation for rocky planets.
# Currently, it's set to a fixed `1 * atm` before `iterate_surface_temp`.
# Evaluate if the commented-out detailed `pressure(...)` calculation
# (based on `volatile_gas_inventory`) should be re-enabled and integrated
# into the iteration loop for a self-consistent atmosphere model.

if planet.surf_pressure == 0 * atm:
planet.boil_point = 0 * K
Expand All @@ -576,8 +610,12 @@ def generate_planet(protoplanet,

# Next we assign a type to the planet.

# TODO(woursler): Figure out if the pressure unit is atm or millibar.
if (planet.surf_pressure < 1 * millibar):
# TODO: Critical: Ensure consistent units for `planet.surf_pressure`.
# It's initialized in `atm` (e.g. `1 * atm` or `INCREDIBLY_LARGE_NUMBER * atm`)
# but compared with `millibar` values in type classification (e.g. `< 1 * millibar`).
# Convert values to a consistent unit (e.g., millibar or atm) before
# comparison to ensure correct planet type assignment.
if (planet.surf_pressure < 1 * millibar): # Assuming surf_pressure is converted to mbar for this block
if (not is_moon) and planet.mass < ASTEROID_MASS_LIMIT:
planet.type = PlanetType.ASTERIODS
else:
Expand Down Expand Up @@ -615,7 +653,10 @@ def generate_planet(protoplanet,
elif (planet.surf_temp < FREEZING_POINT_OF_WATER):
planet.type = PlanetType.ICE
else:
# TODO(woursler): Consider throwing an error here.
# TODO: Implement error handling for unclassified planets.
# If planet type remains `UNKNOWN` after classification logic,
# raise an error or log a critical warning to identify and address
# gaps in the classification rules.
planet.type = PlanetType.UNKNOWN
'''if (flag_verbose & 0x0001)
fprintf (stderr, "%12s\tp=%4.2Lf\tm=%4.2Lf\tg=%4.2Lf\tt=%+.1Lf\t%s\t Unknown %s\n",
Expand Down Expand Up @@ -656,9 +697,12 @@ def generate_planet(protoplanet,
moon.orbit = Orbit(a=moon_a, e=moon_e)

else:
# TODO(woursler): This seems like an error?
# The moon has no stable space to orbit.
moon.orbit = Orbit(a=0, e=0)
# TODO: Address unstable moon orbits. If
# `roche_limit_r * 1.5 >= hill_sphere_r / 2`, setting orbit to a=0, e=0
# is not a physical solution. Consider alternative outcomes: moon is
# destroyed (becomes rings), ejected, or crashes into the planet.
# Log this event clearly.
moon.orbit = Orbit(a=0*km, e=0) # Representing no stable orbit found

planet.moons.append(moon)

Expand Down Expand Up @@ -713,7 +757,10 @@ def generate_planet(protoplanet,

print(system)

# Output SVG TODO(woursler): Configure as flags?
# Output SVG
# TODO: Make SVG output configurable. Add command-line flags
# (e.g., using `argparse`) to control whether the SVG visualization is
# generated and to specify the output path.

max_x = 1500
max_y = 120
Expand Down
69 changes: 51 additions & 18 deletions garnets/stellar_system.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,10 @@ def mass(self):

@property
# Approximates the luminosity of the star.
# TODO: express only as ratio?
# TODO: Clarify if this property should return a pre-calculated ratio or if the
# current calculation (which results in a de facto ratio to Solar luminosity)
# is the intended design. If the latter, consider renaming to
# `calculated_luminosity_ratio` or adding a comment explaining units.
# Source: http://en.wikipedia.org/wiki/Mass%E2%80%93luminosity_relation
def luminosity_ratio(self):
if (self.mass_ratio < .43):
Expand All @@ -55,12 +58,18 @@ def luminosity_ratio(self):
return 3200 * self.mass_ratio

@property
# Source: StarGen, TODO Verify against current data.
# Source: StarGen
# TODO: Verify the formula `200 * (self.mass_ratio**(1/3)) * au` for stellar
# dust limit against current astrophysical models and data, as its source is StarGen.
def stellar_dust_limit(self):
return 200 * (self.mass_ratio**(1/3)) * au

@property
def r_ecosphere(self): # Source: StarGen, TODO Name? Value? Possible habitable zone?
# Source: StarGen
# TODO: Verify the term 'r_ecosphere' and its calculation
# `sqrt(self.luminosity_ratio) * au` against modern definitions of the habitable
# zone center. Confirm if this simplified model is appropriate for the simulation's scope.
def r_ecosphere(self):
return sqrt(self.luminosity_ratio) * au

@property
Expand All @@ -72,7 +81,11 @@ def max_r_ecosphere(self):
return sqrt(self.luminosity_ratio / 0.48) * au

@property
def life(self): # Source: StarGen, TODO Name? Value?
# Source: StarGen
# TODO: Verify the formula `10**10 * (self.mass_ratio / self.luminosity_ratio)`
# for stellar lifetime against current astrophysical models. Ensure the
# implicit unit (years) is clear or made explicit by returning a Quantity.
def life(self):
return 10**10 * (self.mass_ratio / self.luminosity_ratio)

@property
Expand Down Expand Up @@ -112,7 +125,11 @@ def apoapsis(self):
return (1 + self.e) * self.a


STAR_MASS = 1 * solar_mass # TODO(woursler): Use the actual star mass
# TODO: Critical: Replace global `STAR_MASS` constant with the actual mass of
# the star associated with the planetoid (e.g., `self.disk.star.mass` for
# Planetesimals or `self.star.mass` for Protoplanets) in the `reduced_mass`
# calculation for accuracy. This constant is a placeholder.
STAR_MASS = 1 * solar_mass


@attrs
Expand All @@ -130,19 +147,24 @@ def reduced_mass(self):
# To understand what this is all about...
# http://spiff.rit.edu/classes/phys440/lectures/reduced/reduced.html
# But some sort of 3 body case, see dole.
# TODO: Understand better?
# TODO(woursler): Actually use the mass of the star?
# TODO: Investigate and document the origin and physical basis of the
# `reduced_mass` formula `(self.mass / (STAR_MASS + self.mass))**0.25`,
# particularly the 0.25 exponent. Cite Dole or other relevant sources.
# This is also linked to the TODO regarding the global `STAR_MASS`.
return (self.mass / (STAR_MASS + self.mass))**0.25

@property
def inner_effect_limit(self):
# TODO(woursler): Should extend further in, include dust destabilized
# by indirect gravitational effects.
# TODO: Enhance `inner_effect_limit` model. Consider extending the limit
# further inward to account for dust destabilization by indirect gravitational
# effects, beyond the current geometric calculation.
return self.orbit.a * (1.0 - self.orbit.e) / (1.0 + DISK_ECCENTRICITY)

@property
def outer_effect_limit(self):
# TODO(woursler): See note on inner_effect_limit
# TODO: Enhance `outer_effect_limit` model, similar to `inner_effect_limit`.
# Consider effects beyond the current geometric calculation (e.g., indirect
# gravitational perturbations) that might extend the zone of influence.
return self.orbit.a * (1.0 + self.orbit.e) / (1.0 - DISK_ECCENTRICITY)


Expand All @@ -154,11 +176,14 @@ class Planetesimal(Planetoid):
def critical_mass(self):
perihelion_dist = self.orbit.a * (1 - self.orbit.e)
temp = perihelion_dist * sqrt(self.disk.star.luminosity_ratio)
# TODO(woursler): Understand the basis for this.
# TODO: Document the origin and theoretical basis of the `critical_mass`
# formula `B * (dimensionless_with_units(temp, au)**-0.75) * solar_mass`.
# Cite the relevant planetary formation theory or source (e.g., Dole, Hayashi, etc.)
# for this specific calculation involving constant B.
return B * (dimensionless_with_units(temp, au)**-0.75) * solar_mass


# TODO(woursler): Migrate to xatu core...
# TODO: Refactor `mass_repr` by migrating its functionality to `xatu.core`.
# quantity_repr(mass, {lunar_mass, earth_mass, jupiter_mass, solar_mass})
# or even quantity_repr(mass, celestial_mass_units)
def mass_repr(mass) -> str:
Expand Down Expand Up @@ -204,9 +229,20 @@ class Protomoon(Planetoid):
protoplanet = attr()


# TODO(woursler): Go over these with a fine tooth comb. Many are not relevant, or only relevant during initialization.
# Many should be properties.
# Many should be initialized differently.
# TODO: Major refactor for `Planet` class attributes:
# 1. Review each attribute currently defaulted to zero or an empty factory
# (e.g., `moon_a`, `core_radius`, `density`, `orb_period`, temperatures, etc.).
# 2. Identify attributes that should be calculated as properties from other
# more fundamental attributes (e.g., `density` from mass and radius).
# 3. Determine which attributes are true state variables requiring storage
# versus derived values that can be computed on-demand.
# 4. For attributes primarily used during the generation/initialization phase,
# assess if they need to be stored on the `Planet` instance long-term or
# if they can be transient.
# 5. Ensure appropriate default values (e.g., `None` or unit-consistent zeros
# like `0*km` or `0*K`) for attributes that are not immediately known or set
# at instantiation, rather than using plain scalar `0` where it might be
# misleading for quantity types.
@attrs(repr=False)
class Planet():
# Orbital details.
Expand All @@ -219,7 +255,6 @@ class Planet():

moons = attr(factory=list)

# ZEROES start here -- TODO(woursler): A bunch of these should be other Zero-like types.
gas_giant = attr(default=False) # TRUE if the planet is a gas giant
# semi-major axis of lunar orbit
moon_a = attr(factory=lambda: 0*au, repr=quantity_formatter(au))
Expand Down Expand Up @@ -268,8 +303,6 @@ class Planet():
atmosphere = attr(default=None)
type = attr(default=PlanetType.UNKNOWN) # Type code

# ZEROES end here

def __repr__(self):
if self.atmosphere is None:
atmosphere_string = "No Atmosphere"
Expand Down