Running multiple instances of the same scenario at the same time.

[Blyron]

I am trying to run the same scenario multiple times for RL purpouses; I thought this was good idea:

Perpendicular-orbit phasing example for Basilisk.
- Sat1: reference circular orbit, a = 7000 km, i = 28 deg, Omega = 0 deg.
- Sat2: perpendicular orbit, i = 90 deg, Omega = 90 deg (so orbital planes are ~orthogonal).
Goal:
- At t=0 they are ~1000 km apart (the script finds Sat2 initial true anomaly to satisfy this).
- At t = 6 hours they should have a close approach (~50 m). The script adjusts a2 to meet that.
- Run the full Basilisk sim for 12 hours and plot the relative distance and orbits.
All comments and variable names are in English.
"""
from Basilisk.simulation import spacecraft, spiceInterface
from Basilisk.utilities import (SimulationBaseClass, macros, orbitalMotion,
                                simIncludeGravBody)
import numpy as np
import os
from Basilisk import __path__
import matplotlib.pyplot as plt
import math

bskPath = __path__[0]

# ---------- user parameters ----------
stepTime = 1.0            # simulation step [s]
runTime = 12 * 3600.0      # full run duration [s]
target_time = 6 * 3600.0   # desired close approach time [s]
target_close_m = 50.0      # target close approach [m]
init_sep_m = 1000e3        # desired initial separation [m]
tolerance_a = 1e-3         # tolerance for semi-major axis solver [m]
max_iter = 50              # max iterations for solver
# -------------------------------------

def post_process(dataRec1, dataRec2, design_info):
    a2, f2, i2, Omega2 = design_info
    t = dataRec1.times() * macros.NANO2SEC
    r1 = np.array(dataRec1.r_BN_N)
    r2 = np.array(dataRec2.r_BN_N)
    relDist = np.linalg.norm(r1 - r2, axis=1)

    # Print closest approach and time
    idx_min = np.argmin(relDist)
    print(f"Closest approach at t = {t[idx_min]/3600.0:.4f} h")
    print(f"Minimum separation = {relDist[idx_min]:.3f} m")

    # Plot relative distance
    plt.figure(figsize=(8,4))
    plt.plot(t/3600.0, relDist/1000.0)
    plt.xlabel("Time [hours]")
    plt.ylabel("Relative distance [km]")
    plt.title("Relative distance between Sat1 and Sat2")
    plt.grid(True)

    # Zoom on region around 6h
    plt.figure(figsize=(8,4))
    plt.plot(t/3600.0, relDist)
    plt.xlim(5.7, 6.3)   # around 6h
    plt.xlabel("Time [hours]")
    plt.ylabel("Relative distance [m]")
    plt.title("Close approach region (around 6h)")
    plt.grid(True)

    # 3D plot of orbits
    fig = plt.figure(figsize=(7,7))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(r1[:,0]/1000.0, r1[:,1]/1000.0, r1[:,2]/1000.0, label="Sat1")
    ax.plot(r2[:,0]/1000.0, r2[:,1]/1000.0, r2[:,2]/1000.0, label="Sat2")
    ax.scatter(0,0,0, color='r', marker='o', label='Earth')
    ax.set_xlabel("x [km]"); ax.set_ylabel("y [km]"); ax.set_zlabel("z [km]")
    ax.set_title("Orbit trajectories (12h)")
    ax.legend()
    plt.show()

def rot_z(theta):
    c = math.cos(theta); s = math.sin(theta)
    return np.array([[c, -s, 0.0],
                     [s,  c, 0.0],
                     [0., 0., 1.0]])

def rot_x(theta):
    c = math.cos(theta); s = math.sin(theta)
    return np.array([[1.0, 0.0, 0.0],
                     [0.0,   c, -s],
                     [0.0,   s,  c]])

def position_in_eci_from_orbital(a, f, i, Omega):
    """
    For circular orbits (e=0) the position in orbital plane is:
      r_orb = [a * cos(f), a * sin(f), 0]
    Rotate by R = Rz(Omega) * Rx(i) * Rz(0) (omega=0 assumed).
    Returns ECI position vector [m].
    """
    r_orb = np.array([a * math.cos(f), a * math.sin(f), 0.0])
    R = rot_z(Omega).dot(rot_x(i))
    r_eci = R.dot(r_orb)
    return r_eci

def mean_motion(mu, a):
    if a <= 6378e3:   # must be larger than Earth radius
        return float('nan')
    return math.sqrt(mu / a**3)

def find_sat2_initial_true_anomaly(a1, f1, a2, i1, Omega1, i2, Omega2, desired_sep, mu):
    """
    Find f2 in [0, 2pi) such that |r1(f1)-r2(f2)| approx desired_sep.
    Simple approach: sample many points and pick the best, then refine with golden-section.
    """
    # coarse grid search
    gridN = 2000
    f_vals = np.linspace(0.0, 2.0 * math.pi, gridN, endpoint=False)
    r1 = position_in_eci_from_orbital(a1, f1, i1, Omega1)
    best_idx = None
    best_err = 1e99
    for idx, f2 in enumerate(f_vals):
        r2 = position_in_eci_from_orbital(a2, f2, i2, Omega2)
        err = abs(np.linalg.norm(r1 - r2) - desired_sep)
        if err < best_err:
            best_err = err
            best_idx = idx
    # initial guess
    f2_guess = f_vals[best_idx]

    # refine using golden-section search on small window around guess
    def objective(f2):
        r2 = position_in_eci_from_orbital(a2, f2, i2, Omega2)
        return abs(np.linalg.norm(r1 - r2) - desired_sep)

    # bracket around guess
    delta = 0.5  # radians window
    left = f2_guess - delta
    right = f2_guess + delta

    # wrap to [0, 2pi]
    def wrap(x):
        return (x + 2*math.pi) % (2*math.pi)

    left = wrap(left); right = wrap(right)
    # golden section, careful with wrapping -> sample inside window evenly
    iters = 60
    for _ in range(iters):
        c1 = left + (right - left) * 0.382
        c2 = left + (right - left) * 0.618
        # ensure c1,c2 in [0,2pi] wrapped
        c1w, c2w = wrap(c1), wrap(c2)
        if objective(c1w) < objective(c2w):
            right = c2
        else:
            left = c1
    f2_opt = wrap((left + right) / 2.0)
    return f2_opt

def distance_at_time_for_a2(a2, f2_0, a1, f1_0, i1, Omega1, i2, Omega2, t, mu):
    """
    Compute distance between satellites at time t assuming circular orbits,
    mean motion propagation (f(t) = f0 + n*t), using a2 provided.
    """
    n1 = mean_motion(mu, a1)
    n2 = mean_motion(mu, a2)
    f1_t = f1_0 + n1 * t
    f2_t = f2_0 + n2 * t
    r1 = position_in_eci_from_orbital(a1, f1_t, i1, Omega1)
    r2 = position_in_eci_from_orbital(a2, f2_t, i2, Omega2)
    return np.linalg.norm(r1 - r2)

def solve_for_a2_target(a1, f1_0, a2_initial_guess, f2_0, i1, Omega1, i2, Omega2, target_time, target_distance, mu):
    """
    Solve for a2 such that distance at target_time equals target_distance.
    Use a secant-like method with bounds.
    """
    # initial bracket guesses:
    # try small offsets around a1
    da = 50.0  # start +/- 50 meters
    a_lo = a1 - 500.0
    a_hi = a1 + 500.0
    # ensure positive and sane
    a_lo = max(1000.0, a_lo)
    a_hi = max(a_lo + 1.0, a_hi)

    def f(a):
        if a <= 6378e3:
            return 1e9 
        return distance_at_time_for_a2(a, f2_0, a1, f1_0, i1, Omega1, i2, Omega2, target_time, mu) - target_distance

    f_lo = f(a_lo)
    f_hi = f(a_hi)

    # expand bracket if signs the same
    expand = 2
    attempts = 0
    while f_lo * f_hi > 0 and attempts < 30:
        a_lo = max(1000.0, a_lo - da * (expand))
        a_hi = a_hi + da * (expand)
        f_lo = f(a_lo)
        f_hi = f(a_hi)
        expand *= 2
        attempts += 1

    # if still same sign, fallback to secant from two starting points near a1
    if f_lo * f_hi > 0:
        a0 = a1 - 100.0
        a1g = a1 + 100.0
        fa0 = f(a0)
        fa1 = f(a1g)
        for _ in range(60):
            if abs(fa1 - fa0) < 1e-12:
                break
            a_new = a1g - fa1 * (a1g - a0) / (fa1 - fa0)
            fa_new = f(a_new)
            if abs(fa_new) < 1e-6:
                return a_new
            a0, fa0 = a1g, fa1
            a1g, fa1 = a_new, fa_new
        return a1g

    # bracketed bisection
    for _ in range(max_iter):
        a_mid = 0.5 * (a_lo + a_hi)
        f_mid = f(a_mid)
        if abs(f_mid) < 1e-6:
            return a_mid
        if f_lo * f_mid <= 0:
            a_hi = a_mid
            f_hi = f_mid
        else:
            a_lo = a_mid
            f_lo = f_mid
        if abs(a_hi - a_lo) < tolerance_a:
            return 0.5 * (a_lo + a_hi)

    return 0.5 * (a_lo + a_hi)

# =========== design stage (analytic) ===========
def design_sat2(a1, f1_0, i1, Omega1, desired_init_sep, target_time, target_close_m, mu):
    """
    Design Sat2 initial true anomaly and semi-major axis to:
      - start ~desired_init_sep away from Sat1 at t=0
      - have distance target_close_m at t = target_time
    Sat2 plane is perpendicular: i2 = 90 deg, Omega2 = Omega1 + 90 deg
    """
    i2 = math.radians(90.0)
    Omega2 = Omega1 + math.radians(90.0)

    # initial guess for a2 near a1
    a2_guess = a1 * 1.0

    # find f2_0 so initial separation is close to desired_init_sep
    f2_0 = find_sat2_initial_true_anomaly(a1, f1_0, a2_guess, i1, Omega1, i2, Omega2, desired_init_sep, mu)

    # solve for a2 that makes distance at target_time ~= target_close_m
    a2_solution = solve_for_a2_target(a1, f1_0, a2_guess, f2_0, i1, Omega1, i2, Omega2, target_time, target_close_m, mu)

    # after finding a2, recompute f2_0 to ensure initial separation still ~desired (a2 changed slightly)
    f2_0 = find_sat2_initial_true_anomaly(a1, f1_0, a2_solution, i1, Omega1, i2, Omega2, desired_init_sep, mu)

    return a2_solution, f2_0, i2, Omega2

# =========== Basilisk run ===========
def run_basilisk(index_process=0, show_plots=True):
    # Create sim base
    scSim = SimulationBaseClass.SimBaseClass()
    for i in range(0, 20):
        simTaskName = f"simTask_{i}"
        simProcessName = f"simProcess_{i}"


        # Create process and task
        dynProcess = scSim.CreateNewProcess(simProcessName)
        simulationTimeStep = macros.sec2nano(stepTime)
        dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep))

        # Create spacecraft
        sc1 = spacecraft.Spacecraft()
        sc1.ModelTag = f"bsk-Sat1_{i}"
        scSim.AddModelToTask(simTaskName, sc1)

        sc2 = spacecraft.Spacecraft()
        sc2.ModelTag = f"bsk-Sat2_{i}"
        scSim.AddModelToTask(simTaskName, sc2)

        # Gravity factory
        gravFactory = simIncludeGravBody.gravBodyFactory()
        earth = gravFactory.createEarth()
        earth.isCentralBody = True
        moon = gravFactory.createMoon()
        gravFactory.addBodiesTo(sc1)
        gravFactory.addBodiesTo(sc2)

        # SPICE
        spiceObject = spiceInterface.SpiceInterface()
        spiceObject.ModelTag = f"SpiceInterfaceData_{i}"
        spiceObject.zeroBase = 'Earth'
        spiceObject.addPlanetNames(["Earth", "Moon"])
        kernelPath = bskPath + '/supportData/EphemerisData/'
        spiceObject.SPICEDataPath = kernelPath
        spiceObject.loadSpiceKernel("naif0012.tls", kernelPath)
        spiceObject.loadSpiceKernel("de430.bsp", kernelPath)
        scSim.AddModelToTask(simTaskName, spiceObject)
        earth.planetBodyInMsg.subscribeTo(spiceObject.planetStateOutMsgs[0])
        moon.planetBodyInMsg.subscribeTo(spiceObject.planetStateOutMsgs[1])
        
        # Sat1 base elements
        mu = earth.mu
        oe1 = orbitalMotion.ClassicElements()
        oe1.a = 7000e3
        oe1.e = 0.0
        oe1.i = 28.0 * macros.D2R
        oe1.Omega = 0.0
        oe1.omega = 0.0
        oe1.f = 0.0

        # analytic design for Sat2
        a1 = oe1.a
        f1_0 = oe1.f
        i1 = oe1.i
        Omega1 = oe1.Omega

        # Set Sat1 initial state using Basilisk's element->rv
        r1_init, v1_init = orbitalMotion.elem2rv(mu, oe1)
        sc1.hub.r_CN_NInit = r1_init
        sc1.hub.v_CN_NInit = v1_init

        # Sat2 elements (circular)
        oe2 = orbitalMotion.ClassicElements()
        oe2.a = 7000e3
        oe2.e = 0.0
        oe2.i = 90.0 * macros.D2R
        oe2.Omega = 90.0 * macros.D2R
        oe2.omega = 0.0
        oe2.f = 180.0 * macros.D2R 

        r2_init, v2_init = orbitalMotion.elem2rv(mu, oe2)
        sc2.hub.r_CN_NInit = r2_init
        sc2.hub.v_CN_NInit = v2_init

        # Loggers
        #dataRec1 = sc1.scStateOutMsg.recorder()
        #dataRec2 = sc2.scStateOutMsg.recorder()
        #scSim.AddModelToTask(simTaskName, dataRec1)
        #scSim.AddModelToTask(simTaskName, dataRec2)

    # Initialize and run
    scSim.TotalSim.resetThreads(4)
    scSim.InitializeSimulation()
    scSim.ConfigureStopTime(macros.sec2nano(runTime))
    scSim.ExecuteSimulation()

        #post_process(dataRec1, dataRec2, (a2, f2_0, i2, Omega2))


if __name__ == "__main__":
    run_basilisk(index_process=0, show_plots=True)

""""""
Perpendicular-orbit phasing example for Basilisk.
- Sat1: reference circular orbit, a = 7000 km, i = 28 deg, Omega = 0 deg.
- Sat2: perpendicular orbit, i = 90 deg, Omega = 90 deg (so orbital planes are ~orthogonal).
Goal:
- At t=0 they are ~1000 km apart (the script finds Sat2 initial true anomaly to satisfy this).
- At t = 6 hours they should have a close approach (~50 m). The script adjusts a2 to meet that.
- Run the full Basilisk sim for 12 hours and plot the relative distance and orbits.
All comments and variable names are in English.
"""
from Basilisk.simulation import spacecraft, spiceInterface
from Basilisk.utilities import (SimulationBaseClass, macros, orbitalMotion,
                                simIncludeGravBody)
import numpy as np
import os
from Basilisk import __path__
import matplotlib.pyplot as plt
import math

bskPath = __path__[0]

# ---------- user parameters ----------
stepTime = 1.0            # simulation step [s]
runTime = 12 * 3600.0      # full run duration [s]
target_time = 6 * 3600.0   # desired close approach time [s]
target_close_m = 50.0      # target close approach [m]
init_sep_m = 1000e3        # desired initial separation [m]
tolerance_a = 1e-3         # tolerance for semi-major axis solver [m]
max_iter = 50              # max iterations for solver
# -------------------------------------

def post_process(dataRec1, dataRec2, design_info):
    a2, f2, i2, Omega2 = design_info
    t = dataRec1.times() * macros.NANO2SEC
    r1 = np.array(dataRec1.r_BN_N)
    r2 = np.array(dataRec2.r_BN_N)
    relDist = np.linalg.norm(r1 - r2, axis=1)

    # Print closest approach and time
    idx_min = np.argmin(relDist)
    print(f"Closest approach at t = {t[idx_min]/3600.0:.4f} h")
    print(f"Minimum separation = {relDist[idx_min]:.3f} m")

    # Plot relative distance
    plt.figure(figsize=(8,4))
    plt.plot(t/3600.0, relDist/1000.0)
    plt.xlabel("Time [hours]")
    plt.ylabel("Relative distance [km]")
    plt.title("Relative distance between Sat1 and Sat2")
    plt.grid(True)

    # Zoom on region around 6h
    plt.figure(figsize=(8,4))
    plt.plot(t/3600.0, relDist)
    plt.xlim(5.7, 6.3)   # around 6h
    plt.xlabel("Time [hours]")
    plt.ylabel("Relative distance [m]")
    plt.title("Close approach region (around 6h)")
    plt.grid(True)

    # 3D plot of orbits
    fig = plt.figure(figsize=(7,7))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(r1[:,0]/1000.0, r1[:,1]/1000.0, r1[:,2]/1000.0, label="Sat1")
    ax.plot(r2[:,0]/1000.0, r2[:,1]/1000.0, r2[:,2]/1000.0, label="Sat2")
    ax.scatter(0,0,0, color='r', marker='o', label='Earth')
    ax.set_xlabel("x [km]"); ax.set_ylabel("y [km]"); ax.set_zlabel("z [km]")
    ax.set_title("Orbit trajectories (12h)")
    ax.legend()
    plt.show()

def rot_z(theta):
    c = math.cos(theta); s = math.sin(theta)
    return np.array([[c, -s, 0.0],
                     [s,  c, 0.0],
                     [0., 0., 1.0]])

def rot_x(theta):
    c = math.cos(theta); s = math.sin(theta)
    return np.array([[1.0, 0.0, 0.0],
                     [0.0,   c, -s],
                     [0.0,   s,  c]])

def position_in_eci_from_orbital(a, f, i, Omega):
    """
    For circular orbits (e=0) the position in orbital plane is:
      r_orb = [a * cos(f), a * sin(f), 0]
    Rotate by R = Rz(Omega) * Rx(i) * Rz(0) (omega=0 assumed).
    Returns ECI position vector [m].
    """
    r_orb = np.array([a * math.cos(f), a * math.sin(f), 0.0])
    R = rot_z(Omega).dot(rot_x(i))
    r_eci = R.dot(r_orb)
    return r_eci

def mean_motion(mu, a):
    if a <= 6378e3:   # must be larger than Earth radius
        return float('nan')
    return math.sqrt(mu / a**3)

def find_sat2_initial_true_anomaly(a1, f1, a2, i1, Omega1, i2, Omega2, desired_sep, mu):
    """
    Find f2 in [0, 2pi) such that |r1(f1)-r2(f2)| approx desired_sep.
    Simple approach: sample many points and pick the best, then refine with golden-section.
    """
    # coarse grid search
    gridN = 2000
    f_vals = np.linspace(0.0, 2.0 * math.pi, gridN, endpoint=False)
    r1 = position_in_eci_from_orbital(a1, f1, i1, Omega1)
    best_idx = None
    best_err = 1e99
    for idx, f2 in enumerate(f_vals):
        r2 = position_in_eci_from_orbital(a2, f2, i2, Omega2)
        err = abs(np.linalg.norm(r1 - r2) - desired_sep)
        if err < best_err:
            best_err = err
            best_idx = idx
    # initial guess
    f2_guess = f_vals[best_idx]

    # refine using golden-section search on small window around guess
    def objective(f2):
        r2 = position_in_eci_from_orbital(a2, f2, i2, Omega2)
        return abs(np.linalg.norm(r1 - r2) - desired_sep)

    # bracket around guess
    delta = 0.5  # radians window
    left = f2_guess - delta
    right = f2_guess + delta

    # wrap to [0, 2pi]
    def wrap(x):
        return (x + 2*math.pi) % (2*math.pi)

    left = wrap(left); right = wrap(right)
    # golden section, careful with wrapping -> sample inside window evenly
    iters = 60
    for _ in range(iters):
        c1 = left + (right - left) * 0.382
        c2 = left + (right - left) * 0.618
        # ensure c1,c2 in [0,2pi] wrapped
        c1w, c2w = wrap(c1), wrap(c2)
        if objective(c1w) < objective(c2w):
            right = c2
        else:
            left = c1
    f2_opt = wrap((left + right) / 2.0)
    return f2_opt

def distance_at_time_for_a2(a2, f2_0, a1, f1_0, i1, Omega1, i2, Omega2, t, mu):
    """
    Compute distance between satellites at time t assuming circular orbits,
    mean motion propagation (f(t) = f0 + n*t), using a2 provided.
    """
    n1 = mean_motion(mu, a1)
    n2 = mean_motion(mu, a2)
    f1_t = f1_0 + n1 * t
    f2_t = f2_0 + n2 * t
    r1 = position_in_eci_from_orbital(a1, f1_t, i1, Omega1)
    r2 = position_in_eci_from_orbital(a2, f2_t, i2, Omega2)
    return np.linalg.norm(r1 - r2)

def solve_for_a2_target(a1, f1_0, a2_initial_guess, f2_0, i1, Omega1, i2, Omega2, target_time, target_distance, mu):
    """
    Solve for a2 such that distance at target_time equals target_distance.
    Use a secant-like method with bounds.
    """
    # initial bracket guesses:
    # try small offsets around a1
    da = 50.0  # start +/- 50 meters
    a_lo = a1 - 500.0
    a_hi = a1 + 500.0
    # ensure positive and sane
    a_lo = max(1000.0, a_lo)
    a_hi = max(a_lo + 1.0, a_hi)

    def f(a):
        if a <= 6378e3:
            return 1e9 
        return distance_at_time_for_a2(a, f2_0, a1, f1_0, i1, Omega1, i2, Omega2, target_time, mu) - target_distance

    f_lo = f(a_lo)
    f_hi = f(a_hi)

    # expand bracket if signs the same
    expand = 2
    attempts = 0
    while f_lo * f_hi > 0 and attempts < 30:
        a_lo = max(1000.0, a_lo - da * (expand))
        a_hi = a_hi + da * (expand)
        f_lo = f(a_lo)
        f_hi = f(a_hi)
        expand *= 2
        attempts += 1

    # if still same sign, fallback to secant from two starting points near a1
    if f_lo * f_hi > 0:
        a0 = a1 - 100.0
        a1g = a1 + 100.0
        fa0 = f(a0)
        fa1 = f(a1g)
        for _ in range(60):
            if abs(fa1 - fa0) < 1e-12:
                break
            a_new = a1g - fa1 * (a1g - a0) / (fa1 - fa0)
            fa_new = f(a_new)
            if abs(fa_new) < 1e-6:
                return a_new
            a0, fa0 = a1g, fa1
            a1g, fa1 = a_new, fa_new
        return a1g

    # bracketed bisection
    for _ in range(max_iter):
        a_mid = 0.5 * (a_lo + a_hi)
        f_mid = f(a_mid)
        if abs(f_mid) < 1e-6:
            return a_mid
        if f_lo * f_mid <= 0:
            a_hi = a_mid
            f_hi = f_mid
        else:
            a_lo = a_mid
            f_lo = f_mid
        if abs(a_hi - a_lo) < tolerance_a:
            return 0.5 * (a_lo + a_hi)

    return 0.5 * (a_lo + a_hi)

# =========== design stage (analytic) ===========
def design_sat2(a1, f1_0, i1, Omega1, desired_init_sep, target_time, target_close_m, mu):
    """
    Design Sat2 initial true anomaly and semi-major axis to:
      - start ~desired_init_sep away from Sat1 at t=0
      - have distance target_close_m at t = target_time
    Sat2 plane is perpendicular: i2 = 90 deg, Omega2 = Omega1 + 90 deg
    """
    i2 = math.radians(90.0)
    Omega2 = Omega1 + math.radians(90.0)

    # initial guess for a2 near a1
    a2_guess = a1 * 1.0

    # find f2_0 so initial separation is close to desired_init_sep
    f2_0 = find_sat2_initial_true_anomaly(a1, f1_0, a2_guess, i1, Omega1, i2, Omega2, desired_init_sep, mu)

    # solve for a2 that makes distance at target_time ~= target_close_m
    a2_solution = solve_for_a2_target(a1, f1_0, a2_guess, f2_0, i1, Omega1, i2, Omega2, target_time, target_close_m, mu)

    # after finding a2, recompute f2_0 to ensure initial separation still ~desired (a2 changed slightly)
    f2_0 = find_sat2_initial_true_anomaly(a1, f1_0, a2_solution, i1, Omega1, i2, Omega2, desired_init_sep, mu)

    return a2_solution, f2_0, i2, Omega2

# =========== Basilisk run ===========
def run_basilisk(index_process=0, show_plots=True):
    # Create sim base
    scSim = SimulationBaseClass.SimBaseClass()
    for i in range(0, 20):
        simTaskName = f"simTask_{i}"
        simProcessName = f"simProcess_{i}"


        # Create process and task
        dynProcess = scSim.CreateNewProcess(simProcessName)
        simulationTimeStep = macros.sec2nano(stepTime)
        dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep))

        # Create spacecraft
        sc1 = spacecraft.Spacecraft()
        sc1.ModelTag = f"bsk-Sat1_{i}"
        scSim.AddModelToTask(simTaskName, sc1)

        sc2 = spacecraft.Spacecraft()
        sc2.ModelTag = f"bsk-Sat2_{i}"
        scSim.AddModelToTask(simTaskName, sc2)

        # Gravity factory
        gravFactory = simIncludeGravBody.gravBodyFactory()
        earth = gravFactory.createEarth()
        earth.isCentralBody = True
        moon = gravFactory.createMoon()
        gravFactory.addBodiesTo(sc1)
        gravFactory.addBodiesTo(sc2)

        # SPICE
        spiceObject = spiceInterface.SpiceInterface()
        spiceObject.ModelTag = f"SpiceInterfaceData_{i}"
        spiceObject.zeroBase = 'Earth'
        spiceObject.addPlanetNames(["Earth", "Moon"])
        kernelPath = bskPath + '/supportData/EphemerisData/'
        spiceObject.SPICEDataPath = kernelPath
        spiceObject.loadSpiceKernel("naif0012.tls", kernelPath)
        spiceObject.loadSpiceKernel("de430.bsp", kernelPath)
        scSim.AddModelToTask(simTaskName, spiceObject)
        earth.planetBodyInMsg.subscribeTo(spiceObject.planetStateOutMsgs[0])
        moon.planetBodyInMsg.subscribeTo(spiceObject.planetStateOutMsgs[1])
        
        # Sat1 base elements
        mu = earth.mu
        oe1 = orbitalMotion.ClassicElements()
        oe1.a = 7000e3
        oe1.e = 0.0
        oe1.i = 28.0 * macros.D2R
        oe1.Omega = 0.0
        oe1.omega = 0.0
        oe1.f = 0.0

        # analytic design for Sat2
        a1 = oe1.a
        f1_0 = oe1.f
        i1 = oe1.i
        Omega1 = oe1.Omega

        # Set Sat1 initial state using Basilisk's element->rv
        r1_init, v1_init = orbitalMotion.elem2rv(mu, oe1)
        sc1.hub.r_CN_NInit = r1_init
        sc1.hub.v_CN_NInit = v1_init

        # Sat2 elements (circular)
        oe2 = orbitalMotion.ClassicElements()
        oe2.a = 7000e3
        oe2.e = 0.0
        oe2.i = 90.0 * macros.D2R
        oe2.Omega = 90.0 * macros.D2R
        oe2.omega = 0.0
        oe2.f = 180.0 * macros.D2R 

        r2_init, v2_init = orbitalMotion.elem2rv(mu, oe2)
        sc2.hub.r_CN_NInit = r2_init
        sc2.hub.v_CN_NInit = v2_init

        # Loggers
        #dataRec1 = sc1.scStateOutMsg.recorder()
        #dataRec2 = sc2.scStateOutMsg.recorder()
        #scSim.AddModelToTask(simTaskName, dataRec1)
        #scSim.AddModelToTask(simTaskName, dataRec2)

    # Initialize and run
    scSim.TotalSim.resetThreads(4)
    scSim.InitializeSimulation()
    scSim.ConfigureStopTime(macros.sec2nano(runTime))
    scSim.ExecuteSimulation()

        #post_process(dataRec1, dataRec2, (a2, f2_0, i2, Omega2))


if __name__ == "__main__":
    run_basilisk(index_process=0, show_plots=True)

Perpendicular-orbit phasing example for Basilisk.
- Sat1: reference circular orbit, a = 7000 km, i = 28 deg, Omega = 0 deg.
- Sat2: perpendicular orbit, i = 90 deg, Omega = 90 deg (so orbital planes are ~orthogonal).
Goal:
- At t=0 they are ~1000 km apart (the script finds Sat2 initial true anomaly to satisfy this).
- At t = 6 hours they should have a close approach (~50 m). The script adjusts a2 to meet that.
- Run the full Basilisk sim for 12 hours and plot the relative distance and orbits.
All comments and variable names are in English.
"""
from Basilisk.simulation import spacecraft, spiceInterface
from Basilisk.utilities import (SimulationBaseClass, macros, orbitalMotion,
                                simIncludeGravBody)
import numpy as np
import os
from Basilisk import __path__
import matplotlib.pyplot as plt
import math

bskPath = __path__[0]

# ---------- user parameters ----------
stepTime = 1.0            # simulation step [s]
runTime = 12 * 3600.0      # full run duration [s]
target_time = 6 * 3600.0   # desired close approach time [s]
target_close_m = 50.0      # target close approach [m]
init_sep_m = 1000e3        # desired initial separation [m]
tolerance_a = 1e-3         # tolerance for semi-major axis solver [m]
max_iter = 50              # max iterations for solver
# -------------------------------------

def post_process(dataRec1, dataRec2, design_info):
    a2, f2, i2, Omega2 = design_info
    t = dataRec1.times() * macros.NANO2SEC
    r1 = np.array(dataRec1.r_BN_N)
    r2 = np.array(dataRec2.r_BN_N)
    relDist = np.linalg.norm(r1 - r2, axis=1)

    # Print closest approach and time
    idx_min = np.argmin(relDist)
    print(f"Closest approach at t = {t[idx_min]/3600.0:.4f} h")
    print(f"Minimum separation = {relDist[idx_min]:.3f} m")

    # Plot relative distance
    plt.figure(figsize=(8,4))
    plt.plot(t/3600.0, relDist/1000.0)
    plt.xlabel("Time [hours]")
    plt.ylabel("Relative distance [km]")
    plt.title("Relative distance between Sat1 and Sat2")
    plt.grid(True)

    # Zoom on region around 6h
    plt.figure(figsize=(8,4))
    plt.plot(t/3600.0, relDist)
    plt.xlim(5.7, 6.3)   # around 6h
    plt.xlabel("Time [hours]")
    plt.ylabel("Relative distance [m]")
    plt.title("Close approach region (around 6h)")
    plt.grid(True)

    # 3D plot of orbits
    fig = plt.figure(figsize=(7,7))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(r1[:,0]/1000.0, r1[:,1]/1000.0, r1[:,2]/1000.0, label="Sat1")
    ax.plot(r2[:,0]/1000.0, r2[:,1]/1000.0, r2[:,2]/1000.0, label="Sat2")
    ax.scatter(0,0,0, color='r', marker='o', label='Earth')
    ax.set_xlabel("x [km]"); ax.set_ylabel("y [km]"); ax.set_zlabel("z [km]")
    ax.set_title("Orbit trajectories (12h)")
    ax.legend()
    plt.show()

def rot_z(theta):
    c = math.cos(theta); s = math.sin(theta)
    return np.array([[c, -s, 0.0],
                     [s,  c, 0.0],
                     [0., 0., 1.0]])

def rot_x(theta):
    c = math.cos(theta); s = math.sin(theta)
    return np.array([[1.0, 0.0, 0.0],
                     [0.0,   c, -s],
                     [0.0,   s,  c]])

def position_in_eci_from_orbital(a, f, i, Omega):
    """
    For circular orbits (e=0) the position in orbital plane is:
      r_orb = [a * cos(f), a * sin(f), 0]
    Rotate by R = Rz(Omega) * Rx(i) * Rz(0) (omega=0 assumed).
    Returns ECI position vector [m].
    """
    r_orb = np.array([a * math.cos(f), a * math.sin(f), 0.0])
    R = rot_z(Omega).dot(rot_x(i))
    r_eci = R.dot(r_orb)
    return r_eci

def mean_motion(mu, a):
    if a <= 6378e3:   # must be larger than Earth radius
        return float('nan')
    return math.sqrt(mu / a**3)

def find_sat2_initial_true_anomaly(a1, f1, a2, i1, Omega1, i2, Omega2, desired_sep, mu):
    """
    Find f2 in [0, 2pi) such that |r1(f1)-r2(f2)| approx desired_sep.
    Simple approach: sample many points and pick the best, then refine with golden-section.
    """
    # coarse grid search
    gridN = 2000
    f_vals = np.linspace(0.0, 2.0 * math.pi, gridN, endpoint=False)
    r1 = position_in_eci_from_orbital(a1, f1, i1, Omega1)
    best_idx = None
    best_err = 1e99
    for idx, f2 in enumerate(f_vals):
        r2 = position_in_eci_from_orbital(a2, f2, i2, Omega2)
        err = abs(np.linalg.norm(r1 - r2) - desired_sep)
        if err < best_err:
            best_err = err
            best_idx = idx
    # initial guess
    f2_guess = f_vals[best_idx]

    # refine using golden-section search on small window around guess
    def objective(f2):
        r2 = position_in_eci_from_orbital(a2, f2, i2, Omega2)
        return abs(np.linalg.norm(r1 - r2) - desired_sep)

    # bracket around guess
    delta = 0.5  # radians window
    left = f2_guess - delta
    right = f2_guess + delta

    # wrap to [0, 2pi]
    def wrap(x):
        return (x + 2*math.pi) % (2*math.pi)

    left = wrap(left); right = wrap(right)
    # golden section, careful with wrapping -> sample inside window evenly
    iters = 60
    for _ in range(iters):
        c1 = left + (right - left) * 0.382
        c2 = left + (right - left) * 0.618
        # ensure c1,c2 in [0,2pi] wrapped
        c1w, c2w = wrap(c1), wrap(c2)
        if objective(c1w) < objective(c2w):
            right = c2
        else:
            left = c1
    f2_opt = wrap((left + right) / 2.0)
    return f2_opt

def distance_at_time_for_a2(a2, f2_0, a1, f1_0, i1, Omega1, i2, Omega2, t, mu):
    """
    Compute distance between satellites at time t assuming circular orbits,
    mean motion propagation (f(t) = f0 + n*t), using a2 provided.
    """
    n1 = mean_motion(mu, a1)
    n2 = mean_motion(mu, a2)
    f1_t = f1_0 + n1 * t
    f2_t = f2_0 + n2 * t
    r1 = position_in_eci_from_orbital(a1, f1_t, i1, Omega1)
    r2 = position_in_eci_from_orbital(a2, f2_t, i2, Omega2)
    return np.linalg.norm(r1 - r2)

def solve_for_a2_target(a1, f1_0, a2_initial_guess, f2_0, i1, Omega1, i2, Omega2, target_time, target_distance, mu):
    """
    Solve for a2 such that distance at target_time equals target_distance.
    Use a secant-like method with bounds.
    """
    # initial bracket guesses:
    # try small offsets around a1
    da = 50.0  # start +/- 50 meters
    a_lo = a1 - 500.0
    a_hi = a1 + 500.0
    # ensure positive and sane
    a_lo = max(1000.0, a_lo)
    a_hi = max(a_lo + 1.0, a_hi)

    def f(a):
        if a <= 6378e3:
            return 1e9 
        return distance_at_time_for_a2(a, f2_0, a1, f1_0, i1, Omega1, i2, Omega2, target_time, mu) - target_distance

    f_lo = f(a_lo)
    f_hi = f(a_hi)

    # expand bracket if signs the same
    expand = 2
    attempts = 0
    while f_lo * f_hi > 0 and attempts < 30:
        a_lo = max(1000.0, a_lo - da * (expand))
        a_hi = a_hi + da * (expand)
        f_lo = f(a_lo)
        f_hi = f(a_hi)
        expand *= 2
        attempts += 1

    # if still same sign, fallback to secant from two starting points near a1
    if f_lo * f_hi > 0:
        a0 = a1 - 100.0
        a1g = a1 + 100.0
        fa0 = f(a0)
        fa1 = f(a1g)
        for _ in range(60):
            if abs(fa1 - fa0) < 1e-12:
                break
            a_new = a1g - fa1 * (a1g - a0) / (fa1 - fa0)
            fa_new = f(a_new)
            if abs(fa_new) < 1e-6:
                return a_new
            a0, fa0 = a1g, fa1
            a1g, fa1 = a_new, fa_new
        return a1g

    # bracketed bisection
    for _ in range(max_iter):
        a_mid = 0.5 * (a_lo + a_hi)
        f_mid = f(a_mid)
        if abs(f_mid) < 1e-6:
            return a_mid
        if f_lo * f_mid <= 0:
            a_hi = a_mid
            f_hi = f_mid
        else:
            a_lo = a_mid
            f_lo = f_mid
        if abs(a_hi - a_lo) < tolerance_a:
            return 0.5 * (a_lo + a_hi)

    return 0.5 * (a_lo + a_hi)

# =========== design stage (analytic) ===========
def design_sat2(a1, f1_0, i1, Omega1, desired_init_sep, target_time, target_close_m, mu):
    """
    Design Sat2 initial true anomaly and semi-major axis to:
      - start ~desired_init_sep away from Sat1 at t=0
      - have distance target_close_m at t = target_time
    Sat2 plane is perpendicular: i2 = 90 deg, Omega2 = Omega1 + 90 deg
    """
    i2 = math.radians(90.0)
    Omega2 = Omega1 + math.radians(90.0)

    # initial guess for a2 near a1
    a2_guess = a1 * 1.0

    # find f2_0 so initial separation is close to desired_init_sep
    f2_0 = find_sat2_initial_true_anomaly(a1, f1_0, a2_guess, i1, Omega1, i2, Omega2, desired_init_sep, mu)

    # solve for a2 that makes distance at target_time ~= target_close_m
    a2_solution = solve_for_a2_target(a1, f1_0, a2_guess, f2_0, i1, Omega1, i2, Omega2, target_time, target_close_m, mu)

    # after finding a2, recompute f2_0 to ensure initial separation still ~desired (a2 changed slightly)
    f2_0 = find_sat2_initial_true_anomaly(a1, f1_0, a2_solution, i1, Omega1, i2, Omega2, desired_init_sep, mu)

    return a2_solution, f2_0, i2, Omega2

# =========== Basilisk run ===========
def run_basilisk(index_process=0, show_plots=True):
    # Create sim base
    scSim = SimulationBaseClass.SimBaseClass()
    for i in range(0, 20):
        simTaskName = f"simTask_{i}"
        simProcessName = f"simProcess_{i}"


        # Create process and task
        dynProcess = scSim.CreateNewProcess(simProcessName)
        simulationTimeStep = macros.sec2nano(stepTime)
        dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep))

        # Create spacecraft
        sc1 = spacecraft.Spacecraft()
        sc1.ModelTag = f"bsk-Sat1_{i}"
        scSim.AddModelToTask(simTaskName, sc1)

        sc2 = spacecraft.Spacecraft()
        sc2.ModelTag = f"bsk-Sat2_{i}"
        scSim.AddModelToTask(simTaskName, sc2)

        # Gravity factory
        gravFactory = simIncludeGravBody.gravBodyFactory()
        earth = gravFactory.createEarth()
        earth.isCentralBody = True
        moon = gravFactory.createMoon()
        gravFactory.addBodiesTo(sc1)
        gravFactory.addBodiesTo(sc2)

        # SPICE
        spiceObject = spiceInterface.SpiceInterface()
        spiceObject.ModelTag = f"SpiceInterfaceData_{i}"
        spiceObject.zeroBase = 'Earth'
        spiceObject.addPlanetNames(["Earth", "Moon"])
        kernelPath = bskPath + '/supportData/EphemerisData/'
        spiceObject.SPICEDataPath = kernelPath
        spiceObject.loadSpiceKernel("naif0012.tls", kernelPath)
        spiceObject.loadSpiceKernel("de430.bsp", kernelPath)
        scSim.AddModelToTask(simTaskName, spiceObject)
        earth.planetBodyInMsg.subscribeTo(spiceObject.planetStateOutMsgs[0])
        moon.planetBodyInMsg.subscribeTo(spiceObject.planetStateOutMsgs[1])
        
        # Sat1 base elements
        mu = earth.mu
        oe1 = orbitalMotion.ClassicElements()
        oe1.a = 7000e3
        oe1.e = 0.0
        oe1.i = 28.0 * macros.D2R
        oe1.Omega = 0.0
        oe1.omega = 0.0
        oe1.f = 0.0

        # analytic design for Sat2
        a1 = oe1.a
        f1_0 = oe1.f
        i1 = oe1.i
        Omega1 = oe1.Omega

        # Set Sat1 initial state using Basilisk's element->rv
        r1_init, v1_init = orbitalMotion.elem2rv(mu, oe1)
        sc1.hub.r_CN_NInit = r1_init
        sc1.hub.v_CN_NInit = v1_init

        # Sat2 elements (circular)
        oe2 = orbitalMotion.ClassicElements()
        oe2.a = 7000e3
        oe2.e = 0.0
        oe2.i = 90.0 * macros.D2R
        oe2.Omega = 90.0 * macros.D2R
        oe2.omega = 0.0
        oe2.f = 180.0 * macros.D2R 

        r2_init, v2_init = orbitalMotion.elem2rv(mu, oe2)
        sc2.hub.r_CN_NInit = r2_init
        sc2.hub.v_CN_NInit = v2_init

        # Loggers
        #dataRec1 = sc1.scStateOutMsg.recorder()
        #dataRec2 = sc2.scStateOutMsg.recorder()
        #scSim.AddModelToTask(simTaskName, dataRec1)
        #scSim.AddModelToTask(simTaskName, dataRec2)

    # Initialize and run
    scSim.TotalSim.resetThreads(4)
    scSim.InitializeSimulation()
    scSim.ConfigureStopTime(macros.sec2nano(runTime))
    scSim.ExecuteSimulation()

        #post_process(dataRec1, dataRec2, (a2, f2_0, i2, Omega2))


if __name__ == "__main__":
    run_basilisk(index_process=0, show_plots=True)

But spiceinterface do not like multithreading and running different simulation tasks. Is there a way to run multiple dynProcess = scSim.CreateNewProcess(simProcessName) simultaneosusly ?

[Mark2000]

I would recommend checking out BSK-RL, which is our package for doing reinforcement learning with Basilisk. We’ve addressed many of the issues that come up when trying to run many simulations in parallel in that repo: https://avslab.github.io/bsk_rl/