以下は、lsq_latt2.py プログラムのソースコードを解析し、Sphinx（MyST）で問題なくビルド可能なMarkdownドキュメントとして記述したものです。

# ``lsq_latt2.py`` プログラム解説

## 概要

``lsq_latt2.py`` は、X線回折データなどの回折データから結晶の格子定数を導出するために設計されたPythonスクリプトです。このプログラムは、観測データの前処理、重み付け、重み付き最小二乗フィッティング、そして最終的な格子定数（逆格子および実格子）の導出を行います。様々な格子系に対応しており、計算の精度を向上させるための外れ値の自動除去機能も備えています。

入力ファイルは特定のフォーマットに従い、複数の計算ジョブを連続して処理できます。各ジョブでは、格子系、補正フラグ、重み付けモード、回折位置の変換モードなどが設定され、それに従って観測データが処理されます。計算結果は整形されたレポートファイルとして出力され、各ジョブのプロットも生成されます。

## 機能

*   **格子定数計算**: X線回折データなどの回折データに基づき、最小二乗法を用いて結晶の格子定数を計算します。
*   **多様な格子系への対応**: 三斜晶、単斜晶、直方晶、正方晶、立方晶、菱面体晶（三方晶）、六方晶の各格子系に対応しています。
*   **データ前処理**: 生の回折位置データを指定されたモード（例: ``sin(theta)``、``lambda/(2*d)`` 型変換など）に従って変換します。
*   **重み付け**: 観測データに重みを適用し、重み付き最小二乗法を実行します。重み付けモードは、一定重み、シグマ値に基づく重み、ユーザー指定重みなどに対応しています。
*   **補正項**: 特定の系統誤差（例: 試料位置ずれ、ローレンツ・偏光因子など）を考慮するための補正項を最小二乗モデルに含めることができます。
*   **外れ値除去**: 最小二乗フィッティング中に、残差に基づいて統計的に有意な外れ値を自動的に検出し、除去する機能を備えています。
*   **誤差伝播**: 最小二乗法で得られたパラメータの標準偏差から、導出された格子定数の標準偏差を誤差伝播法を用いて計算します。
*   **詳細なレポート出力**: 計算された格子定数、フィットされたパラメータ、補正モデル、および各反射の詳細情報（観測値、計算値、残差、重みなど）を含むテキストレポートを生成します。
*   **結果のプロット**: 各計算ジョブについて、Q (散乱ベクトル) 対 1/d (逆格子距離) のプロットを生成し、観測値と計算値、および誤差棒を表示します。

## 非標準ライブラリ

このプログラムは以下の非標準（Python標準ライブラリ以外）のライブラリを使用しています。

*   ``matplotlib``: データのプロットと可視化に使用されます。
*   ``numpy``: 数値計算、特に配列操作、行列計算、最小二乗法の求解に使用されます。

なお、``argparse``, ``math``, ``sys``, ``dataclasses``, ``pathlib``, ``typing`` はPythonの標準ライブラリに含まれるため、別途インストールする必要はありません。

## インストール

必要な非標準ライブラリは ``pip`` を使ってインストールできます。

```bash
pip install numpy matplotlib

使用方法

コマンドライン引数

lsq_latt2.py は以下のコマンドライン引数をサポートしています。

  • input_file:

    • 型: Path

    • 説明: 入力データファイルへのパスを指定します。

  • output_file:

    • 型: Path

    • 説明: 出力レポートファイルへのパスを指定します。指定がない場合、入力ファイル名に .out 拡張子を付加したファイル名が使用されます。

    • オプション: nargs='?'

  • --silent:

    • 型: store_true

    • 説明: プログラムの開始バナーと終了メッセージのコンソール出力を抑制します。

    • オプション: フラグ

  • --no-show:

    • 型: store_true

    • 説明: プロットファイルを保存しますが、プロットウィンドウを自動的に開きません。

    • オプション: フラグ

使用例:

python lsq_latt2.py data.txt
python lsq_latt2.py my_experiment.dat --output_file results.txt --no-show

入力ファイル形式

入力ファイルは、一連の計算ジョブで構成されます。各ジョブは、タイトル、設定パラメータ、波長、装置半径、そして複数の回折観測データを含みます。

入力ファイルの全体構造

入力ファイルは、以下の要素の繰り返しで構成されます。

  1. ジョブタイトル行

  2. ジョブ設定パラメータ行

  3. 初期波長と装置半径行

  4. 反射データブロック:

    • 複数の反射データ行

    • (位置モード3の場合) ゼロシフト行

    • (重みモード0, 1, 3, 4の場合) 重みまたはシグマ値行

    • 特殊な更新行: 2000 wavelength radius

    • 特殊な終了行: 1000 ...

  5. ジョブ終了コード行 (0 または 1)

ファイル終端までこの繰り返しが続きます。

ジョブタイトル

ファイルの最初の行（または前のジョブの終了コードの後の行）は、そのジョブのタイトルです。 空行はスキップされます。

例:

My First Experiment

ジョブ設定パラメータ行

ジョブタイトル行の次の行は、計算ジョブの設定パラメータです。スペース区切りで9個以上の整数値が並びます。

フォーマット: lattice_system C1_flag C2_flag C3_flag C4_flag C5_flag weight_mode position_mode profile_mode [その他]

  • lattice_system (整数): 格子系の種類。

    • 1: 三斜晶 (Triclinic)

    • 2: 単斜晶 (Monoclinic, unique b axis)

    • 3: 単斜晶 (Monoclinic, unique c axis)

    • 4: 直方晶 (Orthorhombic)

    • 5: 正方晶 (Tetragonal)

    • 6: 立方晶 (Cubic)

    • 7: 菱面体晶 (Rhombohedral / Trigonal)

    • 8: 六方晶 (Hexagonal)

  • C1_flag, C2_flag, C3_flag, C4_flag, C5_flag (整数): 補正項の適用フラグ。

    • 0: 補正なし

    • 1: 補正あり

  • weight_mode (整数): 重み付けモード。

    • 0: sigma(theta) ベース重み

    • 1: sigma(2theta) ベース重み

    • 2: 一定重み (1.0)

    • 3: ユーザー指定重み

    • 4: 1/sigma^2

  • position_mode (整数): 回折位置変換モード。

    • 1: 角度 (度数) から sin(theta) へ変換

    • 2: 生の位置値をそのまま使用

    • 3: ゼロシフト補正付き sin(theta)

    • 4: sin(theta/2)

    • 5: lambda/(2*d) 型変換

    • その他 (例: 6): 幾何条件依存の変換

  • profile_mode (整数): プロファイル分析モード。このプログラム内では使用されていません。

例:

1 1 0 0 0 0 2 1 0

(三斜晶、C1補正有効、一定重み、sin(theta) 変換、プロファイルモード0)

波長/半径行

ジョブ設定パラメータ行の次の行は、初期波長と装置半径です。スペース区切りで2つの浮動小数点数が並びます。

フォーマット: wavelength radius

  • wavelength (浮動小数点数): X線の波長。

  • radius (浮動小数点数): 測定装置の半径 (位置モード3などで使用)。

例:

1.54059 100.0

反射データ行

波長/半径行の次から、反射データのブロックが始まります。各反射は通常、ミラー指数と生の回折位置データを含む1行で指定されます。

フォーマット: h k l raw_position

  • h, k, l (整数): ミラー指数。

  • raw_position (浮動小数点数): 生の回折位置データ。単位は position_mode によって異なります（度数、mmなど）。

特殊な反射データ行:

  • 波長/半径更新行: h の値が 2000 の場合、波長と半径を更新するための行として解釈されます。この行の次の行に、新しい波長と半径が指定されます。 フォーマット: 2000 new_wavelength new_radius 例:

    2000
1.54059 100.0

    (1000 より小さい値の h, k, l を伴う通常の反射データ行の直後に続くべきです。)

  • ジョブ終了シグナル: h の値が 1000 以上で、かつ 2000 でない場合、現在のジョブの反射データブロックの終了を示します。この行自体は反射データとして扱われません。

例:

1 0 0 10.234
1 1 0 14.567

ゼロシフト行（位置モード3の場合）

ジョブ設定の position_mode が 3 の場合、各反射データ行の直後にゼロシフト値を含む行が続きます。

フォーマット: zero_shift

  • zero_shift (浮動小数点数): ゼロ点シフト値。

例 (反射データ行の次に続く):

1 0 0 10.234
0.012

重み/シグマ行（重みモード0, 1, 3, 4の場合）

ジョブ設定の weight_mode が 0, 1, 3, 4 のいずれかの場合、各反射データ行（および位置モード3の場合はゼロシフト行）の直後に、重みまたはシグマ値を含む行が続きます。

フォーマット: sigma_or_weight

  • sigma_or_weight (浮動小数点数):

    • 重みモード 0 または 1 または 4 の場合: 回折角度の標準偏差。

    • 重みモード 3 の場合: 直接指定される重み値。

例 (反射データ行の次に続く、位置モード3でなければ):

1 0 0 10.234
0.005

例 (反射データ行の次に続く、位置モード3の場合):

1 0 0 10.234
0.012  (ゼロシフト値)
0.005  (重み/シグマ値)

ジョブ終了コード行

各ジョブの反射データブロックの後に、そのジョブの終了を示す整数値を含む行が続きます。

フォーマット: job_code

  • job_code (整数):

    • 0: 全てのジョブが終了したことを示します。プログラムが終了します。

    • 1: 次のジョブが続くことを示します。プログラムは次のジョブの解析に進みます。

例:

1

(このジョブの後に別のジョブが続く)

出力ファイル形式

レポートファイル

指定された出力ファイル、または自動生成されたファイル (例: input.out) には、各ジョブの詳細な計算レポートがテキスト形式で書き込まれます。レポートは以下のようなセクションで構成されます。

  • ジョブ概要: 格子系、参照波長、使用されたモード、反射数、自由度、重み付き残差二乗和、RMS残差。

  • 実格子定数: a, b, c, alpha, beta, gamma (度数とコサイン)、体積 V およびそれらの標準偏差。

  • 逆格子定数: a*, b*, c*, alpha*, beta*, gamma* (度数とコサイン)、体積 V* およびそれらの標準偏差。

  • 格子フィットパラメータ: 最小二乗法でフィットされた各格子パラメータ（例: a*^2, 2 b* c* cos(alpha*) など）の値と標準偏差。

  • 補正モデル: 有効化された各補正項 (C1-C5) の式、推定される物理的意味（英語と日本語）、およびフィットされた係数とその標準偏差。

  • 反射一覧: 最小二乗法に使用された各反射のデータ。通し番号、hkl、観測 1/d、計算 1/d、生の回折位置、計算 2θ、y 空間での残差、2θ 空間での残差、重み。

  • 除外された反射: 外れ値として除外された反射のリスト。通し番号、hkl、生の回折位置、除外理由。

プロットファイル

各ジョブに対して、プロットファイルがPNG形式で保存されます。ファイル名は [output_fileのstem]_job[ジョブ番号]_title.png の形式で自動生成されます。 例: results_job01_My_First_Experiment.png

プロットには以下の情報が表示されます。

  • x軸: 散乱ベクトル Q。

  • y軸: 逆格子距離 1/d。

  • 凡例:

    • 理論曲線: 1/d = Q / 2π

    • 計算された 1/d (誤差棒付き): 最小二乗法でフィットされた 1/d 値とその標準偏差。

    • 観測された 1/d: 入力データから計算された 1/d 値。

  • 右y軸: σ(1/d) (1/d の標準偏差)。

  • タイトル: ジョブのタイトル。

データクラス

Observation

一つの回折観測データを保持するデータクラスです。

  • serial_index: int

    • 観測データの通し番号。

  • h: int

    • ミラー指数 h。

  • k: int

    • ミラー指数 k。

  • l: int

    • ミラー指数 l。

  • raw_position: float

    • 生の回折位置データ。

  • transformed_position: float

    • 変換された回折位置データ（最小二乗法に使用）。通常は 1/d^2 に相当する値。

  • weight: float

    • この観測データに対する重み。

  • design_row: np.ndarray

    • この観測データに対応するデザイン行列の行。

  • wavelength: float

    • この観測で使用されたX線の波長。

LeastSquaresResult

最小二乗法の結果を保持するデータクラスです。

  • parameters: np.ndarray

    • フィットされたパラメータの配列。

  • parameter_sigma: np.ndarray

    • 各パラメータの標準偏差の配列。

  • fitted_y: np.ndarray

    • フィットモデルによって計算された y 値の配列（1/d^2 に対応）。

  • residuals: np.ndarray

    • 実測 y 値とフィットされた y 値の残差の配列。

  • observation_sigma: np.ndarray

    • 各観測の標準偏差の配列。

  • kept_observations: list[Observation]

    • 最小二乗法に使用された観測データのリスト。

  • covariance: np.ndarray

    • パラメータの共分散行列。

  • rejected_observations: list[Observation]

    • 外れ値として除外された観測データのリスト。デフォルトは空のリスト。

CellConstants

格子定数（逆格子および実格子）を保持するデータクラスです。

  • reciprocal_lengths: np.ndarray

    • 逆格子の長さ (a*, b*, c*) の配列。デフォルトは全て 0.0。

  • reciprocal_length_sigma: np.ndarray

    • 逆格子の長さの標準偏差の配列。デフォルトは全て 0.0。

  • reciprocal_angles_deg: np.ndarray

    • 逆格子の角度 (alpha*, beta*, gamma*) の配列（度数）。デフォルトは全て 90.0。

  • reciprocal_angle_sigma_deg: np.ndarray

    • 逆格子の角度の標準偏差の配列（度数）。デフォルトは全て 0.0。

  • reciprocal_cosines: np.ndarray

    • 逆格子の角度のコサイン (cos(alpha*), cos(beta*), cos(gamma*)) の配列。デフォルトは全て 0.0。

  • reciprocal_cosine_sigma: np.ndarray

    • 逆格子の角度のコサインの標準偏差の配列。デフォルトは全て 0.0。

  • reciprocal_volume: float

    • 逆格子の体積。デフォルトは 0.0。

  • reciprocal_volume_sigma: float

    • 逆格子の体積の標準偏差。デフォルトは 0.0。

  • direct_lengths: np.ndarray

    • 実格子の長さ (a, b, c) の配列。デフォルトは全て 0.0。

  • direct_length_sigma: np.ndarray

    • 実格子の長さの標準偏差の配列。デフォルトは全て 0.0。

  • direct_angles_deg: np.ndarray

    • 実格子の角度 (alpha, beta, gamma) の配列（度数）。デフォルトは全て 90.0。

  • direct_angle_sigma_deg: np.ndarray

    • 実格子の角度の標準偏差の配列（度数）。デフォルトは全て 0.0。

  • direct_cosines: np.ndarray

    • 実格子の角度のコサイン (cos(alpha), cos(beta), cos(gamma)) の配列。デフォルトは全て 0.0。

  • direct_cosine_sigma: np.ndarray

    • 実格子の角度のコサインの標準偏差の配列。デフォルトは全て 0.0。

  • direct_volume: float

    • 実格子の体積。デフォルトは 0.0。

  • direct_volume_sigma: float

    • 実格子の体積の標準偏差。デフォルトは 0.0。

JobConfig

格子定数計算ジョブの設定を保持するデータクラスです。

  • lattice_system: int

    • 格子系の種類を示す整数コード（例: 1=三斜晶, 4=直方晶など）。

  • correction_flags: list[int]

    • 5つの補正フラグのリスト。各フラグは特定の補正の適用有無を示す（0:なし, 1:あり）。

  • weight_mode: int

    • 観測データに適用する重み付けモード。

  • position_mode: int

    • 生の回折位置データを変換するモード。

  • profile_mode: int

    • プロファイル分析モード（このコードでは未使用）。

JobData

一つの格子定数計算ジョブに関する全てのデータを保持するデータクラスです。

  • title: str

    • ジョブのタイトル。

  • config: JobConfig

    • このジョブの設定情報。

  • observations: list[Observation]

    • このジョブで処理される観測データのリスト。

  • reference_wavelength: float

    • このジョブの参照波長。

例外クラス

LatticeCalculationError

格子定数計算におけるエラーを示す例外です。

この例外は、入力データが無効な場合、または最小二乗問題の解決中に計算上の問題（例: 特異行列）が発生した場合に発生します。

ヘルパークラス

LineReader

テキストストリームから行を読み込むためのヘルパークラスです。

このクラスは、テキストファイルを一行ずつ読み込み、現在の行番号を追跡し、ファイルの予期せぬ終端を検出する機能を提供します。空行でない行を確実に読み込むためのメソッドも備えています。

  • __init__(self, stream: TextIO)

    • LineReader の新しいインスタンスを初期化します。

    • stream: TextIO: 読み込むテキストストリーム。

  • read_line(self, *, allow_eof: bool = False) -> str | None

    • ストリームから一行を読み込みます。

    • ファイルの終端に達した場合、allow_eof が True であれば None を返します。False の場合は EOFError を発生させます。読み込んだ行末の改行文字は除去されます。

    • allow_eof: bool: True の場合、ファイルの終端に達しても EOFError を発生させず、None を返します。デフォルトは False。

    • 戻り値: str | None: 読み込まれた行（末尾の改行文字は除去）。EOFで allow_eof が True の場合は None。

    • 発生しうる例外: EOFError: allow_eof が False のときに、ファイルの予期せぬ終端に達した場合。

  • read_nonempty_line(self) -> str

    • ストリームから空でない行を読み込みます。

    • 空行をスキップし、最初に見つかった空でない行を返します。ファイルの終端に達した場合は EOFError を発生させます。

    • 戻り値: str: 読み込まれた空でない行（末尾の改行文字は除去）。

    • 発生しうる例外: EOFError: ファイルの終端に達した場合。

SigmaPropagator

誤差伝播の計算を行うクラスです。

このクラスは、複数の入力値とそれぞれの標準偏差が与えられたときに、それらの入力値に依存する関数が返す値とその標準偏差を、数値微分（摂動法）を用いて計算します。

  • __init__(self, values: Iterable[float], sigmas: Iterable[float])

    • SigmaPropagator の新しいインスタンスを初期化します。

    • values: Iterable[float]: 関数に渡す入力値のシーケンス。

    • sigmas: Iterable[float]: 各入力値の標準偏差のシーケンス。

  • evaluate(self, func: Callable[[np.ndarray], float]) -> tuple[float, float]

    • 与えられた関数を評価し、その値と標準偏差を計算します。

    • 入力値の各要素を微小に摂動させ、その変化から数値微分を計算します。これらの数値微分と各入力値の標準偏差を用いて、誤差伝播の法則により関数の結果の標準偏差を導出します。

    • func: Callable[[np.ndarray], float]: 誤差伝播を計算する関数。NumPy配列を受け取り、単一の float を返す必要があります。

    • 戻り値: tuple[float, float]: 関数の評価値と、その標準偏差のタプル。

TeeStream

ファイルとコンソールへ同時に書き込む簡易ストリームです。

このクラスは、複数の TextIO ストリーム（例: ファイルオブジェクトと sys.stdout）をラップし、write メソッドが呼び出された際に、登録された全てのストリームにテキストを転送します。これにより、計算結果をファイルに保存しつつ、同時にコンソールにも表示することができます。

  • __init__(self, *streams: TextIO)

    • TeeStream の新しいインスタンスを初期化します。

    • streams: TextIO: 書き込み先の TextIO ストリームを可変長引数で指定。

  • write(self, text: str) -> int

    • 全ての登録ストリームにテキストを書き込みます。

    • write メソッドが呼び出されると、登録された全ての TextIO ストリームに対して引数 text を書き込みます。また、各ストリームの flush メソッドも呼び出されます。

    • text: str: 書き込むテキスト。

    • 戻り値: int: 書き込まれた文字数 (入力 text の長さ)。

  • flush(self) -> None

    • 全ての登録ストリームをフラッシュします。

    • 登録された全ての TextIO ストリームの flush メソッドを呼び出し、バッファリングされたデータを強制的に出力します。

関数

clamp_unit()

値を -1.0 から 1.0 の範囲にクランプします。

主に math.acos の引数が有効な範囲内にあることを保証するために使用されます。入力値が -1.0 未満であれば -1.0 に、1.0 を超えれば 1.0 に丸めます。

  • x: float: クランプする値。

  • 戻り値: float: -1.0 から 1.0 の範囲にクランプされた値。

safe_sqrt()

非負の引数に対して安全に平方根を計算します。

math.sqrt は負の引数に対して ValueError を発生させますが、この関数は引数が負の場合でも 0.0 として処理することでエラーを回避し、常に非負の値を返します。

  • x: float: 平方根を計算する値。

  • 戻り値: float: x の平方根。x が負の場合は 0.0.

acos_deg()

コサイン値から角度を度数で計算します。

入力されたコサイン値が有効な範囲 [-1, 1] 内にあることを保証するために、clamp_unit() 関数で値をクランプしてから math.acos を適用します。結果はラジアンから度数に変換されます。

  • cosine: float: コサイン値。

  • 戻り値: float: 度数で表された角度。

reciprocal_volume_from_cosines()

逆格子のコサイン値から逆格子単位胞体積の計算に用いる係数を計算します。

逆格子単位胞体積 V* は、V* = a*b*c* * sqrt(1 - cos^2(alpha*) - cos^2(beta*) - cos^2(gamma*) + 2cos(alpha*)cos(beta*)cos(gamma*)) で与えられます。この関数は、上記式の sqrt(...) の部分を計算します。

  • cosines: np.ndarray: 逆格子の角度のコサイン値 (cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。

  • 戻り値: float: 逆格子単位胞体積計算に使用される平方根の値。

direct_cosine_from_reciprocal()

逆格子のコサイン値から実格子のコサイン値を計算します。

逆格子の角度 (alpha*, beta*, gamma*) のコサイン値から、対応する実格子の角度 (alpha, beta, gamma) のコサイン値を導出します。例えば、実格子の cos(alpha) は逆格子の cos(alpha*), cos(beta*), cos(gamma*) から計算されます。

  • values: np.ndarray: (cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。

  • 戻り値: float: 対応する実格子のコサイン値。

direct_length_from_reciprocal()

逆格子の長さとコサイン値から実格子の長さを計算します。

逆格子定数 (a*, cos_alpha*, cos_beta*, cos_gamma*) から、対応する実格子の長さ (a) を導出します。 例: 実格子の長さ a は、逆格子の a*, cos(alpha*), cos(beta*), cos(gamma*) から計算されます。

  • values: np.ndarray: (a_star, cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。

  • 戻り値: float: 対応する実格子の長さ。

direct_volume_from_reciprocal()

逆格子の長さとコサイン値から実格子の体積を計算します。

逆格子定数 (a*, b*, c*, cos_alpha*, cos_beta*, cos_gamma*) から実格子の体積を導出します。実格子の体積は逆格子の体積の逆数です。

  • values: np.ndarray: (a_star, b_star, c_star, cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。

  • 戻り値: float: 実格子の体積。

trigonal_direct_length()

三方晶系の逆格子定数から実格子の長さを計算します。

三方晶系において、逆格子の a* と cos(alpha*) から実格子の長さ a を導出します。等方的な格子定数を持つ三方晶系（菱面体晶）の特殊なケースに適用されます。

  • values: np.ndarray: (a_star, cos_alpha_star) を含むNumPy配列。

  • 戻り値: float: 三方晶系における実格子の長さ a。

trigonal_direct_cosine()

三方晶系の逆格子のコサイン値から実格子のコサイン値を計算します。

三方晶系において、逆格子の cos(alpha*) から実格子の cos(alpha) を導出します。

  • values: np.ndarray: (cos_alpha_star,) を含むNumPy配列。

  • 戻り値: float: 三方晶系における実格子のコサイン値 cos(alpha)。

trigonal_direct_volume()

三方晶系の逆格子定数から実格子の体積を計算します。

三方晶系において、逆格子の a* と cos(alpha*) から実格子の体積 V を導出します。

  • values: np.ndarray: (a_star, cos_alpha_star) を含むNumPy配列。

  • 戻り値: float: 三方晶系における実格子の体積 V。

base_design_terms()

回折指数 (h, k, l) からデザイン行列の基本的な項を計算します。

格子定数計算におけるデザイン行列を構築するための基本要素として、h^2, k^2, l^2, 2kl, 2lh, 2hk の6つの項を計算します。これらの項は、一般的な三斜晶系を含む様々な格子系のデザイン行列の基礎となります。

  • h: int: ミラー指数 h。

  • k: int: ミラー指数 k。

  • l: int: ミラー指数 l。

  • 戻り値: np.ndarray: 6つのデザイン行列項を含むNumPy配列。

lattice_design_row()

指定された格子系とミラー指数に基づいて、デザイン行列の行を構築します。

各格子系は、結晶学的に異なる格子定数の制約を持ちます。この関数は、その制約に基づいて、base_design_terms() から必要な項を選択・結合し、デザイン行列の1行を生成します。

格子系のコードと対応するパラメータ:

  • 1: 三斜晶 (Triclinic): a*^2, b*^2, c*^2, 2b*c*cos(alpha*), 2c*a*cos(beta*), 2a*b*cos(gamma*)

  • 2: 単斜晶 (Monoclinic, unique b): a*^2, b*^2, c*^2, 2c*a*cos(beta*)

  • 3: 単斜晶 (Monoclinic, unique c): a*^2, b*^2, c*^2, 2a*b*cos(gamma*)

  • 4: 直方晶 (Orthorhombic): a*^2, b*^2, c*^2

  • 5: 正方晶 (Tetragonal): a*^2 + b*^2, c*^2

  • 6: 立方晶 (Cubic): a*^2 + b*^2 + c*^2

  • 7: 菱面体晶 (Rhombohedral, trigonal): a*^2*(h^2+k^2+l^2), 2*a*^2*cos(alpha*)*(kl+lh+hk) に対応する係数

  • 8: 六方晶 (Hexagonal): a*^2*(h^2+k^2+hk), c*^2*l^2 に対応する係数

  • lattice_system: int: 格子系の種類を示す整数コード（1-8）。

  • h: int: ミラー指数 h。

  • k: int: ミラー指数 k。

  • l: int: ミラー指数 l。

  • 戻り値: np.ndarray: 指定された格子系のデザイン行列の行。

  • 発生しうる例外: LatticeCalculationError: 未対応の格子系が指定された場合。

transformed_position()

生の回折位置データを指定されたモードに基づいて変換します。

回折角度（2θ）や距離などの生の観測データを、最小二乗法で格子定数を決定するために適切な形式（例: sin(θ)^2 や 1/d^2 などに比例する値）に変換します。各モードは異なる回折実験のセットアップや測定方法に対応します。

モードの種類:

  • 1: sin(θ) （θは度数で与えられる2θの半分）

  • 2: 生の位置データそのまま

  • 3: ゼロ点シフトを考慮した sin(θ)（デバイシェラー法などで使用、raw_position は距離）

  • 4: sin(θ/2) （θは度数で与えられる2θの半分）

  • 5: λ/(2*d) または 1/d に関連する逆格子距離（例: raw_position が d または 2θ を表す場合）

  • その他 (例: 6): tan(θ) / sqrt(1 + tan^2(θ)) のような変換（raw_position が距離を表す場合）

  • raw_position: float: 生の回折位置データ。単位はモードによって異なる（度数、mmなど）。

  • mode: int: 位置変換モード（1-6）。

  • wavelength: float: 使用されるX線の波長。

  • radius: float: 測定装置の半径（位置モード3などで使用）。

  • zero_shift: float | None: 位置モード3で使用されるゼロ点シフト値。それ以外では None。

  • 戻り値: float: 変換された位置データ。

  • 発生しうる例外: LatticeCalculationError: 位置モード3でゼロシフト値が提供されない場合。

calc_weight()

指定された重みモードと変換された位置データに基づいて観測の重みを計算します。

最小二乗法では、各観測データの信頼度に応じて重みを設定することで、より信頼性の高いデータに大きな影響を与えることができます。この関数は、異なる重み付けスキーム（一定、測定誤差に基づくなど）を実装しています。

重みモードの種類:

  • 0: シグマ・シータに基づく重み。transformed_position が sin(theta) である場合、1/(sigma_theta^2 * 4 * sin^2(theta) * cos^2(theta)) に比例する項を含む。

  • 1: シグマ・ツーシータに基づく重み。transformed_position が sin(theta) である場合、1/(sigma_2theta^2 * 4 * sin^2(theta) * cos^2(theta)) に比例する項を含む。

  • 2: 一定の重み (1.0)。

  • 3: 直接指定された重み。

  • 4: シグマ値の逆二乗に基づく重み (1/sigma^2)。

  • weight_mode: int: 重み計算モード（0-4）。

  • transformed: float: 変換された回折位置データ（通常は sin(θ)）。

  • sigma_or_weight: float | None: モードによって、回折角度の標準偏差または直接指定される重み値。None の場合もある。

  • 戻り値: float: 計算された観測の重み。

  • 発生しうる例外: LatticeCalculationError: 重みモードに必要な追加データが提供されない場合。

correction_terms()

指定された補正フラグと変換された位置データに基づいて補正項を計算します。

回折データの測定には、様々な系統誤差が影響を与える可能性があります。この関数は、これらの誤差を補正するための特定の補正項を計算し、最小二乗法のデザイン行列に追加できるようにします。transformed は通常 sin(θ) に対応します。

  • flags: list[int]: 5つの補正フラグのリスト。各フラグは特定の補正の有無を示す（0:なし、1:あり）。

  • transformed: float: 変換された回折位置データ（通常は sin(θ)）。

  • 戻り値: list[float]: 計算された補正項のリスト。

read_numeric_line()

リーダーから一行を読み込み、指定された数の浮動小数点数に変換します。

入力ファイルから設定値やパラメータを読み込む際に使用されるヘルパー関数です。読み込んだ行をスペースで分割し、最初の expected 個の要素を浮動小数点数に変換します。

  • reader: LineReader: 入力ファイルを読み込むための LineReader インスタンス。

  • expected: int: 期待される数値の最小数。

  • 戻り値: list[float]: 読み込まれた数値のリスト。

  • 発生しうる例外: LatticeCalculationError: 期待されるよりも少ない数値が読み込まれた場合。

parse_jobs()

入力ストリームから複数の格子定数計算ジョブを解析します。

入力ファイルは、一連の計算ジョブで構成されます。各ジョブはタイトル、計算設定、および複数の回折観測データを含みます。この関数は、ストリームの終端に達するか、特定の終了コード（h=1000 で次のジョブ、job_code=0 で全体終了）が読み込まれるまで、ジョブを順次解析して JobData オブジェクトのリストを生成します。

  • stream: TextIO: 入力データを含むテキストストリーム。

  • 戻り値: list[JobData]: 解析された JobData オブジェクトのリスト。

  • 発生しうる例外: EOFError: ファイルの予期せぬ終端に達した場合。LatticeCalculationError: 入力ファイルの形式が無効な場合（例: ヘッダー情報不足、反射行のデータ不足）。

weighted_least_squares()

観測データに対して重み付き最小二乗法を実行します。

与えられた観測データに対し、重み付き最小二乗法を用いてモデルパラメータを推定します。このプロセスには、計算された残差に基づいて外れ値を自動的に検出し、除去する機能が含まれます。外れ値の除去は、安定したパラメータ推定値が得られるまで繰り返されます。

  • observations: list[Observation]: 最小二乗法に使用する観測データのリスト。

  • output: TextIO: 現在のところ、この関数からは直接出力されません。将来的なログ出力用。

  • 戻り値: LeastSquaresResult: 最小二乗法の計算結果、フィットされたパラメータ、残差、使用および除外された観測データを含むオブジェクト。

  • 発生しうる例外: LatticeCalculationError: 観測データが提供されない場合、またはデザイン行列が特異な場合。

derive_cell_constants()

最小二乗法で得られたパラメータから逆格子および実格子の定数を導出します。

最小二乗法によって決定された格子パラメータとそれらの標準偏差を受け取り、指定された格子系の結晶学的な関係に基づいて、最終的な逆格子定数 (a*, b*, c*, alpha*, beta*, gamma*, V*) および実格子定数 (a, b, c, alpha, beta, gamma, V) を計算します。各定数に対して誤差伝播法を用いて標準偏差も同時に計算されます。

  • lattice_system: int: 格子系の種類を示す整数コード（1-8）。

  • params: np.ndarray: 最小二乗法で得られたパラメータの配列。

  • sigma: np.ndarray: パラメータの標準偏差の配列。

  • 戻り値: CellConstants: 導出された格子定数を含むオブジェクト。

  • 発生しうる例外: LatticeCalculationError: 未対応の格子系が指定された場合。

make_default_output_path()

出力ファイルが未指定の場合、入力ファイルの拡張子を .out に置き換えたパスを生成します。

  • input_path: Path: 入力ファイルのパス。

  • 戻り値: Path: 生成された出力ファイルのパス。

lattice_system_name()

格子系のコードに対応する名前を日本語と英語で返します。

  • code: int: 格子系の整数コード。

  • 戻り値: str: 格子系の名前（例: "Triclinic / 三斜晶"）。

lattice_parameter_labels()

指定された格子系の最小二乗フィットパラメータに対応するラベルのリストを返します。

各格子系における結晶学的な制約に基づいて、フィットされたパラメータが何を表すかを示す適切なラベル（例: "a*^2", "2 c* a* cos(beta*)"）を提供します。

  • lattice_system: int: 格子系の種類を示す整数コード。

  • 戻り値: list[str]: パラメータラベルのリスト。

correction_metadata()

補正項に関するメタデータ（名前、式、説明）のリストを返します。

5種類の補正項それぞれについて、その名称、計算式、および物理的な意味に関する暫定的な解釈を英語と日本語で提供します。ユーザーが解釈を確認する責任があります。

  • 戻り値: list[dict[str, str]]: 補正項のメタデータ辞書のリスト。

format_value_sigma()

値と標準偏差を指定された書式で整形した文字列を生成します。

value と sigma を指定された桁数 (width) と精度 (prec) でフォーマットし、"値 ± 標準偏差" の形式で返します。

  • value: float: フォーマットする値。

  • sigma: float: フォーマットする標準偏差。

  • width: int: 数値部分の全体の幅。デフォルトは 14。

  • prec: int: 小数点以下の精度。デフォルトは 7。

  • 戻り値: str: 整形された文字列。

write_section_title()

レポートのセクションタイトルを英語と日本語で出力します。

  • output: TextIO: 書き込み先のストリーム。

  • title_en: str: 英語のタイトル。

  • title_ja: str: 日本語のタイトル。

parameter_blocks()

フィットされたパラメータを格子定数関連と補正項関連の2つのブロックに分割して返します。

最小二乗法でフィットされた全パラメータを、格子系に特有のパラメータと、有効化された補正項の係数とに分類し、それぞれのラベル、値、標準偏差のタプルリストとして返します。

  • job: JobData: 現在の計算ジョブのデータ。

  • fit: LeastSquaresResult: 最小二乗法のフィット結果。

  • 戻り値: tuple[list[tuple[str, float, float]], list[tuple[str, float, float]]]: 格子パラメータのリストと補正項係数のリストのタプル。

calc_fit_statistics()

フィットの統計情報（重み付き残差二乗和、RMS残差、自由度）を計算します。

  • fit: LeastSquaresResult: 最小二乗法のフィット結果。

  • 戻り値: tuple[float, float, int]: 重み付き残差二乗和、RMS残差、自由度のタプル。

mode_description_position()

位置変換モードのコードに対応する説明を返します。

  • mode: int: 位置変換モードの整数コード。

  • 戻り値: str: モードの説明文字列。

mode_description_weight()

重み付けモードのコードに対応する説明を返します。

  • mode: int: 重み付けモードの整数コード。

  • 戻り値: str: モードの説明文字列。

write_job_report()

単一の格子定数計算ジョブの結果を整形されたレポートとして出力します。

この関数は、ジョブの概要、実格子定数、逆格子定数、フィットされたパラメータ、適用された補正モデル、そして各反射の詳細なテーブル（使用されたものと除外されたもの）を含む包括的なレポートを output ストリームに書き出します。

  • output: TextIO: レポートを書き出すためのテキストストリーム。

  • job: JobData: 計算されたジョブの設定と観測データ。

  • fit: LeastSquaresResult: 最小二乗法のフィット結果。

  • cell: CellConstants: 導出された格子定数。

sanitize_stem()

ファイル名として安全な文字列を生成するため、テキストから無効な文字を除去します。

入力文字列から英数字、ハイフン、アンダースコア以外の文字をアンダースコアに置換し、前後のアンダースコアを削除して返します。空文字列になった場合は 'job' を返します。これにより、プラットフォーム間で互換性のあるファイル名を生成できます。

  • text: str: サニタイズする元の文字列。

  • 戻り値: str: サニタイズされた文字列。

plot_job_result()

各ジョブの Q-1/d プロットを生成し、画像ファイルとして保存します。

観測された 1/d と計算された 1/d を散乱ベクトル Q に対してプロットし、誤差棒も表示します。右軸には 1/d の標準偏差も表示されます。生成されたプロットはPNG画像として保存され、オプションで表示も可能です。

  • job: JobData: プロット対象のジョブデータ。

  • fit: LeastSquaresResult: プロットに使用するフィット結果。

  • output_file: Path: レポート出力ファイル。プロットのファイル名生成に使用される。

  • job_index: int: ジョブのインデックス（ファイル名に使用）。

  • show_plot: bool: True の場合、プロットウィンドウを表示する。デフォルトは True。

  • 戻り値: Path: 保存されたプロットファイルのパス。

build_argument_parser()

コマンドライン引数を解析するためのパーサーを構築します。

  • 戻り値: argparse.ArgumentParser: 構築された ArgumentParser オブジェクト。

main()

スクリプトのエントリポイントです。

コマンドライン引数を解析し、入力ファイルから格子定数計算ジョブを読み込みます。各ジョブに対して重み付き最小二乗法を実行し、格子定数を導出します。計算結果はレポートファイルとコンソールに表示され、各ジョブのプロットも生成・保存されます。

  • 戻り値: int: プログラムの終了コード。成功時は 0。