Custom modules

This page contains all code of the custom modules developed for the manuscript.

derive_predictions.py

"""
This module provides functions to optimise stimulation timing and SSC
parameters for periodic MTC length trajectories using a Hill-type MTC model.

Functions
---------
opt_ssc_par(stim, cf, fts, mle, lmtc_avg, muspar, initial_guess={})
    Optimise free SSC parameters using bounded numerical minimisation.
    
opt_stim(x, stim, cf, fts, mle, lmtc_avg, muspar, initial_guess={})
    Optimise stimulation timing to maximise AMPO.
    
sim_periodic(t_stim, cf, fts, mle, lmtc_avg, muspar)
    Simulate a single periodic cycle of MTC dynamics.
"""

import concurrent.futures
import numpy as np
from scipy import optimize

import hillmodel, trajectories

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def opt_ssc_par(stim, cf, fts, mle, lmtc_avg, muspar, initial_guess={}):
    """
    Optimise SSC parameters using bounded numerical minimisation.

    This function optimises whichever SSC parameters are passed as `None` or
    as bounds. Fixed parameters are held constant, while free parameters are
    optimised together with stimulation timing to maximise AMPO.

    Parameters
    ----------
    stim : int
        Type of stimulation passed to `opt_stim` (1 = single pulse,
        2 = stimulation block with onset and offset).
    cf : float or tuple or None
        Cycle frequency [Hz]. If tuple, interpreted as (lower, upper) bounds.
        If float, fixed value. If None, default bounds are used.
    fts : float or tuple or None
        Fraction of cycle time spent shortening. Same rules as `cf`.
    mle : float or tuple or None
        MTC length excursion [m]. Same rules as `cf`.
    lmtc_avg : float
        Average MTC length [m].
    muspar : dict
        Muscle parameter set passed to the Hill-type MTC model.
    initial_guess : dict, optional
        Dictionary of initial guesses. Keys may include:
        - 'cfGuess'
        - 'ftsGuess'
        - 'mleGuess'
        - 'stimGuess'

    Returns
    -------
    p_mech : float
        Optimised AMPO [W].
    y : tuple
        Simulation output returned by `sim_periodic`.
    x_opt : numpy.ndarray
        Optimised SSC parameter vector.

    Notes
    -----
    Uses Nelder-Mead optimisation with bounds passed to SciPy.

    """

    bounds = []
    x0 = []

    # CF
    if isinstance(cf, tuple):
        cf_bounds = cf
    elif cf is None:
        cf_bounds = (0.2, 12)
    else:
        cf_bounds = (cf, cf)

    if not isinstance(cf, (int, float)):
        bounds.append(cf_bounds)
        cf_guess = initial_guess.get(
            "cfGuess",
            np.random.rand() * 3 + 1  # 1–4 Hz default guess
        )
        x0.append(cf_guess)

    # FTS
    if isinstance(fts, tuple):
        fts_bounds = fts
    elif fts is None:
        fts_bounds = (0.02, 0.98)
    else:
        fts_bounds = (fts, fts)

    if not isinstance(fts, (int, float)):
        bounds.append(fts_bounds)
        fts_guess = initial_guess.get(
            "ftsGuess",
            np.random.uniform(*fts_bounds)
        )
        x0.append(fts_guess)

    # MLE
    if isinstance(mle, tuple):
        mle_bounds = mle
    elif mle is None:
        mle_bounds = (0.5e-3, 16e-3)
    else:
        mle_bounds = (mle, mle)

    if not isinstance(mle, (int, float)):
        bounds.append(mle_bounds)
        mle_guess = initial_guess.get(
            "mleGuess",
            np.random.uniform(*mle_bounds)
        )
        x0.append(mle_guess)
    
    # Objective function (maximize stimOpt -> minimize negative)
    def objective(x):
        return -opt_stim(x, stim, cf, fts, mle, lmtc_avg, muspar, initial_guess)[0]

    result = optimize.minimize(
        objective,
        x0,
        method="Nelder-Mead",
        bounds=bounds,
        options={"xatol": 1e-6, "fatol": 1e-6}
    )

    x_opt = result.x

    p_mech, y = opt_stim(
        x_opt, stim, cf, fts, mle, lmtc_avg, muspar, initial_guess
    )

    return p_mech, y, x_opt

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def opt_stim(x, stim, cf, fts, mle, lmtc_avg, muspar, initial_guess={}):
    """
    Optimise stimulation onset and offset to maximise AMPO.

    Parameters
    ----------
    x : array-like
        Current values of the free SSC parameters being optimised.
    stim : int
        Type of stimulation (1 = single pulse, 2 = constant). Determines bounds.
    cf : float
        Cycle frequency [Hz].
    fts : float
        Fraction of cycle time spent shortening [-].
    mle : float
        MTC length excursion [m].
    lmtc_avg : float
        Average MTC length [m].
    muspar : dict
        Muscle parameters dictionary including gamma parameters.
    
    Returns
    -------
    p_mech : float
        AMPO [W].
    sim_result : tuple
        Full simulation results from `sim_periodic`.
    """
    
    # Print which variables we have
    print('Initial guess SSCpar = ..')
    print(x)
    print('Initial guess Stim = ..')
    print(initial_guess['stimGuess'])
    
    # Unpack SSC parameters from 'x'
    if type(cf) == tuple or cf == None:
        cf = x[0]
        x = x[1:]
        # print('Optimising CF')
    if type(fts) == tuple or fts == None: 
        fts = x[0]
        x = x[1:]
        # print('Optimising FTS')
    if type(mle) == tuple or mle == None: # we are imposing fts
        mle = x[0]
        x = x[1:]
        # print('Optimising MLE')
    
    # Make intial guess for stim
    stim_guess = initial_guess.get('stimGuess', np.nan)
    
    # Define bounds and initial guess
    if stim == 2:
        bounds = ((-np.inf, np.inf), (0, np.inf))
        if np.isnan(stim_guess).any():
            x0 = [0, fts / cf * 0.2]
        else:
            x0 = stim_guess
    elif stim == 1:
        bounds = ((0, fts / cf),)
        if np.isnan(stim_guess).any():
            x0 = [fts / cf * 0.2]
        else:
            x0 = stim_guess
            
    # Objective: negative mechanical power
    objective = lambda x: -sim_periodic(x, cf, fts, mle, lmtc_avg, muspar)[0]

    result = optimize.minimize(
        objective,
        x0,
        method='Nelder-Mead',
        bounds=bounds,
        options={'xatol': 1e-5, 'fatol': 1e-3}
    )

    t_stim_opt = result.x
    p_mech, sim_result = sim_periodic(t_stim_opt, cf, fts, mle, lmtc_avg, muspar)
    initial_guess['stimGuess'] = t_stim_opt # update initialGuess based stimOpt of previous round -> faster convergence at the end!
    
    return p_mech, sim_result

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def sim_periodic(t_stim, cf, fts, mle, lmtc_avg, muspar):
    """
    Simulate a single periodic cycle of muscle-tendon complex (MTC) dynamics.

    Parameters
    ----------
    t_stim : array-like or float
        Stimulation onset and offset times. If single value, onset = 0.
    cf : float
        Cycle frequency [Hz].
    fts : float
        Fraction of cycle time spent shortening [-].
    mle : float
        MTC length excursion [m].
    lmtc_avg : float
        Average MTC length [m].
    muspar : dict
        Muscle parameters dictionary including gamma parameters.
    
    Returns
    -------
    p_mech : float
        AMPO [W].
    sim_result : tuple
        Detailed simulation results.
    """
    
    # Determine stimulation times
    if len(t_stim) > 1:
        t_stim_on, t_stim_off = t_stim
    else:
        t_stim_on, t_stim_off = 0, t_stim[0]
    
    # Time discretization
    n_points = 2000
    time = np.unique(np.hstack((
        np.arange(0, fts / cf, 1 / n_points),
        np.arange(fts / cf, 1 / cf, 1 / n_points),
        [1 / cf]
    )))
    
    # MTC length over time
    lmtc = trajectories.cv(time, cf, fts, mle, lmtc_avg)[0]
    
    # Initial states (gamma and lcerel) at t=0
    gamma0 = hillmodel.anly_gamma(0, cf, t_stim_on, t_stim_off, 1, muspar)[0]
    lcerel0 = min(1.4, hillmodel.force_eq(lmtc[0], gamma0, muspar)[1] - 1e-2)
    lcerel_f = [lcerel0]
    
    # Setup solution dictionary
    t_on = [t_stim_on, t_stim_on+1/cf, t_stim_on+2/cf]
    t_off = [t_stim_off, t_stim_off+1/cf, t_stim_off+2/cf]
    
    inputs = {
        'time': time,
        'lmtc': lmtc,
        't_stim': np.array([t_on,t_off]).T,
        'cf': cf
    }

    # Convergence parameters
    dFsee, dLcerel = 1000, 1
    iRound, iFail = 0, 0
    timeout = 10  # seconds

    # Simulate until difference in SEE force is <10 mN
    dFsee, dLcerel = 1000, 1
    iRound, iFail = 0, 0
    timeout = 10 # [s]

    # ODE solver wrapper
    def solve_ode(gamma0, lcerel0, ode_opts, u):
        W_mech, sim_result = hillmodel.solve_simu_mtc(gamma0, lcerel0, muspar, u, ode_opts)
        
        # Extract key variables
        time, _, _, _, lcerel, _, _, _, _, fsee, *_ = sim_result
        
        # Check if solution is complete and within bounds
        if time[-1] != ode_opts['t_eval'][-1] or lcerel[-1] > 2:
            raise RuntimeError("Incomplete simulation or lcerel blew up")
        
        return W_mech, sim_result
    
    ode_opts = {}   
    ode_opts['method'] = 'Radau'
    ode_opts['rtol'] = 1e-9
    ode_opts['atol'] = 1e-6
    ode_opts['t_eval'] = time
    lcerel_f = []
    
    # DEBUGGING REMOVE
    hillmodel.solve_simu_mtc(gamma0, lcerel0, muspar, inputs, ode_opts)
    
    # Sometimes a simulation does get stuck, so if it takes longer than 10s we abort it and try again with a almost identical initial state.
    while dFsee > muspar['fmax']*0.1/100 or abs(dLcerel) > 1e-3:
        with concurrent.futures.ThreadPoolExecutor(max_workers=1) as executor:
            # Using ThreadPoolExecutor to run solve_ivp with timeout
            future = executor.submit(solve_ode,gamma0,lcerel0,ode_opts,inputs) 
            try:
                # Wait for the result with a timeout
                Wmech,y = future.result(timeout=timeout)
                time, lmtc, stim, gamma, lcerel, q, lsee, lpee, fisomrel, fsee, fpee, fce, fcerel, vcerel = y
                dFsee = np.abs(fsee[0]-fsee[-1])
                dLcerel = lcerel[-1]-lcerel[0]
                lcerel_f.append(lcerel[-1])
                # lcerel0 = np.mean(lcerel_f[-3:])  # smooth over last 3 values
                
                # sometimes we have large outliers, select the one within 2 std.
                lcerel_sel = np.array(lcerel_f[-4:]) # select last 4 values
                lcerel_sel = lcerel_sel[np.abs(lcerel_sel - np.mean(lcerel_sel)) <= 2*np.std(lcerel_sel)]
                lcerel0 = np.mean(lcerel_sel)  # avg.
            except:
                print(f"Timeout at iRound {iRound}, trying again")
                lcerel0 -= 0.1
                iFail += 1
            finally:
                # Use shutdown(wait=False) to avoid blocking while cleaning up
                executor.shutdown(wait=False)
  
            iRound += 1
            if iRound > 20 or iFail > 3:
                dFsee, dLcerel, Wmech, y = 0,0,np.nan,None
    
    Pmech = Wmech*cf
    print(f'AMPO = {Pmech*1e3:0.3f} mW')
    # import matplotlib.pyplot as plt
    # plt.figure(); plt.plot(time,lmtc)
    # plt.figure(); plt.plot(time,stim)
    # breakpoint()
    return Pmech, y

helpers.py

"""
This module provides helper functions for loading simulation results and
computing average mechanical power output (AMPO).

Functions
---------
load_sims(cf_set, fts_set, mle_set, mus, data_dir)
    Load simulation files for a grid of SSC parameters and compute AMPO.

get_ampo(filepaths)
    Compute AMPO from one or multiple simulation CSV files.
"""

import os
import numpy as np
import pandas as pd

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def load_sims(cf_set, fts_set, mle_set, mus, data_dir):
    """
    Load simulation data and compute AMPO values for given SSC parameter grid.

    This function iterates over combinations of cycle frequency, FTS and
    MTC length excursion, loads the corresponding CSV files, computes AMPO and
    returns an array with one value for each parameter combination.

    Parameters
    ----------
    cf_set : array-like
        Cycle frequency values [Hz].
    fts_set : array-like
        Fraction of cycle time spent shortening [-].
    mle_set : array-like
        MTC length excursion values [m].
    mus : str
        Base filename prefix used for constructing file names.
    data_dir : str
        Directory path containing the CSV simulation files.

    Returns
    -------
    np.ndarray
        Array of computed AMPO values with shape:
        (len(cf_set), len(fts_set), len(mle_set)),
        or a squeezed version if dimensions are singleton.

    Notes
    -----
    If a file cannot be read or processed, the corresponding entry
    is set to NaN.

    Examples
    --------
    >>> load_sims([2.0], [0.5, 0.75], [0.004], "GMe1", "./data/")
    array([...])
    """
    
    cf_set = np.atleast_1d(cf_set)
    fts_set = np.atleast_1d(fts_set)
    mle_set = np.atleast_1d(mle_set)

    ampo_set = np.full(
        (len(cf_set), len(fts_set), len(mle_set)),
        np.nan
    )

    for i_cf, cf in enumerate(cf_set):
        for i_fts, fts in enumerate(fts_set):
            for i_mle, mle in enumerate(mle_set):
                try:
                    file_name = (
                        f"{mus}_cf{cf:0.1f}Hz_fts{fts:0.2f}_mle{mle*1e3:0.1f}mm"
                    )
                    filepath = os.path.join(data_dir, file_name + ".csv")
                    df = pd.read_csv(filepath)
                    data = df.to_numpy()

                    time, lmtc, _, fsee = data.T[:4]

                    w_mech = -np.trapezoid(fsee, lmtc)
                    ampo_set[i_cf, i_fts, i_mle] = w_mech * cf

                except Exception:
                    ampo_set[i_cf, i_fts, i_mle] = np.nan

    return np.squeeze(ampo_set)

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def get_ampo(filepaths):
    """
    Compute AMPO from one or more CSV files.

    This function loads time-series data from CSV file(s) and computes AMPO as:

        AMPO = -∫ fsee d(lmtc) / time_end

    The function supports flexible input structures, including a single file,
    1D, 2D, or 3D arrays of file paths. Each file is expected to contain at
    least the following rows after transposition:
    time, lmtc, stim, fsee.

    Parameters
    ----------
    filepaths : str or array-like of str
        Path(s) to CSV file(s). Can be:
        - str: single file path
        - 1D list/array: multiple files
        - 2D list/array: grid of file paths
        - 3D list/array: 3D block of file paths

    Returns
    -------
    float or numpy.ndarray
        Computed AMPO values:
        - float for a single file
        - 1D array for 1D input
        - 2D array for 2D input
        - 3D array for 3D input
        Entries are np.nan when computation fails for a file.

    Notes
    -----
    - CSV files are read using pandas and transposed (`.T.to_numpy()`).
    - The function assumes consistent row ordering across files.
    - Integration is performed using the trapezoidal rule.
    - Division by `time[-1]` normalizes the metric by total duration.
    - Any file that cannot be read or processed is assigned np.nan.

    Raises
    ------
    ValueError
        If the input structure has more than 3 dimensions or is unsupported.

    Examples
    --------
    >>> get_ampo("file.csv")
    0.42

    >>> get_ampo(["f1.csv", "f2.csv"])
    array([0.42, nan])

    >>> get_ampo([[ "a.csv", "b.csv" ],
    ...           [ "c.csv", "d.csv" ]])
    array([[0.41, 0.38],
           [0.44, nan]])
    """

    # Case 0: Single string input
    if isinstance(filepaths, str):
        try:
            data = pd.read_csv(filepaths).T.to_numpy()
            time, lmtc, stim, fsee, *_ = data
            return -np.trapz(fsee, lmtc) / time[-1]
        except Exception:
            return np.nan

    # Convert to NumPy array to check dimensions
    filepaths_arr = np.array(filepaths, dtype=object)

    # Case 12: 1D list of filepaths
    if filepaths_arr.ndim == 1:
        AMPOs = []
        for filepath in filepaths_arr:
            try:
                data = pd.read_csv(filepath).T.to_numpy()
                time, lmtc, stim, fsee, *_ = data
                AMPO = -np.trapz(fsee, lmtc) / time[-1]
                AMPOs.append(AMPO)
            except Exception:
                AMPOs.append(np.nan)
        return np.array(AMPOs)

    # Case 2: 2D list of filepaths
    elif filepaths_arr.ndim == 2:
        shape = filepaths_arr.shape
        AMPOs = np.full(shape, np.nan, dtype=float)
        for i in range(shape[0]):
            for j in range(shape[1]):
                try:
                    data = pd.read_csv(filepaths_arr[i, j]).T.to_numpy()
                    time, lmtc, stim, fsee, *_ = data
                    AMPO = -np.trapz(fsee, lmtc) / time[-1]
                    AMPOs[i, j] = AMPO
                except Exception:
                    continue  # Already NaN
        return AMPOs
    
    # Case 3: 3D list of filepaths
    elif filepaths_arr.ndim == 3:
        shape = filepaths_arr.shape
        AMPOs = np.full(shape, np.nan, dtype=float)
        for i in range(shape[0]):
            for j in range(shape[1]):
                for k in range(shape[2]):
                    try:
                        data = pd.read_csv(filepaths_arr[i, j, k]).T.to_numpy()
                        time, lmtc, stim, fsee, *_ = data
                        AMPO = -np.trapz(fsee, lmtc) / time[-1]
                        AMPOs[i, j, k] = AMPO
                    except Exception:
                        continue  # Already NaN
        return AMPOs
    else:
        raise ValueError("Unsupported input structure for 'filepaths'.")

hillmodel.py

"""
This module provides functions for a Hill-type muscle-tendon complex model.

The module includes activation dynamics, force-length and force-velocity
relations, elastic element properties and forward simulation utilities.

This model is extensively described in the following papers:
    -   van Soest, A.J.K. & Bobbert (1993) 
        The contribution of muscle properties in the control of explosive 
        movements
        https://doi.org/10.1007/BF00198959
    -   Reuvers, E.D.H.M. & Kistemaker, D.A. (2025)
        Accuracy of experimentally estimated muscle properties: Evaluation and 
        improvement using a newly developed toolbox
        https://doi.org/10.1101/2025.09.29.678508 

Functions
---------
simu_mtc(t, state, muspar, inputs)
    Compute state derivatives for the Hill-type MTC model.

solve_simu_mtc(gamma0, lcerel0, muspar, inputs, ode_opts=None)
    Forward simulate the Hill-type MTC model.

gamma_dot(gamma, stim, muspar)
    Compute the time derivative of gamma.

act_state(gamma, lcerel, muspar)
    Compute active state from gamma and relative CE length.

force_length(lcerel, muspar)
    Compute the CE force-length relation.

lee2force(lsee, lpee, muspar)
    Compute SEE and PEE forces from elastic element lengths.

force2lee(fsee, fpee, muspar)
    Compute SEE and PEE lengths from elastic element forces.

fce2vce(fce, q, lcerel, muspar)
    Compute CE velocity from CE force.

vce2fce(vce, q, lcerel, muspar)
    Compute CE force from CE velocity.

force_eq(lmtc, gamma, muspar)
    Find the relative CE length satisfying static force equilibrium.

anly_gamma(time, cf, t_stim_on, t_stim_off, stim, muspar)
    Compute periodic gamma analytically for block stimulation.

gamma_relax(time, gamma0, muspar)
    Compute gamma during relaxation.

gamma_rise(time, gamma0, stim, muspar)
    Compute gamma during activation.
"""

import numpy as np
from scipy.integrate import solve_ivp, trapezoid

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def simu_mtc(t, state, muspar, inputs):
    """
    Compute state derivatives for a Hill-type muscle-tendon complex (MTC) model.

    Parameters
    ----------
    t : float
        Current time [s].
    state : array-like, shape (2,)
        Current state vector:
        - state[0] : gamma, normalized Ca2+ concentration between filaments.
        - state[1] : lcerel, relative contractile element (CE) length
                     (lce / lce_opt).
    muspar : dict
        Muscle parameters including gamma dynamics, CE properties, etc.
    inputs : dict
        Input containing:
        - 'time': time axis [s]
        - 'lmtc': MTC length over time [m]
        - 'stim' or 't_stim': stimulation over time or onset/offset times

    Returns
    -------
    gammad : float
        Time derivative of gamma [1/s].
    vcerel : float
        Time derivative of relative CE length [1/s].
    y : list
        List of additional variables for debugging/analysis:
        [lmtc, stim, q, lsee, lpee, fisomrel, fsee, fpee, fce,
        fcerel, vcerel]
    """
    
    gamma, lcerel = state

    # Muscle-tendon length
    lmtc = np.interp(t, inputs['time'], inputs['lmtc'])
    
    # Determine stimulation
    try: # first try if stim onset- and offset times are present
        tStim = inputs['t_stim'] # [s]
        tStim = np.atleast_2d(tStim)
        stim = np.zeros_like(t)
        for start, end in tStim:
            stim[(t >= start) & (t <= end)] = 1
    except KeyError: # if stim onset- and offset time are not present stim(t) should be present
        stim = np.interp(t,inputs['time'],inputs['stim']) # [ ]
       
    # Activation dynamics
    gamma_0 = muspar['gamma_0']
    gamma = (gamma>gamma_0)*gamma + (gamma<=gamma_0)*gamma_0 # [ ]
    q = act_state(gamma, lcerel, muspar)[0]
    gammad = np.where(
        stim >= gamma,
        (stim * (1 - gamma_0) - gamma + gamma_0) / muspar['tact'],
        (stim * (1 - gamma_0) - gamma + gamma_0) / muspar['tdeact']
    )

    # Contraction dynamics
    lce = lcerel * muspar['lce_opt']
    lsee = lmtc - lce
    lpee = lce
    fisomrel = force_length(lcerel, muspar)[0]
    fsee, fpee = lee2force(lsee, lpee, muspar)[0:2]
    fce = fsee - fpee
    fcerel = fce / muspar['fmax']
    vcerel = fce2vce(fce, q, lcerel, muspar)[1]

    # Debugging
    # if np.isnan(vcerel).any():
    #     breakpoint()

    y = [lmtc, stim, q, lsee, lpee, fisomrel, fsee,
         fpee, fce, fcerel, vcerel]

    return gammad, vcerel, y

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def solve_simu_mtc(gamma0, lcerel0, muspar, inputs, ode_opts=None):
    """
    Forward simulation of a Hill-type muscle-tendon complex (MTC) model.

    Parameters
    ----------
    gamma0 : float
        Initial normalized Ca2+ concentration between filaments [ ].
    lcerel0 : float
        Initial relative CE length (lce / lce_opt) [ ].
    muspar : dict
        Muscle parameters dictionary.
    inputs : dict
        Inputs containing time axis, lmtc(t), stim(t) or t_stim.
    ode_opts : dict, optional
        ODE solver options:
        - 'method' : str, solver method (default: 'Radau')
        - 'max_step' : float, maximum step size
        - 'rtol' : float, relative tolerance
        - 'atol' : float, absolute tolerance
        - 't_eval' : array-like, time points to evaluate solution

    Returns
    -------
    W_mech : float
        Mechanical work done by the SEE [J].
    y : list
        List of simulation outputs: [time, lmtc, stim, gamma, lcerel, q,
        lsee, lpee, fisomrel, fsee, fpee, fce, fcerel, vcerel].
    """
    if ode_opts is None:
        ode_opts = {}
    
    # Solver options
    method = ode_opts.get('method', 'Radau')
    max_step = ode_opts.get('max_step', np.inf)
    rtol = ode_opts.get('rtol', 1e-3)
    atol = ode_opts.get('atol', 1e-6)
    t_eval = ode_opts.get('t_eval', inputs['time'])

    # Initial state and timespan
    state0 = [gamma0, lcerel0]
    t_span = [inputs['time'][0], inputs['time'][-1]]

    # ODE function wrapper
    def ode_fun(t, state):
        gammad, vcerel, _ = simu_mtc(t, state, muspar, inputs)
        return [gammad, vcerel]

    # Solve ODE
    sol = solve_ivp(
        ode_fun, t_span, state0,
        method=method, max_step=max_step,
        rtol=rtol, atol=atol, t_eval=t_eval
    )

    if not sol.success:
        W_mech = np.nan
        # breakpoint()
        # raise RuntimeError(f"ODE solver failed: {sol.message}")

    # Evaluate solution at all time points
    gammad, vcerel, y_list = simu_mtc(sol.t, sol.y, muspar, inputs)
    
    lmtc, stim, q, lsee, lpee, fisomrel, fsee, fpee, fce, fcerel, vcerel = y_list[0:11]
    
    y = [sol.t, lmtc, stim, sol.y[0], sol.y[1], q, lsee, lpee,
         fisomrel, fsee, fpee, fce, fcerel, vcerel]

    # Mechanical work by SEE
    W_mech = -trapezoid(fsee, lmtc)
    return W_mech, y

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def gamma_dot(gamma, stim, muspar):
    """
    Compute rate of change of Ca²⁺ concentration.
    
    Parameters
    ----------
    gamma : float or numpy.ndarray
        Current Ca²⁺ level [-].
    stim : float or numpy.ndarray
        Neural stimulation [-].
    muspar : dict
        Parameters: gamma_0, tact, tdeact
    
    Returns
    -------
    gamma_dot : float or numpy.ndarray
        Time derivative of gamma [1/s].
    """
    
    # Unravel parameter values
    gamma_0 = muspar['gamma_0']
    tact    = muspar['tact']
    tdeact  = muspar['tdeact']
    
    # Computations
    gammadot = (stim>=gamma)*((stim*(1-gamma_0)-gamma + gamma_0)/tact) + (stim<gamma)*((stim*(1-gamma_0)-gamma + gamma_0)/tdeact) # [1/s]
    
    # Output
    return gammadot

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def act_state(gamma, lcerel, muspar):
    """
    Computes the active state based on the relative amount of Ca2+ 
    between the filaments and the relative CE length. This version is based on 
    Hatze (1981), p. 37-41, but slightly modified such that the parameter
    values has physiological relevance and that it is easier to use in 
    Optimal Control.
        
    Parameters
    ----------
    gamma : float or numpy.ndarray
        Relative amount Ca2+ between the filaments [-]
    lcerel : float or numpy.ndarray
        Relative CE length (Lce / Lce_opt) [-].
    muspar : dict
        Muscle parameters including q0, kCa, a_act and b_act.
    
    Returns
    -------
    q : float or numpy.ndarray
        Active state (relative amount of Ca2+ bound to troponin C) [-].
    dqdlcerel : float or numpy.ndarray
        Partial derivative of q with respect to lcerel [-].
    gamma_05 : float or numpy.ndarray
        Gamma value at which q = 0.5 [-].
    """
    
    # Unravel parameters values
    q0 = muspar['q0']
    kCa = muspar['kCa']
    a_act = muspar['a_act']
    b_act = muspar['b_act']
    
    # Computations
    a1_act = np.log10(np.exp(a_act))
    B_act = b_act[0] + b_act[1] * lcerel + b_act[2] * lcerel**2
    
    q = q0 + (1 - q0) / (1 + (kCa * gamma)**a1_act * np.exp(a_act * B_act))
    
    try:
        dqdlcerel = (a_act*np.exp(a_act*(B_act))*(gamma*kCa)**(np.log(np.exp(a_act))/np.log(10))*(q0-1)*(b_act[1] + 2*b_act[2]*lcerel))/(np.exp(a_act*(B_act))*(gamma*kCa)**(np.log(np.exp(a_act))/np.log(10)) + 1)**2
        gamma_05 = ((1-0.5)/(kCa**a1_act*np.exp(a_act*B_act)*(0.5-q0)))**(1/a1_act)
    except:
        dqdlcerel = np.nan*gamma
        gamma_05 = np.nan*gamma
    
    return q, dqdlcerel, gamma_05

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def force_length(lcerel, muspar):
    """
    Computes the relative isometric CE force based on the relative CE length.
    
    Parameters
    ----------
    lcerel : float or numpy.ndarray
        Relative CE length (Lce / Lceopt) [-].
    muspar : dict
        Muscle parameters.
    
    Returns
    -------
    fisomrel : float or numpy.ndarray
        Relative isometric CE force (CE isometric force / Fcemax) [-].
    kce : float or numpy.ndarray
        Derivative fisomrel with respect to lcerel [-].
    """
    
    # Unravel parameter values
    n = muspar['n'] # [ ]
    C = -1/muspar['w']**n # [ ]
    
    # Compute tails of parabola (exponential tails)
    Fp = 0.1 # exp function kicks in a 10% fisomrel
    xp = ((Fp-1)/C)**(1/2) # intercept of function with y=Fp
    xp = np.array([-xp+1, xp+1]) # two solutions because of root
    dFdx = 2*C*(xp-1) # first derivative of F at xp
    # function has the form: y=a*exp(b*x), so:
    b = dFdx/Fp
    a = Fp/(np.exp(b*xp))
    
    # Compute fisomrel
    fisomrel = np.clip(
        np.piecewise(
            lcerel, 
            [lcerel < xp[0], (lcerel >= xp[0]) & (lcerel <= xp[1]), lcerel > xp[1]], 
            [lambda x: a[0]*np.exp(b[0]*x),           # Exponential tail for lcerel < xp[0]
             lambda x: C*(x-1)**n + 1,                # Middle section
             lambda x: a[1]*np.exp(b[1]*x)]           # Exponential tail for lcerel > xp[1]
        ),
        1e-9, # Minimum value
        None  # No maximum value
    )
    
    # Compute derivative (i.e., fisomrel/dlcerel)
    kce = np.piecewise(lcerel, 
                        [lcerel < xp[0], (lcerel >= xp[0]) & (lcerel <= xp[1]), lcerel > xp[1]], 
                        [lambda x: a[0]*b[0]*np.exp(b[0]*x), # [ ] dfisomrel/dlcerel for lcerel < xp[0]
                         lambda x: 2*C*(x-1), # [ ] dfisomrel/dlcerel
                         lambda x: a[1]*b[1]*np.exp(b[1]*x)]) # [ ] dfisomrel/dlcerel for lcerel < xp[0]
    
    return fisomrel,kce

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def lee2force(lsee, lpee, muspar):
    """
    Compute forces in SEE and PEE.
    
    Parameters
    ----------
    lsee : float or numpy.ndarray
        SEE length [m]
    lpee : float or numpy.ndarray
        PEE length [m].
    muspar : dict
        Muscle parameters
    
    Returns
    -------
    fsee : float or numpy.ndarray
        SEE force [N]
    fpee : float or numpy.ndarray
        PEE force [N]
    fce : float or numpy.ndarray
        CE force [N]
    fcerel : float or numpy.ndarray
        CE force normalised to maximal isometric CE force [-]
    Ksee : float or numpy.ndarray
        dFsee/dLsee [N/m]
    Kpee : float or numpy.ndarray
        dFPee/dLpee [N/m]
    """
    
    # Unravel parameter values
    lpee0   = muspar['lpee0']   # [m]       PEE slack length
    kpee    = muspar['kpee']    # [N/m^2]   shape parameter that scales the stiffness of PEE
    lsee0   = muspar['lsee0']   # [m]       SEE slack length
    ksee    = muspar['ksee']    # [N/m^2]   shape parameter that scales the stiffness of SEE
    fmax    = muspar['fmax']    # [N]       maximal isometric CE force
    
    # SEE
    esee    = lsee-lsee0 # [m] SEE elongation
    fsee    = (esee<0)*0 + (esee>=0)*(ksee*esee**2) # [N] SEE force
    Ksee    = (esee<0)*0 + (esee>=0)*(2*ksee*esee) # [N/m] dFsee/dLsee
    
    # PEE
    epee    = lpee-lpee0 # PEE elongation
    fpee    = (epee<0)*0 + (epee>=0)*(kpee*epee**2) # [N] PEE force
    Kpee    = (epee<0)*0 + (epee>=0)*(2*kpee*epee) # [N/m] dFpee/dLpee
    
    # CE
    fce     = fsee-fpee # [N]
    fcerel  = fce/fmax # [ ]
    
    # Output
    return fsee,fpee,fce,fcerel,Ksee,Kpee

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def force2lee(fsee,fpee,muspar):
    """
        Compute SEE and PEE length from SEE and PEE force.
    
    Parameters
    ----------
    fsee : float or numpy.ndarray
        SEE force [N]
    fpee : float or numpy.ndarray
        PEE force [N]
        
    muspar : dict
        Muscle parameters.

    Returns
    -------
    lsee : float or numpy.ndarray
        SEE length [m]
    lpee : float or numpy.ndarray
        PEE length [m]
    """
    
    # Unravel parameter values
    lpee0       = muspar['lpee0']       # [m]       PEE slack length
    lsee0       = muspar['lsee0']       # [m]       SEE slack length 
    ksee        = muspar['ksee']        # [N/m^2]   shape parameter that scales the stiffness of PEE
    kpee        = muspar['kpee']        # [N/m^2]   shape parameter that scales the stiffness of SEE
    
    # SEE
    esee    = (fsee/ksee)**(1/2)                    # [m] SEE elongation
    lsee    = esee+lsee0                            # [m] SEE length
    
    # PEE
    epee    = (fsee/kpee)**(1/2)                    # [m] PEE elongation
    lpee    = epee+lpee0                            # [m] SEE length
                                  
    # Output
    return lsee,lpee

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def fce2vce(fce, q, lcerel, muspar):
    """
    Compute CE velocity from force–velocity relationship.
    
    Parameters
    ----------
    fce : float or numpy.ndarray
        Contractile element force [N]
    q : float or numpy.ndarray
        Active state [-]
    lcerel : float or numpy-array
        Relative CE length [-]
    muspar : dict
        Muscle parameters
    
    Returns
    -------
    vce : float or numpy.ndarray
        CE velocity [m/s]
    vcerel : float or numpy.ndarray
        Relative CE velocity (CE velocity / Lceopt) [1/s]
        Note: computed only if optimum CE length is known, else value is None
    """
    
    # Unravel muscle parameters
    a_c, b_c            = muspar['a'], muspar['b']
    fasymp, fmax        = muspar['fasymp'], muspar['fmax']
    slopfac, vfactmin   = muspar['slopfac'], muspar['vfactmin']
    q0                  = muspar['q0']  
    sloplin = vfactmin*b_c/(slopfac*0.005*0.0975*(fmax+a_c))
    
    # Scale arel and brel if necessary
    fisomrel = force_length(lcerel,muspar)[0]
    
    # Smooth version of KvS brel(q)
    q0_b = (np.log(1/vfactmin-1)+q0*22)/22
    b = b_c/(1+np.exp(-22*(q-q0_b)))
    
    # Scale a
    a = a_c
    a = (lcerel>1)*a*fisomrel + (lcerel<=1)*a

    # Variables for various part of vce-fce relation
    dvdf_isom_con = b/(q*(fisomrel*fmax+a)) # slope in the isometric point at wrt concentric part
    dvdf_isom_ecc = dvdf_isom_con/slopfac # slope in the isometric point at wrt eccentric part
    dFdvcon0      = 1/dvdf_isom_con
    s_as          = 1/sloplin
    p1 = -(fisomrel*q*(fasymp*fmax - fmax))/(s_as - dFdvcon0*slopfac) 
    p2 =  (fisomrel**2*q**2*(fasymp*fmax - fmax)**2)/(s_as - dFdvcon0*slopfac)
    p3 =  -fasymp*fisomrel*q*fmax;
    p4 =  -s_as

    # Compute different regions
    r_c1 = (((fce/q) <= fisomrel*fmax) * (dvdf_isom_con<=sloplin)) # Concentric, dvdf_isom_con<=sloplin (normal)    
    r_c2 = (((fce/q) <= fisomrel*fmax) * (dvdf_isom_con>sloplin)) # Concentric dvdf_isom_con>sloplin (defective case)
    r_e1 = (((fce/q) > fisomrel*fmax) * (dvdf_isom_ecc<=(sloplin/slopfac))) # Eccentric, dvdf_isom_ecc<=sloplin (normal) 
    r_e2 = (((fce/q) > fisomrel*fmax) * (dvdf_isom_ecc>(sloplin/slopfac))) # Eccentric, dvdf_isom_ecc>sloplin (defective case)
    
    #  Compute CE velocity
    vce_c1 = (b*(fce-q*fisomrel*fmax)/(fce+q*a)) # Concentric, dvdf_isom_con<=sloplin (normal)  
    vce_c2 = (sloplin*(fce-q*fisomrel*fmax)) # Concentric dvdf_isom_con>sloplin (defective case)        
    vce_e1 = ((-(fce + p3 + p1*p4 + (fce**2 - 2*fce*p1*p4 + 2*fce*p3 + p1**2*p4**2 - 2*p1*p3*p4 + p3**2 + 4*p2*p4)**(1/2))/(2*p4))) # Eccentric, dvdf_isom_ecc<=sloplin (normal)
    vce_e2 = ((sloplin/slopfac)*(fce-q*fisomrel*fmax)) # Eccentric, dvdf_isom_ecc>sloplin (defective case)
    
    # Output   
    vce = r_c1*vce_c1 + r_c2*vce_c2 + r_e1*vce_e1 + r_e2*vce_e2        
        
    if 'lce_opt' in muspar:
        vcerel = vce/muspar['lce_opt']
    else:
        vcerel = None
    
    return vce, vcerel, [vce_c1, vce_c2, vce_e1, vce_e2], [r_c1, r_c2, r_e1, r_e2]

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def vce2fce(vce,q,lcerel,muspar):
    """
    Compute CE force from force–velocity relationship.
    
    Parameters
    ----------
    vce : float or numpy.ndarray
        CE velocity [m/s]
    q : float or numpy.ndarray
        Active state [-]
    lcerel : float or numpy-array
        Relative CE length [-]
    muspar : dict
        Muscle parameters
    
    Returns
    -------
    fce : float or numpy.ndarray
        Contractile element force [N]
    fcerel : float or numpy.ndarray
        Relative CE force (CE force / Fcemax) [ ]
        Note: computed only if optimum CE length is known, else value is None
    """
            
    # Unravel parameter values
    a_c, b_c            = muspar['a'], muspar['b']
    fasymp, fmax        = muspar['fasymp'], muspar['fmax']
    slopfac, vfactmin   = muspar['slopfac'], muspar['vfactmin']
    q0                  = muspar['q0']  
    sloplin = vfactmin*b_c/(slopfac*0.005*0.0975*(fmax+a_c))
    
    # Scale arel and brel if necessary
    fisomrel = force_length(lcerel,muspar)[0]
    
    # Smooth version of KvS brel(q)
    q0_b = (np.log(1/vfactmin-1)+q0*22)/22
    b = b_c/(1+np.exp(-22*(q-q0_b)))
    
    # Scale a
    a = a_c
    a = (lcerel>1)*a*fisomrel + (lcerel<=1)*a
    
    # Variables for various part of vce-fce relation
    dvdf_isom_con = b/(q*(fisomrel*fmax+a)) # slope in the isometric point at wrt concentric part
    dvdf_isom_ecc = dvdf_isom_con/slopfac # slope in the isometric point at wrt eccentric part
    dFdvcon0      = 1/dvdf_isom_con
    s_as          = 1/sloplin
    p1 = -(fisomrel*q*(fasymp*fmax - fmax))/(s_as - dFdvcon0*slopfac) 
    p2 =  (fisomrel**2*q**2*(fasymp*fmax - fmax)**2)/(s_as - dFdvcon0*slopfac)
    p3 =  -fasymp*fisomrel*q*fmax;
    p4 =  -s_as
    
    # Compute different regions
    r_c1 = ((vce<=0) * (dvdf_isom_con<=sloplin)) # Concentric, dvdf_isom_con<=sloplin (normal)  
    r_c2 = ((vce<=0) * (dvdf_isom_con>sloplin)) # Concentric dvdf_isom_con>sloplin (defective case)
    r_e1 = ((vce>0) * (dvdf_isom_ecc<=(sloplin/slopfac))) # Eccentric, dvdf_isom_ecc<=sloplin (normal) 
    r_e2 = ((vce>0) * (dvdf_isom_ecc>(sloplin/slopfac))) # Eccentric, dvdf_isom_ecc>sloplin (defective case)
    
    # Compute CE force
    fce_c1 = (q*(b*fisomrel*fmax+a*vce)) / (b-vce) # Concentric, dvdf_isom_con<=sloplin (normal)
    fce_c2 = q*fisomrel*fmax+vce/sloplin # Concentric dvdf_isom_con>sloplin (defective case)
    fce_e1 = (p2-(p3+p4*vce)*(p1+vce))/(p1+vce) # Eccentric, dvdf_isom_con<=sloplin (normal)
    fce_e2 = q*fisomrel*fmax+((vce*slopfac)/sloplin) # Eccentric, dvdf_isom_con>sloplin (defective case)

    # Output
    fce = r_c1*fce_c1 + r_c2*fce_c2 + r_e1*fce_e1 + r_e2*fce_e2
    if 'fmax' in muspar:
        fcerel = fce/muspar['fmax']
    else:
        fcerel = None
    
    return fce, fcerel, [fce_c1, fce_c2, fce_e1, fce_e2], [r_c1, r_c2, r_e1, r_e2]

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def force_eq(lmtc,gamma,muspar):
    """
    Finds relative CE length such that SEE force equals the sum of CE and PEE
    force, for any given MTC-length and gamma (i.e., normalised concentration
    Ca2+ between the filaments)
    
    Parameters
    ----------
    lmtc : float or numpy.ndarray
        muscle-tendon-complex length [m]
    gamma : float or numpy.ndarray
        normalised concentration Ca2+ between the filaments [-]
    muspar : dict
        Muscle parameters
    
    Returns
    -------
    fsee : float or numpy.ndarray
        SEE force [N]
    lcerel : float or numpy.ndarray
        relative CE length, i.e. CE length divided by optimum CE length) [-]
    fce : float or numpy.ndarray
        CE force [N]
    fpee : float or numpy.ndarray
        PEE force [N]
    """
        
    # Unravel parameter values
    lce_opt = muspar['lce_opt']     # [m] CE optimum length (i.e., CE length at which isometric CE force is maximal)
    lsee0   = muspar['lsee0']       # [m] SEE slack length
    fmax    = muspar['fmax']        # [N] maximal CE force
    C       = -1/muspar['w']**2     # [ ] parameter based on width of 2nd order polynomial of isometric CE force-length relation
    ksee    = muspar['ksee']        # [N/m^2] shape parameter that scales the stiffness of SEE
    
    # Get initial guess of relative CE length   
    q = (gamma>muspar['q0'])*act_state(gamma,1,muspar)[0] + (gamma<=muspar['q0'])*muspar['q0'] # [ ]
    # lcerel =  -(ksee*lce_opt*lsee0 - ksee*lce_opt*lmtc + fmax**(1/2)*q**(1/2)*(ksee*C*lce_opt**2 + ksee*C*lmtc**2 + ksee*C*lsee0**2 + ksee*lce_opt**2 - C*fmax*q - 2*ksee*C*lce_opt*lmtc + 2*ksee*C*lce_opt*lsee0 - 2*ksee*C*lmtc*lsee0)**(1/2) + C*fmax*q)/(ksee*lce_opt**2 - C*fmax*q)  # [ ]
    a = ksee*lce_opt*lsee0 - ksee*lce_opt*lmtc + C*fmax*q
    b1 = fmax**(1/2)*q**(1/2)
    b2 = ksee*C*lce_opt**2 + ksee*C*lmtc**2 + ksee*C*lsee0**2 + ksee*lce_opt**2 - C*fmax*q - 2*ksee*C*lce_opt*lmtc + 2*ksee*C*lce_opt*lsee0 - 2*ksee*C*lmtc*lsee0
    b2 = (b2>=0)*b2 + (b2<0)*0   
    b = b1*b2**(1/2) 
    c = ksee*lce_opt**2 - C*fmax*q
    lcerel = -(a+b)/c
    
    # Set values for first round
    fce, fpee, fsee, dlcerel = 1e6, 0, 0, 0
    # Set tolerance
    tolF = 1e-5*muspar['fmax'] # until 0.01% of Fmax
    
    # Newton root finding to find relative CE length
    while (np.max(np.abs(fce+fpee-fsee))>tolF):
        lcerel = lcerel+dlcerel  # [ ]
        fisomrel,Kcerel = force_length(lcerel,muspar)  # [ , ]
        q,Kq = act_state(gamma,lcerel,muspar)[0:2] # [ , ]
        fce = q*fisomrel*muspar['fmax'] # [N]
        Kce = fisomrel*Kq*muspar['fmax']+q*Kcerel*muspar['fmax'] # [N]
        lce = lcerel*muspar['lce_opt']
        lsee = lmtc-lce
        lpee = lce
        fsee,fpee,_,_,Ksee,Kpee = lee2force(lsee,lpee,muspar)[0:6] # [N, N, N/m, N/m]
        Ksee = -Ksee*muspar['lce_opt'] # [N]
        Kpee = Kpee*muspar['lce_opt'] # [N]
        dlcerel = (fce+fpee-fsee)/(Ksee-Kce-Kpee) # [ ]
        pass
    
    # Output
    return fsee,lcerel,fce,fpee

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def anly_gamma(time, cf, t_stim_on, t_stim_off, stim, muspar):
    """
    Calculate gamma over time for constant muscle stimulation with periodic behavior.

    Parameters
    ----------
    time : np.ndarray
        Time axis.
    cf : float
        Cycle frequency in Hz.
    t_stim_on : float
        Time at which muscle stimulation switches to the constant value of `stim`.
    t_stim_off : float
        Time at which muscle stimulation switches to 0.
    stim : float
        Value of constant muscle stimulation.
    muspar : dict
        Dictionary containing:
        - 'gamma_0': Minimum value of gamma.
        - 'tact': Time constant for activation.
        - 'tdeact': Time constant for deactivation.

    Returns
    -------
    gamma : np.ndarray
        Gamma as a function of time.
    stim_t : np.ndarray
        Stimulation as a function of time.
    """
    gamma_0 = muspar['gamma_0']
    tact = muspar['tact']
    tdeact = muspar['tdeact']

    # Cycle duration
    Tcycle = 1 / cf

    # Shift time with t_stim_on and take modulus
    time_mod = np.mod(time - t_stim_on, Tcycle)
    t_stim_off_mod = np.mod(t_stim_off - t_stim_on, Tcycle)
    t_stim_on_mod = 0

    # Calculate gamma at t=0 for periodic behavior
    gamma0_1 = (
        gamma_0 +
        np.exp(-Tcycle / tdeact) *
        np.exp(t_stim_off_mod / tdeact) *
        (gamma_0 - 1) *
        (stim * np.exp(-t_stim_off_mod / tact) - stim +
         (gamma_0 * stim * np.exp(-t_stim_off_mod / tact)) / (stim - gamma_0 * stim))
    ) / (
        (stim * np.exp(-Tcycle / tdeact) *
         np.exp(-t_stim_off_mod / tact) *
         np.exp(t_stim_off_mod / tdeact) *
         (gamma_0 - 1)) / (stim - gamma_0 * stim) + 1
    )

    # Calculate gamma at the end of stimulation phase
    gamma0_2 = gamma_rise(t_stim_off_mod, gamma0_1, stim, muspar)

    # Compute gamma(t)
    gamma = np.where(
        (time_mod >= 0) & (time_mod < t_stim_off_mod),
        gamma_rise(time_mod, gamma0_1, stim, muspar),
        gamma_relax(time_mod - t_stim_off_mod, gamma0_2, muspar)
    )

    # Compute stimulation over time
    stim_t = np.where((time_mod >= 0) & (time_mod < t_stim_off_mod), stim, 0)

    return gamma, stim_t


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def gamma_relax(time, gamma0, muspar):
    """
    Calculate gamma over time during relaxation (stim = 0).

    Parameters
    ----------
    time : np.ndarray
        Time in seconds.
    gamma0 : float
        Initial gamma at t=0.
    muspar : dict
        Dictionary containing:
        - 'gamma_0': Minimum value of gamma.
        - 'tdeact': Time constant for deactivation.

    Returns
    -------
    gamma : np.ndarray
        Relative amount of Ca2+ between filaments during relaxation.
    """
    gamma_0 = muspar['gamma_0']
    tdeact = muspar['tdeact']

    t_shift = -np.log((gamma_0 - gamma0) / (gamma_0 - 1)) * tdeact
    gamma = (1 - gamma_0) * np.exp(-(time + t_shift) / tdeact) + gamma_0

    return gamma


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def gamma_rise(time, gamma0, stim, muspar):
    """
    Calculate gamma over time during activation (stim > 0).

    Parameters
    ----------
    time : np.ndarray
        Time in seconds.
    gamma0 : float
        Initial gamma at t=0.
    stim : float
        Normalized muscle stimulation.
    muspar : dict
        Dictionary containing:
        - 'gamma_0': Minimum value of gamma.
        - 'tact': Time constant for activation.

    Returns
    -------
    gamma : np.ndarray
        Relative amount of Ca2+ between filaments during activation.
    """
    gamma_0 = muspar['gamma_0']
    tact = muspar['tact']

    t_shift = -np.log((gamma_0 - gamma0) / (stim - gamma_0 * stim) + 1) * tact
    gamma = stim * (1 - gamma_0) * (1 - np.exp(-(time + t_shift) / tact)) + gamma_0

    return gamma

interpolation.py

"""
This module provides interpolation functions for structured data grids.

The functions interpolate 1D curves, 2D grids and 3D grids onto finer grids.
Missing values (NaNs) are removed before interpolation.

Functions
---------
do_2d(x, y, **kwargs)
    Interpolate 1D data (x, y) onto a finer grid.

do_3d(data, grid, **kwargs)
    Interpolate 2D grid data onto a finer 2D grid.

do_4d(data, grid, **kwargs)
    Interpolate 3D grid data onto a finer 3D grid.
"""

import numpy as np
from scipy.interpolate import interp1d, griddata

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def do_2d(x, y, **kwargs):
    """
    Interpolate 1D data onto a finer grid.

    Parameters
    ----------
    x : array-like
        1D array of x-coordinates.
    y : array-like
        1D array of y-values corresponding to `x`.
    **kwargs : dict, optional
        Additional keyword arguments:
        
        N : int, optional
            Number of interpolation points (default is 100).
        method : str, optional
            Interpolation method passed to `scipy.interpolate.interp1d`
            (default is 'cubic').

    Returns
    -------
    x_fine : ndarray
        Interpolated x-coordinates.
    y_fine : ndarray
        Interpolated y-values.

    Notes
    -----
    NaN values in `y` are removed prior to interpolation.
    """
    num_points = kwargs.get("N", 100)
    method = kwargs.get("method", "cubic")

    x = np.asarray(x)
    y = np.asarray(y)

    # Remove NaNs
    valid_mask = ~np.isnan(y)
    x_valid = x[valid_mask]
    y_valid = y[valid_mask]

    # Fine grid
    x_fine = np.linspace(x[0], x[-1], num_points)

    interpolator = interp1d(
        x_valid,
        y_valid,
        kind=method,
        bounds_error=False,
    )
    y_fine = interpolator(x_fine)

    return x_fine, y_fine

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def do_3d(data, grid, **kwargs):
    """
    Interpolate 2D grid data onto a finer grid.

    Parameters
    ----------
    data : ndarray
        2D array of data values.
    grid : tuple of array-like
        Tuple (x, y) defining the grid coordinates.
    **kwargs : dict, optional
        Additional keyword arguments:
        
        N : int, optional
            Number of interpolation points along x (default is 100).
        method : str, optional
            Interpolation method for `scipy.interpolate.griddata`
            (default is 'linear').

    Returns
    -------
    data_fine : ndarray
        Interpolated data on the finer grid.
    grid_fine : tuple of ndarray
        Tuple (x_fine, y_fine) defining the new grid.

    Notes
    -----
    NaN values in `data` are removed prior to interpolation.
    """
    num_points = kwargs.get("N", 100)
    method = kwargs.get("method", "linear")

    x, y = grid

    x = np.asarray(x)
    y = np.asarray(y)
    data = np.asarray(data)

    # Meshgrid
    x_mesh, y_mesh = np.meshgrid(x, y, indexing="ij")

    # Flatten
    x_flat = x_mesh.ravel()
    y_flat = y_mesh.ravel()
    data_flat = data.ravel()

    # Remove NaNs
    valid_mask = ~np.isnan(data_flat)
    x_valid = x_flat[valid_mask]
    y_valid = y_flat[valid_mask]
    data_valid = data_flat[valid_mask]

    # Fine grid
    x_fine = np.linspace(x[0], x[-1], num_points)
    y_fine = np.linspace(y[0], y[-1], num_points + 1)
    x_fine_mesh, y_fine_mesh = np.meshgrid(x_fine, y_fine, indexing="ij")

    interpolated_values = griddata(
        points=np.vstack([x_valid, y_valid]).T,
        values=data_valid,
        xi=np.vstack([x_fine_mesh.ravel(), y_fine_mesh.ravel()]).T,
        method=method,
    )

    data_fine = interpolated_values.reshape((num_points, num_points + 1))
    grid_fine = (x_fine, y_fine)

    return data_fine.T, grid_fine

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def do_4d(data, grid, **kwargs):
    """
    Interpolate 3D grid data onto a finer grid.

    Parameters
    ----------
    data : ndarray
        3D array of data values.
    grid : tuple of array-like
        Tuple (x, y, z) defining the grid coordinates.
    **kwargs : dict, optional
        Additional keyword arguments:
        
        N : int, optional
            Number of interpolation points along x (default is 100).
        method : str, optional
            Interpolation method for `scipy.interpolate.griddata`
            (default is 'linear').

    Returns
    -------
    data_fine : ndarray
        Interpolated 3D data.
    grid_fine : tuple of ndarray
        Tuple (x_fine, y_fine, z_fine) defining the new grid.

    Notes
    -----
    NaN values in `data` are removed prior to interpolation.
    """
    num_points = kwargs.get("N", 100)
    method = kwargs.get("method", "linear")

    x, y, z = grid

    x = np.asarray(x)
    y = np.asarray(y)
    z = np.asarray(z)
    data = np.asarray(data)

    # Meshgrid
    x_mesh, y_mesh, z_mesh = np.meshgrid(x, y, z, indexing="ij")

    # Flatten
    x_flat = x_mesh.ravel()
    y_flat = y_mesh.ravel()
    z_flat = z_mesh.ravel()
    data_flat = data.ravel()

    # Remove NaNs
    valid_mask = ~np.isnan(data_flat)
    x_valid = x_flat[valid_mask]
    y_valid = y_flat[valid_mask]
    z_valid = z_flat[valid_mask]
    data_valid = data_flat[valid_mask]

    # Fine grid
    x_fine = np.linspace(x[0], x[-1], num_points)
    y_fine = np.linspace(y[0], y[-1], num_points + 1)
    z_fine = np.linspace(z[0], z[-1], num_points + 2)

    x_fine_mesh, y_fine_mesh, z_fine_mesh = np.meshgrid(
        x_fine, y_fine, z_fine, indexing="ij"
    )

    interpolated_values = griddata(
        points=np.vstack([x_valid, y_valid, z_valid]).T,
        values=data_valid,
        xi=np.vstack(
            [
                x_fine_mesh.ravel(),
                y_fine_mesh.ravel(),
                z_fine_mesh.ravel(),
            ]
        ).T,
        method=method,
    )

    data_fine = interpolated_values.reshape(
        (num_points, num_points + 1, num_points + 2)
    )
    grid_fine = (x_fine, y_fine, z_fine)

    return data_fine, grid_fine

stats.py

# -*- coding: utf-8 -*-
"""
This module provides small statistical and formatting helper functions.

Functions
---------
ceil(a, precision=0)
    Round values up to a given decimal precision.

floor(a, precision=0)
    Round values down to a given decimal precision.

rmse(x, y)
    Compute the root mean square error between two arrays.

pdiff(x, y)
    Compute the percentage difference of `x` relative to `y`.

zscore(x, axis=-1)
    Compute z-score normalisation along a given axis.

str_round(value, n)
    Round a value to a fixed number of significant digits and return a string.

find_max(data, grid)
    Find the maximum value in an array and return its grid coordinates.

analyse_3similar(lst, tolerance=1.0)
    Compare three values and identify whether two or more are similar.
"""

import numpy as np

def ceil(a, precision=0):
    """
    Ceil a number or array to the given decimal precision.

    Parameters
    ----------
    a : float or array-like
        Input value(s).
    precision : int, optional
        Number of decimal places. Default is 0.

    Returns
    -------
    float or numpy.ndarray
        Value(s) with ceiling applied at the requested precision.
    """
    return np.true_divide(np.ceil(a * 10**precision), 10**precision)

def floor(a, precision=0):
    """
    Floor a number or array to the given decimal precision.

    Parameters
    ----------
    a : float or array-like
        Input value(s).
    precision : int, optional
        Number of decimal places. Default is 0.

    Returns
    -------
    float or numpy.ndarray
        Value(s) with floor applied at the requested precision.
    """
    return np.true_divide(np.floor(a * 10**precision), 10**precision)

def rmse(x, y):
    """
    Compute the root mean square error between two arrays.

    Parameters
    ----------
    x, y : array-like
        Input arrays with compatible shapes.

    Returns
    -------
    float
        Root mean square error.
    """
    mse = np.mean((x - y) ** 2)
    return float(np.sqrt(mse))

def pdiff(x,y):
    """
    Compute the percentage difference of `x` relative to `y`.

    Parameters
    ----------
    x, y : array-like or scalar
        Values with compatible shapes.

    Returns
    -------
    float or numpy.ndarray
        Percentage difference, computed as `(x - y) / y * 100`.
    """
    return ((x-y)/y)*100

def zscore(x, axis=-1):
    """
    Compute the z-score normalization along a given axis.

    Parameters
    ----------
    x : np.ndarray
        Input array.
    axis : int, optional
        Axis along which to compute the mean and std. Default is -1.

    Returns
    -------
    np.ndarray
        Z-score normalized array.
    """
    m = np.mean(x, axis=axis, keepdims=True)
    s = np.std(x, axis=axis, keepdims=True)
    return (x - m) / s

def str_round(value: float, n: int) -> str:
    """
    Round a number to n significant digits and return a string representation.

    Parameters
    ----------
    value : float
        Value to round.
    n : int
        Number of significant digits.

    Returns
    -------
    str
        Rounded value formatted as a string. Returns '-' for NaN.
    """
    
    if np.isnan(value):
        return '-'
    elif value == 0:
        return '0'
    
    # Determine the number of digits before the decimal point
    num_digits = int(np.floor(np.log10(abs(value)))) + 1
    
    # Calculate the decimal places to round
    decimal_places = max(n - num_digits, 0)
    
    rounded_value = round(value, decimal_places)
    
    # Format as string with fixed decimal places if needed
    if decimal_places > 0:
        return f"{rounded_value:.{decimal_places}f}"
    else:
        return str(int(rounded_value))
    
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def find_max(data, grid):
    """
    Find the maximum value in an N-dimensional array and its corresponding values 
    from the provided grid arrays.

    Parameters
    ----------
    data : ndarray
        N-dimensional array containing the data.
    grid : tuple of ndarray
        Tuple of 1D arrays representing the coordinate axes for each dimension of `data`.
        Length must match the number of dimensions of `data`.

    Returns
    -------
    data_max : float
        Maximum value found in `data`.
    grid_max : tuple
        Coordinates corresponding to the maximum value, taken from `grid`.
        The order matches the order of `grid`.
    
    Notes
    -----
    Uses `np.nanargmax` to ignore NaN values in `data`.
    """
    max_idx = np.nanargmax(data)
    
    if len(grid) == 2:
        x, y = grid
        row_idx, col_idx = np.unravel_index(max_idx, data.shape)
        grid_max = (x[col_idx], y[row_idx])
        data_max = data[row_idx, col_idx]
        
    elif len(grid) == 3:
        x, y, z = grid
        row_idx, col_idx, dep_idx = np.unravel_index(max_idx, data.shape)
        grid_max = (x[row_idx], y[col_idx], z[dep_idx])
        data_max = data[row_idx, col_idx, dep_idx]
        
    else:
        raise ValueError("Grid must have length 2 or 3 corresponding to data dimensions.")
    
    return data_max, grid_max
    
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def analyse_3similar(lst, tolerance=1.0):
    """
    Analyze a list of three values to determine similarity within a tolerance.

    The function compares up to three values and identifies whether any values
    are approximately equal within a specified tolerance. It handles missing
    values (`None` or `np.nan`) and returns the most representative value,
    a differing value (if any), and an indicator of similarity.

    Parameters
    ----------
    lst : list of float or None
        A list containing exactly three elements. Elements may be numeric,
        `None`, or `np.nan`.
    tolerance : float, optional
        The maximum absolute difference for two values to be considered
        similar. Default is 1.0.

    Returns
    -------
    most_frequent : float or None
        The value considered most representative among the inputs. If no
        conclusion can be drawn, returns None.
    different_value : float or None
        A value that differs from the most representative value. Returns
        None if all valid values are similar or insufficient data exists.
    flag : bool or int
        - If bool:
            * True  -> a consistent or usable result was found
            * False -> no meaningful similarity detected
        - If int:
            Index (0, 1, or 2) of the value considered different when
            exactly one value deviates.

    Raises
    ------
    ValueError
        If `lst` does not contain exactly three elements.

    Notes
    -----
    - Missing values (`None` or `np.nan`) are ignored in comparisons.
    - If only one valid value is present, it is returned as the most
      representative value.
    - If two valid values are present, they are compared directly.
    - If all three values are valid, pairwise comparisons determine whether
      a majority agreement exists.

    Examples
    --------
    >>> analyse_3similar([1.0, 1.1, 0.9], tolerance=0.2)
    (1.0, None, True)

    >>> analyse_3similar([1.0, 5.0, 1.1], tolerance=0.2)
    (1.0, 5.0, 1)

    >>> analyse_3similar([None, 2.0, np.nan])
    (2.0, None, True)

    >>> analyse_3similar([1.0, 2.0, 3.0])
    (None, None, False)
    """
    
    if len(lst) != 3:
        raise ValueError("List must have exactly 3 items.")
    
    # Store values with their indices
    valid = [(i, v) for i, v in enumerate(lst) if not (v is None or isinstance(v, float) and np.isnan(v))]
    
    if len(valid) == 0:
        return None, None, False  # No valid data
    if len(valid) == 1:
        return valid[0][1], None, True  # Only one valid value, treat as "most frequent"

    def close(x, y):
        return abs(x - y) <= tolerance

    # Unpack valid items
    (i1, v1), (i2, v2) = valid[0], valid[1]
    
    # If only two valid values
    if len(valid) == 2:
        if close(v1, v2):
            return v1, None, True
        else:
            return v1, v2, i2  # Arbitrarily say v1 is "most frequent"

    # All three are valid
    a, b, c = lst
    if close(a, b) and close(a, c) and close(b, c):
        return a, None, True
    if close(a, b):
        return a, c, 2
    elif close(a, c):
        return a, b, 1
    elif close(b, c):
        return b, a, 0

    return None, None, False

stimulation.py

"""
This module provides functions for analysing stimulation over time.

The functions detect stimulation onset and offset times from time-series
signals and compute stimulation durations from simulation or experimental CSV
files.

Functions
---------
get_stim_timing(time, signal)
    Detect stimulation pulse trains in a signal and return their onset and
    offset time.

get_stim_dur(filepaths)
    Compute mean stimulation duration(s) from one or multiple CSV files
    containing time-series stimulation data.
"""

import numpy as np
import pandas as pd

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def get_stim_timing(time, signal):
    """
    Detect stimulation onset and offset time in a signal.

    This function identifies the onset and offset times of pulse trains in a
    time-series signal by thresholding and grouping consecutive pulses into
    trains based on temporal gaps.

    Parameters
    ----------
    time : array-like
        1D array of time values corresponding to `signal`.
    signal : array-like
        1D array of signal values representing stimulation over time.

    Returns
    -------
    t_stim_on : numpy.ndarray
        Array of onset times for each detected pulse train.
    t_stim_off : numpy.ndarray
        Array of offset times for each detected pulse train.

    Notes
    -----
    - A threshold is computed as the midpoint between the minimum and maximum
      of the signal to binarize it.
    - Pulse starts and stops are detected using differences in the binarized
      signal.
    - Pulse trains are defined by grouping pulses separated by gaps smaller
      than a characteristic pulse width.
    - If unusually large pulse widths are detected (e.g., >100 samples),
      the grouping threshold is reset.
    - The helper function `createBlockSignal` is used to construct the output
      block signal.

    Examples
    --------
    >>> t_on, t_off = get_stim_timing(time, signal)

    >>> t_on
    array([0.5, 2.0])

    >>> t_off
    array([1.0, 2.5])
    """
    # Determine threshold as midpoint between min and max
    threshold = (np.max(signal) + np.min(signal)) / 2

    # Binarize signal
    binary_signal = np.where(signal > threshold, 1, 0)

    # Find pulse start and stop indices
    i_start_pulse = np.where(np.diff(binary_signal, prepend=0) == 1)[0]
    i_stop_pulse = np.where(np.diff(binary_signal, prepend=0) == -1)[0]

    if len(i_start_pulse) == 0 or len(i_stop_pulse) == 0:
        return [], [], []

    # Estimate characteristic pulse width
    i_start_pulse_same = i_start_pulse
    if i_start_pulse_same[0] == 0:
        i_start_pulse_same = i_start_pulse_same[1:]

    pulse_widths = i_stop_pulse - i_start_pulse_same
    pulse_width_median = np.median(pulse_widths)

    # Handle unusually large pulse widths
    if any(abs(pulse_widths) > 100):
        pulse_width_median = 0

    # Group pulses into trains
    i_start_train = [int(i_start_pulse[0])]
    i_stop_train = []

    for i in range(1, len(i_start_pulse)):
        if i_start_pulse[i] - i_stop_pulse[i - 1] > pulse_width_median:
            i_stop_train.append(int(i_stop_pulse[i - 1]))
            i_start_train.append(int(i_start_pulse[i]))

    # Add final stop index
    i_stop_train.append(int(i_stop_pulse[-1]))

    # Output
    t_stim_on = time[i_start_train]
    t_stim_off = time[i_stop_train]

    return t_stim_on, t_stim_off

def get_stim_dur(filepaths):
    """
    Compute mean stimulus duration(s) from one or more CSV files.

    Each file is expected to contain time-series data where stimulation onset
    and offset can be detected using `get_stim_timing`. The function extracts
    stimulus durations and returns their mean for each file. If a file cannot
    be processed, NaN is returned for that entry.

    Parameters
    ----------
    filepaths : str or list of str
        Path or list of paths to CSV file(s) containing the data.

    Returns
    -------
    stim_durations : numpy.ndarray
        Array of mean stimulus durations for each input file. If a file
        cannot be processed, the corresponding value is `np.nan`.

    Notes
    -----
    - The CSV file is read using `pandas.read_csv`, transposed, and converted
      to a NumPy array.
    - The function assumes the data contains at least three rows corresponding
      to time and stimulus signal.
    - The helper function `get_stim_timing` must return stimulation onset and
      offset times.

    Examples
    --------
    >>> get_stim_dur("data.csv")
    array([0.52])

    >>> get_stim_dur(["file1.csv", "file2.csv"])
    array([0.52, nan])
    """
    
    if isinstance(filepaths, str):
        filepaths = [filepaths]

    stim_durations = []

    for filepath in filepaths:
        try:
            data = pd.read_csv(filepath).T.to_numpy()
            time, _, stim, *_ = data

            t_stim_on, t_stim_off = get_stim_timing(
                time[0:-1], stim[0:-1]
            )

            stim_dur = t_stim_off - t_stim_on
            stim_durations.append(np.mean(stim_dur))

        except Exception:
            stim_durations.append(np.nan)

    return np.array(stim_durations)

trajectories.py

"""
This module provides functions to generate prescribed MTC length trajectories.

Functions
---------
cv(time, cf, fts, mle, lmtc_avg)
    Generate a cyclic trajectory with constant shortening and lengthening
    velocities.

scv(time, cf, fts, amp, lmtc_avg, acc)
    Generate a cyclic trajectory with smoothed acceleration around the turning
    points and constant velocity between them.
"""

import numpy as np

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def cv(time, cf, fts, mle, lmtc_avg):
    """
    Generate muscle-tendon complex length and velocity over time.

    This function computes the length (`lmtc`) and velocity (`vmtc`) of a
    muscle-tendon complex assuming a cyclic motion with constant shortening
    and lengthening velocities within each cycle.

    Parameters
    ----------
    time : array_like
        Time vector [s].
    cf : float
        Cycle frequency [Hz].
    fts : float
        Fraction of the cycle spent in shortening phase [-].
    mle : float
        Total muscle-tendon length excursion [m].
    lmtc_avg : float
        Average muscle-tendon complex length [m].

    Returns
    -------
    lmtc : ndarray
        Muscle-tendon complex length over time [m].
    vmtc : ndarray
        Muscle-tendon complex velocity over time [m/s].

    Raises
    ------
    ZeroDivisionError
        If `cf` is zero.
    ValueError
        If `fts` is not between 0 and 1.

    Examples
    --------
    >>> import numpy as np
    >>> t = np.linspace(0, 1, 100)
    >>> lmtc, vmtc = cv(t, 1.0, 0.4, 0.1, 1.0)
    """
    if not 0 <= fts <= 1:
        raise ValueError("fraction_time_shortening must be between 0 and 1.")

    t_mod = np.mod(time, 1 / cf)

    # Calculate shortening and lengthening times
    t_short = fts / cf  # [s]
    t_length = (1 - fts) / cf  # [s]

    # Calculate constant shortening and lengthening velocities
    v_short = -mle / t_short
    v_length = mle / t_length

    # Compute length
    lmtc = (
        (t_mod < t_short) *
        (lmtc_avg + mle / 2 + v_short * t_mod)
        + (t_mod >= t_short) *
        (lmtc_avg - mle / 2 + v_length * (t_mod - t_short))
    )

    # Compute velocity
    vmtc = (
        (t_mod < t_short) * v_short
        + (t_mod >= t_short) * v_length
    )

    return lmtc, vmtc

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def scv(time, cf, fts, amp, lmtc_avg, acc):
    """
    Compute muscle–tendon complex (MTC) kinematics for a prescribed trajectory.

    The trajectory consists of shortening and lengthening phases with
    sinusoidal acceleration profiles near the turning points and constant
    velocity in between.

    Parameters
    ----------
    time : float or ndarray
        Time point(s) at which to evaluate the trajectory [s].
    cf : float
        Cycle frequency [Hz].
    fts : float
        Fraction of the cycle spent shortening [-].
    amp : float
        MTC amplitude [m].
    lmtc_avg : float
        Mean MTC length around which the motion is centered [m].
    acc : float
        Peak acceleration magnitude near turning points [m/s^2].

    Returns
    -------
    lmtc : float or ndarray
        MTC length at `time` [m].
    vmtc : float or ndarray
        MTC velocity at `time` [m/s].
    amtc : float or ndarray
        MTC acceleration at `time` [m/s^2].

    Notes
    -----
    The motion is periodic with period 1 / cf and consists of six phases:
    acceleration–constant–deceleration for both shortening and lengthening.
    """

    # Calculate shortening and lengthening times
    t_short = fts / cf  # [s] shortening time
    t_length = (1 - fts) / cf  # [s] lengthening time

    # Check if inputs are feasible
    # if acc * (acc * np.pi**2 * tLeng**2 + 64 * amp - 16 * amp * np.pi**2) < 0 or acc * np.pi**2 * tShort**2 + 64 * amp - 16 * amp * np.pi**2 < 0:
    #     acc = amp*(16 * np.pi**2 - 64)/(np.pi**2 * tLeng**2)
    #     return None, None, None  # Exit early if the condition is not feasible
        
    # Calculate constant shortening and lengthening velocity
    v_short = (np.pi * np.sqrt(acc * (acc * np.pi**2 * t_short**2 + 64 * amp - 16 * amp * np.pi**2)) - acc * t_short * np.pi**2) / (4 * (np.pi**2 - 4))
    v_length = -(np.pi * np.sqrt(acc * (acc * np.pi**2 * t_length**2 + 64 * amp - 16 * amp * np.pi**2)) - acc * t_length * np.pi**2) / (4 * (np.pi**2 - 4))

    # Calculate acceleration times
    t_acc1 = -v_short / (acc / 2)
    t_acc2 = v_length / (acc / 2)

    # Calculate constant velocity times
    t_con1 = t_short - 2 * t_acc1
    t_con2 = t_length - 2 * t_acc2
        
    # Check if inputs are feasible
    # if tCon1 < 0 or tCon2 < 0:
    #     return None, None, None  # Exit early if the condition is not feasible
    #     #  tCon2 <0: acc = (np.pi**2*(64*amp-16*amp*np.pi**2))/((16*tLeng**2-np.pi**4*tLeng**2))
    #     # maybe give it 0.1% more because of numerical precision
    
    # Calculate points in time
    t_points = np.cumsum([t_acc1, t_con1, t_acc1, t_acc2, t_con2, t_acc2])
    t1, t2, t3, t4, t5, t6 = t_points

    # Current time cycle position
    tc = np.mod(time, t6)

    # Calculate acceleration
    amtc = np.where(tc <= t1, -acc / 2 * (np.sin(np.pi / 2 + (np.pi * tc) / t1) + 1),
                    np.where(tc <= t2, 0,
                             np.where(tc <= t3, (np.sin((tc - t2) * np.pi / (t3 - t2) - 0.5 * np.pi) + 1) * acc / 2,
                                      np.where(tc <= t4, (np.sin((tc - t3) * np.pi / (t4 - t3) + 0.5 * np.pi) + 1) * acc / 2,
                                               np.where(tc <= t5, 0, 
                                                        (np.sin((tc - t5) * np.pi / (t6 - t5) - 0.5 * np.pi) + 1) * -acc / 2)))))

    # Calculate velocity
    vmtc = np.where(tc <= t1, -acc / 2 * tc + (-acc / 2 * t1 * np.sin(np.pi * tc / t1)) / np.pi,
                    np.where(tc <= t2, v_short,
                             np.where(tc <= t3, acc / 2 * (tc - t3) - (acc / 2 * np.sin(np.pi * (t2 - tc) / (t2 - t3) - np.pi) * (t2 - t3)) / np.pi,
                                      np.where(tc <= t4, acc / 2 * (tc - t3) + (acc / 2 * np.sin(np.pi * (t3 - tc) / (t3 - t4) - np.pi) * (t3 - t4)) / np.pi,
                                               np.where(tc <= t5, v_length, 
                                                        -acc / 2 * (tc - t6) + (-acc / 2 * np.sin(np.pi * (t5 - tc) / (t5 - t6)) * (t5 - t6)) / np.pi)))))

    # Initial position conditions
    lmtcP0 = lmtc_avg + amp - ((-acc / 2 * 0**2) / 2 - (-acc / 2 * t1**2 * np.cos((np.pi * 0) / t1)) / np.pi**2)
    lmtcP1 = lmtcP0 + (-acc / 2 * t1**2) / 2 - (-acc / 2 * t1**2 * np.cos((np.pi * t1) / t1)) / np.pi**2
    lmtcP2 = lmtcP1 + (t2 - t1) * v_short - ((acc / 2 * t2**2) / 2 - acc / 2 * t3 * t2 + (acc / 2 * np.cos((np.pi * (t2 - t2)) / (t2 - t3)) * (t2 - t3)**2) / np.pi**2)
    lmtcP3 = lmtcP2 + ((acc / 2 * t3**2) / 2 - acc / 2 * t3 * t3 + (acc / 2 * np.cos((np.pi * (t2 - t3)) / (t2 - t3)) * (t2 - t3)**2) / np.pi**2) - ((acc / 2 * (t3 - t3)**2) / 2 - (acc / 2 * np.cos((np.pi * (t3 - t3)) / (t3 - t4)) * (t3 - t4)**2) / np.pi**2)
    lmtcP4 = lmtcP3 + ((acc / 2 * (t3 - t4)**2) / 2 - (acc / 2 * np.cos((np.pi * (t3 - t4)) / (t3 - t4)) * (t3 - t4)**2) / np.pi**2)
    lmtcP5 = lmtcP4 + (t5 - t4) * v_length - (acc / 2 * t6 * t5 - (acc / 2 * t5**2) / 2 - (acc / 2 * np.cos((np.pi * (t5 - t5)) / (t5 - t6)) * (t5 - t6)**2) / np.pi**2)

    # Calculate position
    lmtc = np.where(tc <= t1, lmtcP0 + (-acc / 2 * tc**2) / 2 - (-acc / 2 * t1**2 * np.cos((np.pi * tc) / t1)) / np.pi**2,
                    np.where(tc <= t2, lmtcP1 + (tc - t1) * v_short,
                             np.where(tc <= t3, lmtcP2 + ((acc / 2 * tc**2) / 2 - acc / 2 * t3 * tc + (acc / 2 * np.cos((np.pi * (t2 - tc)) / (t2 - t3)) * (t2 - t3)**2) / np.pi**2),
                                      np.where(tc <= t4, lmtcP3 + ((acc / 2 * (t3 - tc)**2) / 2 - (acc / 2 * np.cos((np.pi * (t3 - tc)) / (t3 - t4)) * (t3 - t4)**2) / np.pi**2),
                                               np.where(tc <= t5, lmtcP4 + (tc - t4) * v_length,
                                                        lmtcP5 + (acc / 2 * t6 * tc - (acc / 2 * tc**2) / 2 - (acc / 2 * np.cos((np.pi * (t5 - tc)) / (t5 - t6)) * (t5 - t6)**2) / np.pi**2))))))


    return lmtc, vmtc, amtc