calc_kelvin_voigt_force｜弾性力＋粘性抵抗力を計算（多次元対応フォークトモデル、Python関数）

# NumPy配列（ndarray）を入力する場合
import numpy as np                                # 配列計算に使用するライブラリ

# 警告を出力する場合
import warnings                                   # 警告の表示に使用する標準ライブラリ

# 別のpyファイルからインポートする関数
from stiffness_force import calc_stiffness_force  # ばね成分の計算に使用する関数
from viscous_drag import calc_viscous_drag        # 粘性成分の計算に使用する関数

def calc_kelvin_voigt_force(k, c, pos_a, v_a, pos_b=None, v_b=None,
                            pos_a_init=None, pos_b_init=None, f_init=None):

    """
    フォークトモデル（ばねと粘性の並列接続）に基づく粘弾性力の計算（多次元対応）
    ばね成分は calc_stiffness_force、粘性成分は calc_viscous_drag を用いて計算する
    Args:
        k          : float | tuple[float] | list[float] | ndarray         ばね定数 [N/m] (>= 0, 有限値)
        c          : float | tuple[float] | list[float] | ndarray         粘性係数 [N・s/m] (>= 0, 有限値)
        pos_a      : float | tuple[float] | list[float] | ndarray         質点Aの現在位置 [m]（有限値）
        v_a        : float | tuple[float] | list[float] | ndarray         質点Aの現在速度 [m/s]（有限値）
        pos_b      : float | tuple[float] | list[float] | ndarray | None  質点Bの現在位置 [m]（有限値、省略時は固定壁として0。pos_b_initのみ指定した場合はpos_b_initの位置に固定）
        v_b        : float | tuple[float] | list[float] | ndarray | None  質点Bの現在速度 [m/s]（有限値、省略時は固定壁として0）
        pos_a_init : float | tuple[float] | list[float] | ndarray | None  質点Aの初期位置 [m]（有限値、省略時は0）
        pos_b_init : float | tuple[float] | list[float] | ndarray | None  質点Bの初期位置 [m]（有限値、省略時は0）
        f_init     : float | tuple[float] | list[float] | ndarray | None  質点Aがばねから受ける初期荷重 [N]（有限値、省略時は0）
        ※ ndarray は1次元配列（例：np.array([1.0, 2.0])）で指定してください。
        ※ 各要素が各成分（x, y, z）に対応します。
        ※ 0次元配列（例：np.array(1.0)）はスカラーとして扱います。
        ※ 本関数のみを単体で使用する場合は、float または ndarray のみ対応します。
        ※ list・tuple を使う場合やエラーチェックを行う場合は、次の内部関数も必要です。
           ・_to_array_1dsim_com      : 入力値（数値・list・tuple・ndarray）を1次元ndarrayに変換
           ・_format_output_1dsim_com : 計算結果の桁あふれを確認し、基準入力と同じ形式に変換
    Returns:
        f_a : float | tuple[float] | list[float] | ndarray  質点Aに働く粘弾性力 [N]
        f_b : float | tuple[float] | list[float] | ndarray  質点Bに働く粘弾性力 [N]
        ※ いずれも pos_a と同じ形式（float / tuple / list / ndarray）で返します。
    Raises:
        TypeError  : 各引数の要素が実数でない場合、または各引数に数値・list・tuple・ndarray 以外を渡した場合
        ValueError : k, c が0以上の有限値でない場合、各引数が有限値でない場合、各引数に2次元以上の ndarray を渡した場合、
                     各引数の要素数が不一致・空の場合、pos_b と v_b の片方のみを指定した場合、または計算結果が桁あふれした場合
        ※ 値の妥当性チェックは calc_stiffness_force / calc_viscous_drag で行うため、エラーメッセージの先頭にはそれらの関数名が表示されます。
    Warns:
        UserWarning : pos_b, v_b を省略して pos_b_init のみを指定した場合（質点Bを pos_b_init の位置に固定された壁として計算します。pos_b_init がすべて0の場合は出しません）

    ========================================================
    Copyright (c) 2026 1Dsim.com

    本関数の利用・改変・再配布の条件は、
    1Dsim.com の「プライバシーポリシー・免責事項」
    （https://www.1dsim.com/privacy-policy/）に従います。
    このCopyright表記は削除しないでください。
    ========================================================
    """

    # エラーメッセージの先頭に表示する関数名を定義
    func_name = "calc_kelvin_voigt_force"

    # --- 1. 入力値の整合性チェック（pos_b/v_bの指定方法・要素数・固定壁の位置） ---
    # pos_bとv_bの片方だけが省略され、「固定壁なのに動いている」などの矛盾状態になるのを防ぎます。
    if (pos_b is None) != (v_b is None):
        raise ValueError(f"{func_name}: pos_bとv_bは両方指定するか、両方省略する必要があります。")

    # 内部関数（_to_array_1dsim_com）が定義されている場合は、ばね側と粘性側の要素数が一致しているか確認します。
    # （calc_stiffness_force と calc_viscous_drag は、それぞれの引数の中でしか要素数を確認しないため）
    if "_to_array_1dsim_com" in globals():
        dims = {
            "k": len(_to_array_1dsim_com(k, func_name)),
            "c": len(_to_array_1dsim_com(c, func_name)),
            "pos_a": len(_to_array_1dsim_com(pos_a, func_name)),
            "v_a": len(_to_array_1dsim_com(v_a, func_name)),
        }

        # pos_b・v_b・pos_a_init・pos_b_init・f_initは省略可能なため、
        # 指定されているものだけを要素数チェックの対象に加えます。
        if pos_b is not None:
            dims["pos_b"] = len(_to_array_1dsim_com(pos_b, func_name))
            dims["v_b"] = len(_to_array_1dsim_com(v_b, func_name))
        if pos_a_init is not None:
            dims["pos_a_init"] = len(_to_array_1dsim_com(pos_a_init, func_name))
        if pos_b_init is not None:
            dims["pos_b_init"] = len(_to_array_1dsim_com(pos_b_init, func_name))
        if f_init is not None:
            dims["f_init"] = len(_to_array_1dsim_com(f_init, func_name))

        # 各入力の要素数（成分の数）が一致しているか確認します。
        if len(set(dims.values())) != 1:
            raise ValueError(f"{func_name}: 入力値の要素数が一致していません。各引数の要素数: {dims}")

    # pos_b, v_b を省略して pos_b_init のみが指定された場合、そのままでは
    # 「壁が pos_b_init から 0 へ動いた」と扱われ、意図しないばね力が生じます。
    # これを防ぐため、壁は初期位置から動いていない（pos_b = pos_b_init）とみなします。
    is_fixed_wall_at_init = (pos_b is None and pos_b_init is not None)
    if is_fixed_wall_at_init:
        pos_b = pos_b_init

    # --- 2. ばね成分・粘性成分の計算 ---
    # ばね成分：各引数を calc_stiffness_force に渡し、初期位置や初期荷重も含めたばね力を計算します。
    f_a_spring, f_b_spring = calc_stiffness_force(k, pos_a, pos_b, pos_a_init, pos_b_init, f_init)

    # 粘性成分：v_bを省略した場合は、calc_viscous_drag側でも固定壁（速度0）として扱われるため、そのまま渡します。
    f_a_damp, f_b_damp = calc_viscous_drag(c, v_a, v_b)

    # --- 3. ばね成分と粘性成分の合算 ---
    # 内部関数（_to_array_1dsim_com）が定義されている場合は、ばね成分（pos_aの形式）と粘性成分（v_aの形式）を
    # いったん1次元ndarrayにそろえてから足し合わせます。
    if "_to_array_1dsim_com" in globals():
        f_a_spring_calc = _to_array_1dsim_com(f_a_spring, func_name)
        f_b_spring_calc = _to_array_1dsim_com(f_b_spring, func_name)
        f_a_damp_calc   = _to_array_1dsim_com(f_a_damp, func_name)
        f_b_damp_calc   = _to_array_1dsim_com(f_b_damp, func_name)
    else:
        # _to_array_1dsim_com を使用しない場合は、変換を行わず、数値またはndarrayとしてそのまま足し合わせます。
        # （+ 0.0 により、整数型のndarrayが渡された場合でもfloat型に変換し、常に浮動小数点数で計算します）
        f_a_spring_calc = f_a_spring + 0.0
        f_b_spring_calc = f_b_spring + 0.0
        f_a_damp_calc   = f_a_damp + 0.0
        f_b_damp_calc   = f_b_damp + 0.0

    # 質点A・Bに働く粘弾性力（ばね成分 + 粘性成分）を計算します [N]
    # （+ 0.0 により、力が0のとき -0.0 ではなく 0.0 を返します）
    f_a = f_a_spring_calc + f_a_damp_calc + 0.0
    f_b = f_b_spring_calc + f_b_damp_calc + 0.0

    # --- 4. 出力形式の調整 ---
    # 出力形式は pos_a の形式（型）のみで決定します。
    # 内部関数（_format_output_1dsim_com）が定義されている場合は、桁あふれを確認し、pos_a と同じ形式に変換します。
    if "_format_output_1dsim_com" in globals():
        f_a, f_b = _format_output_1dsim_com(pos_a, func_name, f_a, f_b)

    # すべてのエラーチェックが完了した後に、
    # pos_b, v_b を省略して pos_b_init（0以外）のみが指定されていた場合は警告を出します。
    # （pos_b_init がすべて0の場合は置き換えても結果が変わらないため、警告は出しません）
    if "warnings" in globals() and is_fixed_wall_at_init:
        # _to_array_1dsim_com の有無にかかわらず判定できるよう、pos_b_init の各要素をリストに取り出します。
        raw_b_init = pos_b_init.tolist() if hasattr(pos_b_init, "tolist") else pos_b_init
        b_init_vals = raw_b_init if isinstance(raw_b_init, (list, tuple)) else [raw_b_init]
        if any(x != 0 for x in b_init_vals):
            warnings.warn(
                f"{func_name}: pos_b, v_bを省略してpos_b_initのみが指定されたため、"
                "質点Bをpos_b_initの位置に固定された壁（pos_b = pos_b_init）として計算します。",
                stacklevel=2,  # 警告の発生箇所として、この関数の呼び出し元の行を表示します
            )

    # 計算結果を返します（_format_output_1dsim_com を使用しない場合は、数値またはndarrayのまま返ります）。
    return f_a, f_b

# ===== 1次元のテスト（スカラー入力） =====
res_1d = calc_kelvin_voigt_force(10.0, 2.0, 0.5, 0.2)
print("1D Result:", res_1d, "\n")

# ===== 多次元のテスト（ndarray入力） =====
res_nd = calc_kelvin_voigt_force(
    np.array([10.0, 20.0]), np.array([1.0, 2.0]),
    np.array([1.0, -1.0]), np.array([0.5, 0.0]),
    np.array([2.0, 0.0]), np.array([-0.5, 0.5])
)
print("ndarray Result:", res_nd)