resonant-chain/resonance_chain_planetary_system.py at main · amusecode/resonant-chain · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#!/usr/bin/env python


import glob
import numpy as np

from amuse.ext.orbital_elements import orbital_elements
from amuse.io import write_set_to_file, read_set_from_file
from amuse.lab import Particles
from amuse.units import units, constants

from generate_resonant_chain import bring_planet_pair_in_resonance
from generate_resonant_chain import semi_to_orbital_period


def resonant_chain_planetary_system(bodies, tau_a_factor, t_evol, n_steps):
    """
    Iterate over adjacent planets and shake them into resonance.
    Args:
        bodies (Particles): The particles representing the planetary system.
        tau_a_factor (float): Migration parameter in terms of the outer orbital period.
        t_evol (float): Integration time in units of the outer orbital period.
        n_steps (int): Number of steps for the integration.
    Returns:
        Particles: The updated particles after creating the resonant chain.
    """
    for pi in range(len(bodies)-2):
        bodies = resonant_pair_planetary_system(
                        bodies=bodies,
                        inner_planet_id=pi,
                        tau_a_factor=tau_a_factor,
                        t_evol=t_evol,
                        n_steps=n_steps
                        )
    return bodies

def resonant_pair_planetary_system(bodies, inner_planet_id=0, outer_planet_id=1,
                                   tau_a_factor=-1e5, t_evol=100, n_steps=100):
    """
    Create a resonant pair of planets in a planetary system.
    Assumes host to be most massive particle in the system, planets to be all other
    massive particles.
    Args:
        bodies (Particles): The particles representing the planetary system.
        inner_planet_id (int): Index of the inner planet.
        outer_planet_id (int): Index of the outer planet.
        tau_a_factor (float): Migration parameter in terms of the outer orbital period.
        t_evol (float): Integration time in units of the outer orbital period.
        n_steps (int): Number of steps for the integration.
    """
    star = bodies[bodies.mass.argmax()]
    planets = bodies[bodies.type=="PLANET"]
    for pi in planets:
        orbital_elements = orbital_elements(star + pi)
        pi.sma = orbital_elements[2]

    planet_a = planets[inner_planet_id]
    planet_b = planets[outer_planet_id]
    Porb_a = semi_to_orbital_period(planet_a.sma, star.mass + planet_a.mass)[0]
    Porb_b = semi_to_orbital_period(planet_b.sma, star.mass + planet_b.mass)[0]
    if np.isnan(Porb_a.value_in(units.s)) or np.isnan(Porb_b.value_in(units.s)):
        print("Ill-defined orbital periods, returning original bodies")
        return bodies

    bring_planet_pair_in_resonance(
        planetary_system=bodies,
        outer_planet=planet_b,
        tau_a_factor=tau_a_factor,
        t_evol=t_evol,
        n_steps=n_steps
    )
    return bodies

def new_option_parser():
    from amuse.units.optparse import OptionParser
    result = OptionParser()
    result.add_option("--inner", dest="inner_planet_id",
                      type="int", default = 0,
                      help="inner planet id [%default]")
    result.add_option("--outer", dest="outer_planet_id",
                      type="int", default = 1,
                      help="outer planet id [%default]")
    result.add_option("-f", dest="infilename",
                      default = "input_filename.amuse",
                      help="input infilename [%default]")
    result.add_option("-F", dest="outfilename",
                      default = "output_filename.amuse",
                      help="output infilename [%default]")
    result.add_option("--n_steps", dest="n_steps",
                      default = 100, type="int",
                      help="number of steps [%default]")
    result.add_option("--t_evol", dest="t_evol",
                      default = 1000, type="float",
                      help="integration time in units of the outer orbital period [%default]")
    result.add_option("--tau", dest="tau_a_factor",
                      type="float", default = -1e5,
                      help="migration parameter (in terms of outer orbital period) [%default]")
    return result

if __name__ in ('__main__', '__plot__'):
    o, arguments  = new_option_parser().parse_args()

    system = read_set_from_file("planetary_system.hdf5")
    if len(system) < 3:
        print(f"System has 1 planet, skipping")
        exit(-1)

    host = system[system.mass.argmax()]
    system = system.sorted_by_attribute("sma")
    p = system - host


    print(f"Original attributes:")
    print(f"Masses: {system.mass}")
    for pl in p:
        ke = orbital_elements(pl + host, G=constants.G)
        print(f"ecc={ke[3]}, ", end=" ")
        print(f"sma={ke[2].in_(units.au)}", end=" ")
        print(f"inc={ke[5].in_(units.deg)}")

    system = resonant_chain_planetary_system(system, o.tau_a_factor, o.t_evol, o.n_steps)
    host = system[system.mass.argmax()]
    p = system - host
    print(f"After shaking:")
    print(f"Masses: {system.mass}")
    for pl in p:
        ke = orbital_elements(pl + host, G=constants.G)
        print(f"ecc={ke[3]}, ", end=" ")
        print(f"sma={ke[2].in_(units.au)}", end=" ")
        print(f"inc={ke[5].in_(units.deg)}")

    if np.isnan(host.x.value_in(units.pc)):
        print("!!! Task unsuccessful: Resonance not found. !!!")
        exit(-1)

    write_set_to_file(system, o.outfilename, "hdf5", append_to_file=False)