""" 最小二乗法を用いた格子定数計算スクリプト。 このスクリプトは、X線回折データなどの回折データから結晶の格子定数を導出します。 様々な格子系に対応し、観測データの前処理、重み付け、最小二乗フィッティング、 そして最終的な格子定数（逆格子および実格子）の導出を行います。 外れ値の自動除去機能も備えています。 詳細説明: 入力ファイルは特定のフォーマットに従い、複数の計算ジョブを連続して処理できます。 各ジョブでは、格子系、補正フラグ、重み付けモード、回折位置の変換モードなどが 設定され、それに従って観測データが処理されます。 重み付き最小二乗法によりパラメータが決定され、その結果から逆格子定数と実格子定数、 およびそれらの標準偏差が計算されます。計算結果は整形されたレポートファイルとして出力されます。 関連リンク: :doc:`lsq_latt2_usage` """ from __future__ import annotations import argparse import math from dataclasses import dataclass, field from pathlib import Path from typing import Callable, Iterable, TextIO import numpy as np DEG_PER_RAD = 180.0 / math.pi RAD_PER_DEG = math.pi / 180.0 ANGLE_EPS = 1e-12 PERTURB_SCALE = 1e-4 OUTLIER_SIGMA = 3.0 class LatticeCalculationError(RuntimeError): """ 格子定数計算におけるエラーを示す例外。 この例外は、入力データが無効な場合、または最小二乗問題の解決中に 計算上の問題（例: 特異行列）が発生した場合に発生します。 """ @dataclass(slots=True) class Observation: """ 一つの回折観測データを保持するデータクラス。 詳細説明: このクラスは、特定の回折面(hkl)に対応する生の回折位置、変換された位置、 計算で使用される重み、デザイン行列の行、および観測時の波長などの 詳細な情報を含みます。 :param serial_index: int: 観測データの通し番号。 :param h: int: ミラー指数 h。 :param k: int: ミラー指数 k。 :param l: int: ミラー指数 l。 :param raw_position: float: 生の回折位置データ。 :param transformed_position: float: 変換された回折位置データ（最小二乗法に使用）。 :param weight: float: この観測データに対する重み。 :param design_row: np.ndarray: この観測データに対応するデザイン行列の行。 :param wavelength: float: この観測で使用されたX線の波長。 """ serial_index: int h: int k: int l: int raw_position: float transformed_position: float weight: float design_row: np.ndarray wavelength: float @dataclass(slots=True) class LeastSquaresResult: """ 最小二乗法の結果を保持するデータクラス。 詳細説明: 最小二乗法によって計算された格子定数パラメータとその標準偏差、 観測データへのフィット値、残差、観測の標準偏差、および 計算に残された観測データのリストと共分散行列を格納します。 :param parameters: np.ndarray: 最小二乗法によって決定されたパラメータの配列。 :param parameter_sigma: np.ndarray: 各パラメータの標準偏差の配列。 :param fitted_y: np.ndarray: 各観測データに対してフィットされた予測値。 :param residuals: np.ndarray: 各観測データに対する残差 (観測値 - フィット値)。 :param observation_sigma: np.ndarray: 各観測データの標準偏差。 :param kept_observations: list[Observation]: 外れ値除去後に残った観測データのリスト。 :param covariance: np.ndarray: パラメータ間の共分散行列。 """ parameters: np.ndarray parameter_sigma: np.ndarray fitted_y: np.ndarray residuals: np.ndarray observation_sigma: np.ndarray kept_observations: list[Observation] covariance: np.ndarray @dataclass(slots=True) class CellConstants: """ 格子定数（逆格子および実格子）を保持するデータクラス。 詳細説明: 逆格子と実格子の両方について、長さ、角度（度数）、コサイン、 体積およびそれぞれの標準偏差を格納します。 デフォルト値としてゼロまたは標準的な値（角度90度）が設定されています。 :param reciprocal_lengths: np.ndarray: 逆格子の長さ (a*, b*, c*) の配列。デフォルトは全て0.0。 :param reciprocal_length_sigma: np.ndarray: 逆格子の長さの標準偏差の配列。デフォルトは全て0.0。 :param reciprocal_angles_deg: np.ndarray: 逆格子の角度 (alpha*, beta*, gamma*) の配列（度数）。デフォルトは全て90.0。 :param reciprocal_angle_sigma_deg: np.ndarray: 逆格子の角度の標準偏差の配列（度数）。デフォルトは全て0.0。 :param reciprocal_cosines: np.ndarray: 逆格子の角度のコサイン (cos(alpha*), cos(beta*), cos(gamma*)) の配列。デフォルトは全て0.0。 :param reciprocal_cosine_sigma: np.ndarray: 逆格子の角度のコサインの標準偏差の配列。デフォルトは全て0.0。 :param reciprocal_volume: float: 逆格子の体積。デフォルトは0.0。 :param reciprocal_volume_sigma: float: 逆格子の体積の標準偏差。デフォルトは0.0。 :param direct_lengths: np.ndarray: 実格子の長さ (a, b, c) の配列。デフォルトは全て0.0。 :param direct_length_sigma: np.ndarray: 実格子の長さの標準偏差の配列。デフォルトは全て0.0。 :param direct_angles_deg: np.ndarray: 実格子の角度 (alpha, beta, gamma) の配列（度数）。デフォルトは全て90.0。 :param direct_angle_sigma_deg: np.ndarray: 実格子の角度の標準偏差の配列（度数）。デフォルトは全て0.0。 :param direct_cosines: np.ndarray: 実格子の角度のコサイン (cos(alpha), cos(beta), cos(gamma)) の配列。デフォルトは全て0.0。 :param direct_cosine_sigma: np.ndarray: 実格子の角度のコサインの標準偏差の配列。デフォルトは全て0.0。 :param direct_volume: float: 実格子の体積。デフォルトは0.0。 :param direct_volume_sigma: float: 実格子の体積の標準偏差。デフォルトは0.0。 """ reciprocal_lengths: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) reciprocal_length_sigma: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) reciprocal_angles_deg: np.ndarray = field(default_factory=lambda: np.full(3, 90.0, dtype=float)) reciprocal_angle_sigma_deg: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) reciprocal_cosines: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) reciprocal_cosine_sigma: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) reciprocal_volume: float = 0.0 reciprocal_volume_sigma: float = 0.0 direct_lengths: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) direct_length_sigma: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) direct_angles_deg: np.ndarray = field(default_factory=lambda: np.full(3, 90.0, dtype=float)) direct_angle_sigma_deg: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) direct_cosines: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) direct_cosine_sigma: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float)) direct_volume: float = 0.0 direct_volume_sigma: float = 0.0 @dataclass(slots=True) class JobConfig: """ 格子定数計算ジョブの設定を保持するデータクラス。 詳細説明: このクラスは、計算に使用する格子系の種類、適用する補正の種類、 観測データの重み付け方法、回折位置の変換方法、およびプロファイルモードの 設定情報をカプセル化します。 :param lattice_system: int: 格子系の種類を示す整数コード（例: 1=三斜晶, 4=直方晶など）。 :param correction_flags: list[int]: 5つの補正フラグのリスト。各フラグは特定の補正の適用有無を示す。 :param weight_mode: int: 観測データに適用する重み付けモード。 :param position_mode: int: 生の回折位置データを変換するモード。 :param profile_mode: int: プロファイル分析モード（このコードでは未使用だが、入力形式の一部）。 """ lattice_system: int correction_flags: list[int] weight_mode: int position_mode: int profile_mode: int @dataclass(slots=True) class JobData: """ 一つの格子定数計算ジョブに関する全てのデータを保持するデータクラス。 詳細説明: ジョブのタイトル、設定、処理する観測データのリスト、および 参照として使用される波長を含みます。 :param title: str: ジョブのタイトル。 :param config: JobConfig: このジョブの設定情報。 :param observations: list[Observation]: このジョブで処理される観測データのリスト。 :param reference_wavelength: float: このジョブの参照波長。 """ title: str config: JobConfig observations: list[Observation] reference_wavelength: float class LineReader: """ テキストストリームから行を読み込むためのヘルパークラス。 詳細説明: このクラスは、テキストファイルを一行ずつ読み込み、現在の行番号を追跡し、 ファイルの予期せぬ終端を検出する機能を提供します。 空行でない行を確実に読み込むためのメソッドも備えています。 :param stream: TextIO: 読み込むテキストストリーム。 """ def __init__(self, stream: TextIO): """ LineReaderの新しいインスタンスを初期化する。 :param stream: TextIO: 読み込むテキストストリーム。 """ self.stream = stream self.line_number = 0 def read_line(self, *, allow_eof: bool = False) -> str | None: """ ストリームから一行を読み込む。 :param allow_eof: bool: Trueの場合、ファイルの終端に達してもEOFErrorを発生させず、Noneを返す。デフォルトはFalse。 :returns: str | None: 読み込まれた行（末尾の改行文字は除去）。EOFでallow_eofがTrueの場合はNone。 :raises EOFError: `allow_eof`がFalseのときに、ファイルの予期せぬ終端に達した場合。 """ line = self.stream.readline() if line == "": if allow_eof: return None raise EOFError("unexpected end of file") self.line_number += 1 return line.rstrip("\n") def read_nonempty_line(self) -> str: """ ストリームから空でない行を読み込む。 :returns: str: 読み込まれた空でない行（末尾の改行文字は除去）。 :raises EOFError: ファイルの終端に達した場合。 :raises AssertionError: 読み込まれた行がNoneであった場合 (内部的なガード)。 """ line = self.read_line() assert line is not None return line class SigmaPropagator: """ 誤差伝播の計算を行うクラス。 詳細説明: このクラスは、複数の入力値とそれぞれの標準偏差が与えられたときに、 それらの入力値に依存する関数が返す値とその標準偏差を、 数値微分（摂動法）を用いて計算します。 :param values: Iterable[float]: 関数に渡す入力値のシーケンス。 :param sigmas: Iterable[float]: 各入力値の標準偏差のシーケンス。 """ def __init__(self, values: Iterable[float], sigmas: Iterable[float]): """ SigmaPropagatorの新しいインスタンスを初期化する。 :param values: Iterable[float]: 関数に渡す入力値のシーケンス。 :param sigmas: Iterable[float]: 各入力値の標準偏差のシーケンス。 """ self.values = np.array(list(values), dtype=float) self.sigmas = np.array(list(sigmas), dtype=float) def evaluate(self, func: Callable[[np.ndarray], float]) -> tuple[float, float]: """ 与えられた関数を評価し、その値と標準偏差を計算する。 詳細説明: 入力値の各要素を微小に摂動させ、その変化から数値微分を計算します。 これらの数値微分と各入力値の標準偏差を用いて、誤差伝播の法則により 関数の結果の標準偏差を導出します。 :param func: Callable[[np.ndarray], float]: 誤差伝播を計算する関数。 NumPy配列を受け取り、単一のfloatを返す必要があります。 :returns: tuple[float, float]: 関数の評価値と、その標準偏差のタプル。 """ base = func(self.values.copy()) variance = 0.0 for idx, sigma in enumerate(self.sigmas): if abs(sigma) < ANGLE_EPS: continue shifted = self.values.copy() shifted[idx] += PERTURB_SCALE * sigma deriv_sigma = (func(shifted) - base) / PERTURB_SCALE variance += deriv_sigma * deriv_sigma return base, math.sqrt(max(variance, 0.0)) # ---------- elementary math helpers ---------- def clamp_unit(x: float) -> float: """ 値を -1.0 から 1.0 の範囲にクランプする。 詳細説明: 主に `math.acos` の引数が有効な範囲内にあることを保証するために使用されます。 入力値が -1.0 未満であれば -1.0 に、1.0 を超えれば 1.0 に丸めます。 :param x: float: クランプする値。 :returns: float: -1.0 から 1.0 の範囲にクランプされた値。 """ return max(-1.0, min(1.0, x)) def safe_sqrt(x: float) -> float: """ 非負の引数に対して安全に平方根を計算する。 詳細説明: `math.sqrt` は負の引数に対して `ValueError` を発生させますが、 この関数は引数が負の場合でも0.0として処理することでエラーを回避し、 常に非負の値を返します。 :param x: float: 平方根を計算する値。 :returns: float: xの平方根。xが負の場合は0.0。 """ return math.sqrt(max(0.0, x)) def acos_deg(cosine: float) -> float: """ コサイン値から角度を度数で計算する。 詳細説明: 入力されたコサイン値が有効な範囲 [-1, 1] 内にあることを保証するために、 `clamp_unit` 関数で値をクランプしてから `math.acos` を適用します。 結果はラジアンから度数に変換されます。 :param cosine: float: コサイン値。 :returns: float: 度数で表された角度。 """ return math.degrees(math.acos(clamp_unit(cosine))) def reciprocal_volume_from_cosines(cosines: np.ndarray) -> float: """ 逆格子のコサイン値から逆格子単位胞体積の計算に用いる係数を計算する。 詳細説明: 逆格子単位胞体積 V* は、V* = a*b*c* * sqrt(1 - cos^2(alpha*) - cos^2(beta*) - cos^2(gamma*) + 2cos(alpha*)cos(beta*)cos(gamma*)) で与えられます。この関数は、上記式の `sqrt(...)` の部分を計算します。 :param cosines: np.ndarray: 逆格子の角度のコサイン値 (cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。 :returns: float: 逆格子単位胞体積計算に使用される平方根の値。 """ ca, cb, cg = cosines term = 1.0 - ca * ca - cb * cb - cg * cg + 2.0 * ca * cb * cg return safe_sqrt(term) def direct_cosine_from_reciprocal(values: np.ndarray) -> float: """ 逆格子のコサイン値から実格子のコサイン値を計算する。 詳細説明: 逆格子の角度 (alpha*, beta*, gamma*) のコサイン値から、 対応する実格子の角度 (alpha, beta, gamma) のコサイン値を導出します。 例えば、実格子の cos(alpha) は逆格子の cos(alpha*), cos(beta*), cos(gamma*) から計算されます。 :param values: np.ndarray: (cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。 :returns: float: 対応する実格子のコサイン値。 """ cos_alpha_star, cos_beta_star, cos_gamma_star = values sin_beta_star = safe_sqrt(1.0 - cos_beta_star * cos_beta_star) sin_gamma_star = safe_sqrt(1.0 - cos_gamma_star * cos_gamma_star) denom = sin_beta_star * sin_gamma_star if denom <= ANGLE_EPS: return 0.0 return (cos_beta_star * cos_gamma_star - cos_alpha_star) / denom def direct_length_from_reciprocal(values: np.ndarray) -> float: """ 逆格子の長さとコサイン値から実格子の長さを計算する。 詳細説明: 逆格子定数 (a*, cos_alpha*, cos_beta*, cos_gamma*) から、 対応する実格子の長さ (a) を導出します。 例: 実格子の長さaは、逆格子のa*, cos(alpha*), cos(beta*), cos(gamma*) から計算されます。 :param values: np.ndarray: (a_star, cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。 :returns: float: 対応する実格子の長さ。 """ a_star, cos_alpha_star, cos_beta_star, cos_gamma_star = values if abs(a_star) <= ANGLE_EPS: return 0.0 sin_alpha_star = safe_sqrt(1.0 - cos_alpha_star * cos_alpha_star) v_star = reciprocal_volume_from_cosines(np.array([cos_alpha_star, cos_beta_star, cos_gamma_star])) if v_star <= ANGLE_EPS: return 0.0 return sin_alpha_star / (a_star * v_star) def direct_volume_from_reciprocal(values: np.ndarray) -> float: """ 逆格子の長さとコサイン値から実格子の体積を計算する。 詳細説明: 逆格子定数 (a*, b*, c*, cos_alpha*, cos_beta*, cos_gamma*) から 実格子の体積を導出します。実格子の体積は逆格子の体積の逆数です。 :param values: np.ndarray: (a_star, b_star, c_star, cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。 :returns: float: 実格子の体積。 """ a_star, b_star, c_star, cos_alpha_star, cos_beta_star, cos_gamma_star = values v_star = a_star * b_star * c_star * reciprocal_volume_from_cosines( np.array([cos_alpha_star, cos_beta_star, cos_gamma_star]) ) if v_star <= ANGLE_EPS: return 0.0 return 1.0 / v_star def trigonal_direct_length(values: np.ndarray) -> float: """ 三方晶系の逆格子定数から実格子の長さを計算する。 詳細説明: 三方晶系において、逆格子のa*とcos(alpha*)から実格子の長さaを導出します。 等方的な格子定数を持つ三方晶系（菱面体晶）の特殊なケースに適用されます。 :param values: np.ndarray: (a_star, cos_alpha_star) を含むNumPy配列。 :returns: float: 三方晶系における実格子の長さa。 """ a_star, cos_alpha_star = values if abs(a_star) <= ANGLE_EPS: return 0.0 numerator = 1.0 - cos_alpha_star * cos_alpha_star denominator = 1.0 - 3.0 * cos_alpha_star * cos_alpha_star + 2.0 * cos_alpha_star**3 if denominator <= ANGLE_EPS or numerator < 0.0: return 0.0 return safe_sqrt(numerator / denominator) / a_star def trigonal_direct_cosine(values: np.ndarray) -> float: """ 三方晶系の逆格子のコサイン値から実格子のコサイン値を計算する。 詳細説明: 三方晶系において、逆格子のcos(alpha*)から実格子のcos(alpha)を導出します。 :param values: np.ndarray: (cos_alpha_star,) を含むNumPy配列。 :returns: float: 三方晶系における実格子のコサイン値cos(alpha)。 """ (cos_alpha_star,) = values denominator = 1.0 - cos_alpha_star * cos_alpha_star if abs(denominator) <= ANGLE_EPS: return 0.0 return cos_alpha_star * (cos_alpha_star - 1.0) / denominator def trigonal_direct_volume(values: np.ndarray) -> float: """ 三方晶系の逆格子定数から実格子の体積を計算する。 詳細説明: 三方晶系において、逆格子のa*とcos(alpha*)から実格子の体積Vを導出します。 :param values: np.ndarray: (a_star, cos_alpha_star) を含むNumPy配列。 :returns: float: 三方晶系における実格子の体積V。 """ a_star, cos_alpha_star = values term = 1.0 - 3.0 * cos_alpha_star * cos_alpha_star + 2.0 * cos_alpha_star**3 if term <= 0.0: return 0.0 v_star = a_star**3 * safe_sqrt(term) if v_star <= ANGLE_EPS: return 0.0 return 1.0 / v_star # ---------- model construction ---------- def base_design_terms(h: int, k: int, l: int) -> np.ndarray: """ 回折指数 (h, k, l) からデザイン行列の基本的な項を計算する。 詳細説明: 格子定数計算におけるデザイン行列を構築するための基本要素として、 h^2, k^2, l^2, 2kl, 2lh, 2hk の6つの項を計算します。 これらの項は、一般的な三斜晶系を含む様々な格子系のデザイン行列の基礎となります。 :param h: int: ミラー指数h。 :param k: int: ミラー指数k。 :param l: int: ミラー指数l。 :returns: np.ndarray: 6つのデザイン行列項を含むNumPy配列。 """ return np.array( [ float(h * h), float(k * k), float(l * l), 2.0 * float(k * l), 2.0 * float(l * h), 2.0 * float(h * k), ], dtype=float, ) def lattice_design_row(lattice_system: int, h: int, k: int, l: int) -> np.ndarray: """ 指定された格子系とミラー指数に基づいて、デザイン行列の行を構築する。 詳細説明: 各格子系（三斜晶、単斜晶、直方晶、正方晶、六方晶、立方晶、菱面体晶）は、 結晶学的に異なる格子定数の制約を持ちます。この関数は、その制約に基づいて、 `base_design_terms` から必要な項を選択・結合し、デザイン行列の1行を生成します。 格子系のコードと対応するパラメータ: - 1: 三斜晶 (triclinic): a*, b*, c*, cos(alpha*), cos(beta*), cos(gamma*) - 2: 単斜晶 (monoclinic, unique b): a*, b*, c*, cos(beta*) - 3: 単斜晶 (monoclinic, unique c): a*, b*, c*, cos(gamma*) (※コードではgamma*として実装されているが、入力形式によってはalpha*を意味する場合もある) - 4: 直方晶 (orthorhombic): a*, b*, c* - 5: 正方晶 (tetragonal): a*, c* - 6: 立方晶 (cubic): a* - 7: 菱面体晶 (rhombohedral, trigonal): a*, cos(alpha*) (等軸菱面体) - 8: 六方晶 (hexagonal): a*, c* :param lattice_system: int: 格子系の種類を示す整数コード（1-8）。 :param h: int: ミラー指数h。 :param k: int: ミラー指数k。 :param l: int: ミラー指数l。 :returns: np.ndarray: 指定された格子系のデザイン行列の行。 :raises LatticeCalculationError: 未対応の格子系が指定された場合。 """ x = base_design_terms(h, k, l) if lattice_system == 1: # Triclinic return x.copy() if lattice_system == 2: # Monoclinic (unique b axis, beta*) return np.array([x[0], x[1], x[2], x[4]], dtype=float) # h^2, k^2, l^2, 2lh if lattice_system == 3: # Monoclinic (unique c axis, gamma*) or (unique a axis, alpha*) return np.array([x[0], x[1], x[2], x[5]], dtype=float) # h^2, k^2, l^2, 2hk if lattice_system == 4: # Orthorhombic return np.array([x[0], x[1], x[2]], dtype=float) # h^2, k^2, l^2 if lattice_system == 5: # Tetragonal return np.array([x[0] + x[1], x[2]], dtype=float) # h^2+k^2, l^2 if lattice_system == 6: # Cubic return np.array([x[0] + x[1] + x[2]], dtype=float) # h^2+k^2+l^2 if lattice_system == 7: # Rhombohedral (trigonal) return np.array([x[0] + x[1] + x[2], x[3] + x[4] + x[5]], dtype=float) # h^2+k^2+l^2, 2(kl+lh+hk) if lattice_system == 8: # Hexagonal return np.array([x[0] + x[1] + 0.5 * x[5], x[2]], dtype=float) # h^2+k^2+hk, l^2 raise LatticeCalculationError(f"unsupported lattice system: {lattice_system}") def transformed_position(raw_position: float, mode: int, wavelength: float, radius: float, zero_shift: float | None) -> float: """ 生の回折位置データを指定されたモードに基づいて変換する。 詳細説明: 回折角度（2θ）や距離などの生の観測データを、最小二乗法で格子定数を決定するために 適切な形式（例: sin(θ)^2 や 1/d^2 などに比例する値）に変換します。 各モードは異なる回折実験のセットアップや測定方法に対応します。 モードの種類: - 1: sin(θ) （θは度数） - 2: 生の位置データそのまま - 3: ゼロ点シフトを考慮したsin(θ)（デバイシェラー法などで使用） - 4: sin(θ/2) - 5: λ/(2*d) または 1/d に関連する逆格子距離（例: 2θが回折距離を表す場合） - その他: tan(θ) / (1 + tan^2(θ)) のような変換（通常はtan(θ)/sqrt(1+tan^2(θ)) = sin(θ)） :param raw_position: float: 生の回折位置データ。単位はモードによって異なる（度数、mmなど）。 :param mode: int: 位置変換モード（1-6）。 :param wavelength: float: 使用されるX線の波長。 :param radius: float: 測定装置の半径（位置モード3などで使用）。 :param zero_shift: float | None: 位置モード3で使用されるゼロ点シフト値。それ以外ではNone。 :returns: float: 変換された位置データ。 :raises LatticeCalculationError: 位置モード3でゼロシフト値が提供されない場合。 """ if mode == 1: return math.sin(raw_position * RAD_PER_DEG) if mode == 2: return raw_position if mode == 3: if zero_shift is None: raise LatticeCalculationError("position mode 3 requires a following zero-shift line") xyz = math.cos(raw_position / radius) abc = math.atan(zero_shift / radius) ghi = xyz * math.cos(abc) pqr = math.acos(clamp_unit(ghi)) return math.sin(0.5 * pqr) if mode == 4: return math.sin(0.5 * raw_position * RAD_PER_DEG) if mode == 5: # Assuming raw_position is 2θ in degrees, and calculating 1/d (or related) # 2d sin(θ) = nλ => 1/d = 2 sin(θ) / nλ # Here, it computes λ * 0.5 / raw_position which seems to imply raw_position is d (in some unit) # or related to 2θ in some other way. # Original code comment indicates 'transformed = lambda * 0.5 / D (or 2theta)' # Let's assume raw_position is effectively D in some form or a value from which 1/D is derived. return 0.0 if raw_position == 0.0 else wavelength * 0.5 / raw_position # Default behavior for other modes (e.g., mode 6 in some contexts) # This might correspond to tan(theta) / sqrt(1 + tan^2(theta)) = sin(theta) if raw_position/radius is tan(theta) abc = raw_position / radius xyz = 1.0 + abc * abc return 0.0 if xyz <= 0.0 else abc / math.sqrt(xyz) def calc_weight(weight_mode: int, transformed: float, sigma_or_weight: float | None) -> float: """ 指定された重みモードと変換された位置データに基づいて観測の重みを計算する。 詳細説明: 最小二乗法では、各観測データの信頼度に応じて重みを設定することで、 より信頼性の高いデータに大きな影響を与えることができます。 この関数は、異なる重み付けスキーム（一定、測定誤差に基づくなど）を実装しています。 重みモードの種類: - 0: シグマ・シータに基づく重み。 - 1: シグマ・ツーシータに基づく重み。 - 2: 一定の重み (1.0)。 - 3: 直接指定された重み。 - 4: シグマ値の逆二乗に基づく重み (1/sigma^2)。 :param weight_mode: int: 重み計算モード（0-4）。 :param transformed: float: 変換された回折位置データ。 :param sigma_or_weight: float | None: モードによって、回折角度の標準偏差または直接指定される重み値。Noneの場合もある。 :returns: float: 計算された観測の重み。 :raises LatticeCalculationError: 重みモードに必要な追加データが提供されない場合、または未対応のモードが指定された場合。 """ bas = transformed * transformed # transformed position squared if weight_mode == 0: if sigma_or_weight is None: raise LatticeCalculationError("weight mode 0 requires an additional sigma line") sigma_theta = 0.25 * sigma_or_weight sigma_theta = max(sigma_theta, 0.25) # Ensure minimum sigma_theta denom = bas * (1.0 - bas) # Original Fortran: 1/(sigma_theta^2 * 4 * transformed^2 * (1-transformed^2)) # This implementation: sigma_theta / denom, which might be a simplified constant, or a typo from 1/(sigma_theta^2 * X) # However, following the existing code logic: # The weight should be proportional to 1/variance, so usually 1/(sigma^2). # It's possible that `sigma_theta` here represents the numerator factor for 1/variance. # Given the existing code is `sigma_theta / denom` return 0.0 if abs(denom) <= ANGLE_EPS else sigma_theta / denom if weight_mode == 1: if sigma_or_weight is None: raise LatticeCalculationError("weight mode 1 requires an additional sigma line") sigma_two_theta = 2.0 * sigma_or_weight # sigma_two_theta is related to two_theta denom = bas * (1.0 - bas) * sigma_two_theta * sigma_two_theta # Original Fortran: 1/(sigma_2theta^2 * 4 * transformed^2 * (1-transformed^2)) # This implementation: 0.25 / denom, which is consistent with the Fortran `1 / (4 * denom)` if denom already includes sigma^2 return 0.0 if abs(denom) <= ANGLE_EPS else 0.25 / denom if weight_mode == 2: return 1.0 if weight_mode == 3: if sigma_or_weight is None: raise LatticeCalculationError("weight mode 3 requires an additional weight line") return sigma_or_weight if sigma_or_weight > 0.0 else 1.0 if weight_mode == 4: if sigma_or_weight is None: raise LatticeCalculationError("weight mode 4 requires an additional sigma line") # Weight = 1/sigma^2 return 0.0 if sigma_or_weight == 0.0 else 1.0 / (sigma_or_weight * sigma_or_weight) raise LatticeCalculationError(f"unsupported weight mode: {weight_mode}") def correction_terms(flags: list[int], transformed: float) -> list[float]: """ 指定された補正フラグと変換された位置データに基づいて補正項を計算する。 詳細説明: 回折データの測定には、様々な系統誤差（例:吸収、分極、ローレンツ因子など） が影響を与える可能性があります。この関数は、これらの誤差を補正するための 特定の補正項を計算し、最小二乗法のデザイン行列に追加できるようにします。 :param flags: list[int]: 5つの補正フラグのリスト。各フラグは特定の補正の有無を示す（0:なし、1:あり）。 :param transformed: float: 変換された回折位置データ（通常はsin(θ)）。 :returns: list[float]: 計算された補正項のリスト。 """ bas = transformed * transformed # sin(theta)^2 theta = math.asin(clamp_unit(transformed)) # theta in radians terms: list[float] = [] if flags[0]: # Correction 1: absorption? or 1/sin(2θ)? The original Fortran code suggests a zero shift correction # A.k.a. "dθ_c" in some contexts terms.append((math.pi / 2.0 - theta) * math.tan(math.pi / 2.0 - theta)) if flags[1]: # Correction 2: polarization factor, or L.P. factor component # sin^2(2θ) term, or 4 * sin^2(θ) * cos^2(θ) terms.append(4.0 * bas * (1.0 - bas)) if flags[2]: # Correction 3: usually related to sample transparency or beam divergence # The term 1/transformed + 1/theta is unusual. # If transformed = sin(theta), then 1/sin(theta) + 1/theta. # This term is only added if transformed and theta are not extremely small. if abs(transformed) <= ANGLE_EPS or abs(theta) <= ANGLE_EPS: terms.append(0.0) else: terms.append(4.0 * bas * (1.0 - bas) * (1.0 / transformed + 1.0 / theta)) if flags[3]: # Correction 4: Often related to preferred orientation or sample height error # sin(2θ) term terms.append(math.sin(2.0 * theta)) if flags[4]: # Correction 5: Another component, possibly related to preferred orientation or asymmetric peak broadening # sin(2θ) * cos(θ) term terms.append(math.sin(2.0 * theta) * math.cos(theta)) return terms def read_numeric_line(reader: LineReader, expected: int) -> list[float]: """ リーダーから一行を読み込み、指定された数の浮動小数点数に変換する。 詳細説明: 入力ファイルから設定値やパラメータを読み込む際に使用されるヘルパー関数。 読み込んだ行をスペースで分割し、最初の`expected`個の要素を浮動小数点数に変換します。 :param reader: LineReader: 入力ファイルを読み込むためのLineReaderインスタンス。 :param expected: int: 期待される数値の最小数。 :returns: list[float]: 読み込まれた数値のリスト。 :raises LatticeCalculationError: 期待されるよりも少ない数値が読み込まれた場合。 """ raw = reader.read_nonempty_line().split() if len(raw) < expected: raise LatticeCalculationError( f"line {reader.line_number}: expected at least {expected} numeric values, got {len(raw)}" ) return [float(x) for x in raw[:expected]] def parse_jobs(stream: TextIO) -> list[JobData]: """ 入力ストリームから複数の格子定数計算ジョブを解析する。 詳細説明: 入力ファイルは、一連の計算ジョブで構成されます。各ジョブはタイトル、 計算設定、および複数の回折観測データを含みます。この関数は、 ストリームの終端に達するか、特定の終了コード（`h=1000`で次のジョブ、`job_code=0`で全体終了） が読み込まれるまで、ジョブを順次解析して `JobData` オブジェクトのリストを生成します。 入力ファイルの主要なセクション: 1. ジョブタイトル 2. ジョブ設定パラメータ（格子系、補正フラグ、モードなど） 3. 波長と装置半径 4. 回折観測データのブロック - h k l raw_position - (position_mode=3の場合) zero_shift - (weight_mode in {0,1,3,4}の場合) sigma_or_weight - 特殊行: `2000 wavelength radius` で波長と半径を更新 - 特殊行: `1000 ...` で現在のジョブの終了を示す 5. ジョブ終了コード (0: 全体の終了, 1: 次のジョブへ) :param stream: TextIO: 入力データを含むテキストストリーム。 :returns: list[JobData]: 解析されたJobDataオブジェクトのリスト。 :raises EOFError: ファイルの予期せぬ終端に達した場合。 :raises LatticeCalculationError: 入力ファイルの形式が無効な場合（例: ヘッダー情報不足、反射行のデータ不足）。 """ reader = LineReader(stream) jobs: list[JobData] = [] while True: title = reader.read_line(allow_eof=True) if title is None: break # EOF encountered title = title.strip() if title == "": continue # Skip empty title lines # Read job configuration parameters params = reader.read_nonempty_line().split() if len(params) < 9: raise LatticeCalculationError(f"line {reader.line_number}: invalid job header") lattice_system, *raw_flags, weight_mode, position_mode, profile_mode = map(int, params[:9]) correction_flags = raw_flags[:5] # Take first 5 flags config = JobConfig( lattice_system=lattice_system, correction_flags=correction_flags, weight_mode=weight_mode, position_mode=position_mode, profile_mode=profile_mode, ) # Read initial wavelength and radius wavelength, radius = read_numeric_line(reader, 2) reference_wavelength = wavelength observations: list[Observation] = [] serial_index = 0 # Read reflection data line = reader.read_nonempty_line() while True: parts = line.split() if len(parts) < 4: raise LatticeCalculationError(f"line {reader.line_number}: invalid reflection line") h, k, l = map(int, parts[:3]) raw_position = float(parts[3]) if not (h < 1000 or h == 2000): # Check for job termination (h >= 1000 and h != 2000) break if h == 2000: # Special line: update wavelength and radius wavelength, radius = read_numeric_line(reader, 2) if len(observations) == 0: # If this is the first wavelength for the job reference_wavelength = wavelength line = reader.read_nonempty_line() continue serial_index += 1 zero_shift = None if position_mode == 3: zero_shift = float(reader.read_nonempty_line().strip()) transformed = transformed_position(raw_position, position_mode, wavelength, radius, zero_shift) sigma_or_weight = None if weight_mode in {0, 1, 3, 4}: sigma_or_weight = float(reader.read_nonempty_line().strip()) weight = calc_weight(weight_mode, transformed, sigma_or_weight) # Construct design row row = lattice_design_row(lattice_system, h, k, l) row = np.concatenate([row, np.array(correction_terms(correction_flags, transformed), dtype=float)]) # The y-value for least squares, typically proportional to 1/d^2 # Here it's 4.0 / (lambda^2) * sin^2(theta) = 1/d^2, if transformed is sin(theta) # The original code comment says "transformed_position (for least square)" and is (4.0/lambda^2)*sin^2(theta) # If transformed is already sin(theta), then this is 1/d^2. y_value = (4.0 / (wavelength * wavelength)) * transformed * transformed observations.append( Observation( serial_index=serial_index, h=h, k=k, l=l, raw_position=raw_position, transformed_position=y_value, weight=weight, design_row=row, wavelength=wavelength, ) ) line = reader.read_nonempty_line() jobs.append(JobData(title=title, config=config, observations=observations, reference_wavelength=reference_wavelength)) # Read job termination code job_code = int(reader.read_nonempty_line().strip()) if job_code == 0: break # All jobs finished return jobs # ---------- least squares ---------- def weighted_least_squares(observations: list[Observation], output: TextIO) -> LeastSquaresResult: """ 観測データに対して重み付き最小二乗法を実行する。 詳細説明: この関数は、与えられた観測データに対し、反復的に重み付き最小二乗法を適用します。 各反復では、現在のパラメータで残差を計算し、`OUTLIER_SIGMA` を超える残差を持つ観測データを 外れ値として除去します。外れ値がなくなるか、自由度が0になるまでこのプロセスを繰り返します。 最終的に、最適なパラメータ、その標準偏差、適合値、残差、および共分散行列が返されます。 :param observations: list[Observation]: 観測データのリスト。 :param output: TextIO: 外れ値の除去情報を書き込むための出力ストリーム。 :returns: LeastSquaresResult: 最小二乗計算の結果を含むオブジェクト。 :raises LatticeCalculationError: 観測データが提供されない場合、観測数よりもパラメータ数が多い場合、または正規方程式が特異な場合。 """ if not observations: raise LatticeCalculationError("no observations were provided") current = observations[:] n_params = current[0].design_row.size while True: y = np.array([obs.transformed_position for obs in current], dtype=float) w = np.array([obs.weight for obs in current], dtype=float) x = np.vstack([obs.design_row for obs in current]) if len(current) < n_params: raise LatticeCalculationError("number of observations is smaller than number of parameters") # Construct normal equations: (X^T W X) beta = X^T W y # X^T W X is 'normal' # X^T W y is 'rhs' normal = x.T @ (w[:, None] * x) rhs = x.T @ (w * y) try: covariance_unscaled = np.linalg.inv(normal) except np.linalg.LinAlgError as exc: raise LatticeCalculationError("normal matrix is singular") from exc params = covariance_unscaled @ rhs fitted = x @ params residuals = y - fitted dof = len(current) - n_params # Degrees of freedom if dof == 0: # If no degrees of freedom, cannot calculate scale factor for sigma sigma_obs = np.zeros(len(current), dtype=float) sigma_params = np.zeros(n_params, dtype=float) covariance = covariance_unscaled.copy() # No scaling return LeastSquaresResult(params, sigma_params, fitted, residuals, sigma_obs, current, covariance) weighted_ss = float(np.sum(w * residuals * residuals)) # Weighted sum of squares of residuals scale = math.sqrt(weighted_ss / dof) # Scale factor for standard deviation # Calculate observation sigma: sigma_obs = scale / sqrt(weight) sigma_obs = np.where(w > 0.0, scale / np.sqrt(w), 0.0) # Identify outliers outlier_mask = np.abs(residuals) > OUTLIER_SIGMA * sigma_obs if np.any(outlier_mask): # Outliers detected, remove them and re-run least squares removed = [obs for obs, flag in zip(current, outlier_mask) if flag] for obs, pred in zip(current, fitted): # The original code has an empty loop here, suggesting debug output # or a placeholder for detailed logging of all observations. pass for obs in removed: # Log removed observations pred_d = 1.0 / math.sqrt(max(obs.transformed_position, ANGLE_EPS)) # d_obs from original transformed position output.write( "This observation is eliminated.\n" f"No. = {obs.serial_index:2d} hkl=({obs.h:d} {obs.k:d} {obs.l:d}) raw={obs.raw_position:10.3f}\n" ) current = [obs for obs, flag in zip(current, outlier_mask) if not flag] continue # Repeat the loop with cleaned observations # No outliers, final result covariance = covariance_unscaled * (scale * scale) # Scale covariance matrix sigma_params = np.sqrt(np.maximum(np.diag(covariance), 0.0)) # Parameter standard deviations return LeastSquaresResult(params, sigma_params, fitted, residuals, sigma_obs, current, covariance) # ---------- post-processing ---------- def derive_cell_constants(lattice_system: int, params: np.ndarray, sigma: np.ndarray) -> CellConstants: """ 最小二乗法で得られたパラメータから逆格子および実格子の定数を導出する。 詳細説明: 最小二乗法によって決定された格子パラメータとそれらの標準偏差を受け取り、 指定された格子系の結晶学的な関係に基づいて、最終的な逆格子定数 (a*, b*, c*, alpha*, beta*, gamma*, V*) および実格子定数 (a, b, c, alpha, beta, gamma, V) を計算します。 各定数に対して誤差伝播法を用いて標準偏差も同時に計算されます。 :param lattice_system: int: 格子系の種類を示す整数コード（1-8）。 :param params: np.ndarray: 最小二乗法で得られたパラメータの配列。 :param sigma: np.ndarray: パラメータの標準偏差の配列。 :returns: CellConstants: 導出された格子定数を含むオブジェクト。 :raises LatticeCalculationError: 未対応の格子系が指定された場合。 """ cell = CellConstants() def reciprocal_to_direct_general(rec_lengths: np.ndarray, rec_length_sigma: np.ndarray, rec_cos: np.ndarray, rec_cos_sigma: np.ndarray) -> None: """ 一般的な逆格子定数から実格子定数を導出し、CellConstantsオブジェクトに格納するヘルパー関数。 """ cell.reciprocal_lengths[:] = rec_lengths cell.reciprocal_length_sigma[:] = rec_length_sigma cell.reciprocal_cosines[:] = rec_cos cell.reciprocal_cosine_sigma[:] = rec_cos_sigma for i in range(3): # Calculate reciprocal angles and their sigmas cell.reciprocal_angles_deg[i] = acos_deg(rec_cos[i]) cell.reciprocal_angle_sigma_deg[i] = SigmaPropagator([rec_cos[i]], [rec_cos_sigma[i]]).evaluate(lambda v: acos_deg(v[0]))[1] # Propagate error for direct volume prop = SigmaPropagator( [*rec_lengths, *rec_cos], [*rec_length_sigma, *rec_cos_sigma], ) cell.direct_volume, cell.direct_volume_sigma = prop.evaluate(direct_volume_from_reciprocal) # Calculate reciprocal volume (inverse of direct volume) if cell.direct_volume > 0.0: cell.reciprocal_volume = 1.0 / cell.direct_volume # Error propagation for 1/x: d(1/x) = -1/x^2 dx => sigma(1/x) = dx/x^2 cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume) for i in range(3): j = (i + 1) % 3 k = (i + 2) % 3 # Calculate direct cosines and their sigmas cell.direct_cosines[i], cell.direct_cosine_sigma[i] = SigmaPropagator( [rec_cos[i], rec_cos[j], rec_cos[k]], [rec_cos_sigma[i], rec_cos_sigma[j], rec_cos_sigma[k]], ).evaluate(direct_cosine_from_reciprocal) # Calculate direct angles and their sigmas cell.direct_angles_deg[i] = acos_deg(cell.direct_cosines[i]) cell.direct_angle_sigma_deg[i] = SigmaPropagator( [cell.direct_cosines[i]], [cell.direct_cosine_sigma[i]] ).evaluate(lambda v: acos_deg(v[0]))[1] # Calculate direct lengths and their sigmas cell.direct_lengths[i], cell.direct_length_sigma[i] = SigmaPropagator( [rec_lengths[i], rec_cos[i], rec_cos[j], rec_cos[k]], [rec_length_sigma[i], rec_cos_sigma[i], rec_cos_sigma[j], rec_cos_sigma[k]], ).evaluate(direct_length_from_reciprocal) if lattice_system == 1: # Triclinic # Params: a*^2, b*^2, c*^2, 2b*c*cos(alpha*), 2c*a*cos(beta*), 2a*b*cos(gamma*) rec_lengths = np.sqrt(np.maximum(params[:3], 0.0)) # Error propagation for sqrt(x): d(sqrt(x)) = 0.5 * dx / sqrt(x) rec_length_sigma = np.divide(sigma[:3], 2.0 * np.maximum(rec_lengths, ANGLE_EPS), out=np.zeros(3), where=rec_lengths > ANGLE_EPS) rec_cos = np.array([ params[3] / (rec_lengths[1] * rec_lengths[2]) if rec_lengths[1] * rec_lengths[2] > ANGLE_EPS else 0.0, params[4] / (rec_lengths[2] * rec_lengths[0]) if rec_lengths[2] * rec_lengths[0] > ANGLE_EPS else 0.0, params[5] / (rec_lengths[0] * rec_lengths[1]) if rec_lengths[0] * rec_lengths[1] > ANGLE_EPS else 0.0, ]) rec_cos_sigma = np.zeros(3, dtype=float) for i in range(3): j = (i + 1) % 3 k = (i + 2) % 3 # Error propagation for x/(y*z) rec_cos_sigma[i] = SigmaPropagator( [params[3 + i], rec_lengths[j], rec_lengths[k]], [sigma[3 + i], rec_length_sigma[j], rec_length_sigma[k]], ).evaluate(lambda v: 0.0 if abs(v[1] * v[2]) <= ANGLE_EPS else v[0] / (v[1] * v[2]))[1] reciprocal_to_direct_general(rec_lengths, rec_length_sigma, rec_cos, rec_cos_sigma) return cell if lattice_system == 2: # Monoclinic (unique b axis, beta*) # Parameters: a*^2, b*^2, c*^2, 2c*a*cos(beta*) # Missing 2b*c*cos(alpha*) (0), 2a*b*cos(gamma*) (0) expanded_p = np.array([params[0], params[1], params[2], 0.0, params[3], 0.0], dtype=float) expanded_s = np.array([sigma[0], sigma[1], sigma[2], 0.0, sigma[3], 0.0], dtype=float) return derive_cell_constants(1, expanded_p, expanded_s) if lattice_system == 3: # Monoclinic (unique c axis, gamma*) (as implemented) # Parameters: a*^2, b*^2, c*^2, 2a*b*cos(gamma*) # Missing 2b*c*cos(alpha*) (0), 2c*a*cos(beta*) (0) expanded_p = np.array([params[0], params[1], params[2], 0.0, 0.0, params[3]], dtype=float) expanded_s = np.array([sigma[0], sigma[1], sigma[2], 0.0, 0.0, sigma[3]], dtype=float) return derive_cell_constants(1, expanded_p, expanded_s) if lattice_system == 4: # Orthorhombic # Parameters: a*^2, b*^2, c*^2 (all angles are 90 deg, cosines are 0) rec_lengths = np.sqrt(np.maximum(params[:3], 0.0)) rec_length_sigma = np.divide(sigma[:3], 2.0 * np.maximum(rec_lengths, ANGLE_EPS), out=np.zeros(3), where=rec_lengths > ANGLE_EPS) reciprocal_to_direct_general(rec_lengths, rec_length_sigma, np.zeros(3), np.zeros(3)) # Reciprocal cosines and their sigmas are 0 return cell if lattice_system == 5: # Tetragonal # Parameters: a*^2, c*^2 (a*=b*, all angles 90 deg) a_star = safe_sqrt(params[0]) c_star = safe_sqrt(params[1]) da_star = 0.0 if a_star <= ANGLE_EPS else 0.5 * sigma[0] / a_star dc_star = 0.0 if c_star <= ANGLE_EPS else 0.5 * sigma[1] / c_star cell.reciprocal_lengths[:] = [a_star, a_star, c_star] cell.reciprocal_length_sigma[:] = [da_star, da_star, dc_star] # Direct lengths a = 1/a*, c = 1/c* a = 0.0 if a_star <= ANGLE_EPS else 1.0 / a_star c = 0.0 if c_star <= ANGLE_EPS else 1.0 / c_star cell.direct_lengths[:] = [a, a, c] cell.direct_length_sigma[:] = [da_star / max(a_star * a_star, ANGLE_EPS), da_star / max(a_star * a_star, ANGLE_EPS), dc_star / max(c_star * c_star, ANGLE_EPS)] cell.direct_volume = a * a * c # For tetragonal, direct_volume = a^2 * c cell.direct_volume_sigma = SigmaPropagator([a, a, c], [da_star / max(a_star * a_star, ANGLE_EPS), da_star / max(a_star * a_star, ANGLE_EPS), dc_star / max(c_star * c_star, ANGLE_EPS)]).evaluate(lambda v: v[0] * v[1] * v[2])[1] # Propagate error for a*a*c if cell.direct_volume > 0.0: cell.reciprocal_volume = 1.0 / cell.direct_volume cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume) return cell if lattice_system == 6: # Cubic # Parameters: a*^2 (a*=b*=c*, all angles 90 deg) a_star = safe_sqrt(params[0]) da_star = 0.0 if a_star <= ANGLE_EPS else 0.5 * sigma[0] / a_star cell.reciprocal_lengths[:] = [a_star, a_star, a_star] cell.reciprocal_length_sigma[:] = [da_star, da_star, da_star] # Direct length a = 1/a* a = 0.0 if a_star <= ANGLE_EPS else 1.0 / a_star da = da_star / max(a_star * a_star, ANGLE_EPS) cell.direct_lengths[:] = [a, a, a] cell.direct_length_sigma[:] = [da, da, da] cell.direct_volume = a**3 cell.direct_volume_sigma = 3.0 * a * a * da # For a^3, sigma is 3*a^2*da if cell.direct_volume > 0.0: cell.reciprocal_volume = 1.0 / cell.direct_volume cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume) return cell if lattice_system == 7: # Rhombohedral (trigonal) # Parameters: (a*^2 + b*^2 + c*^2), 2(b*c*cos(alpha*) + c*a*cos(beta*) + a*b*cos(gamma*)) # With a*=b*=c* and alpha*=beta*=gamma*, this simplifies: # P0 = 3*a*^2 # P1 = 6*a*^2*cos(alpha*) # So a*^2 = P0/3 => a* = sqrt(P0/3) # cos(alpha*) = P1/P0 # a_star = safe_sqrt(params[0] / 3.0) # params[0] is 3*a*^2 for rhombohedral a_star = safe_sqrt(params[0]) da_star = SigmaPropagator([params[0]], [sigma[0]]).evaluate(lambda v: safe_sqrt(v[0] / 3.0))[1] cos_alpha_star = 0.0 if abs(params[0]) <= ANGLE_EPS else params[1] / params[0] dcos_alpha_star = SigmaPropagator([params[0], params[1]], [sigma[0], sigma[1]]).evaluate( lambda v: 0.0 if abs(v[0]) <= ANGLE_EPS else v[1] / v[0] )[1] cell.reciprocal_lengths[:] = [a_star, a_star, a_star] cell.reciprocal_length_sigma[:] = [da_star, da_star, da_star] cell.reciprocal_cosines[:] = [cos_alpha_star, cos_alpha_star, cos_alpha_star] cell.reciprocal_cosine_sigma[:] = [dcos_alpha_star, dcos_alpha_star, dcos_alpha_star] angle_star = acos_deg(cos_alpha_star) dangle_star = SigmaPropagator([cos_alpha_star], [dcos_alpha_star]).evaluate(lambda v: acos_deg(v[0]))[1] cell.reciprocal_angles_deg[:] = [angle_star, angle_star, angle_star] cell.reciprocal_angle_sigma_deg[:] = [dangle_star, dangle_star, dangle_star] # Direct lengths, cosines, angles, and volume for trigonal cell.direct_lengths[0], cell.direct_length_sigma[0] = SigmaPropagator([a_star, cos_alpha_star], [da_star, dcos_alpha_star]).evaluate(trigonal_direct_length) cell.direct_lengths[:] = [cell.direct_lengths[0]] * 3 cell.direct_length_sigma[:] = [cell.direct_length_sigma[0]] * 3 cell.direct_cosines[0], cell.direct_cosine_sigma[0] = SigmaPropagator([cos_alpha_star], [dcos_alpha_star]).evaluate(trigonal_direct_cosine) cell.direct_cosines[:] = [cell.direct_cosines[0]] * 3 cell.direct_cosine_sigma[:] = [cell.direct_cosine_sigma[0]] * 3 angle = acos_deg(cell.direct_cosines[0]) dangle = SigmaPropagator([cell.direct_cosines[0]], [cell.direct_cosine_sigma[0]]).evaluate(lambda v: acos_deg(v[0]))[1] cell.direct_angles_deg[:] = [angle, angle, angle] cell.direct_angle_sigma_deg[:] = [dangle, dangle, dangle] cell.direct_volume, cell.direct_volume_sigma = SigmaPropagator([a_star, cos_alpha_star], [da_star, dcos_alpha_star]).evaluate(trigonal_direct_volume) if cell.direct_volume > 0.0: cell.reciprocal_volume = 1.0 / cell.direct_volume cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume) return cell if lattice_system == 8: # Hexagonal # Parameters: (a*^2 + b*^2 + a*b*), c*^2 # Which is 3*a*^2 for hexagonal a* = b* and gamma*=120deg. # So P0 = a*^2 + a*^2 + a*^2*(-0.5)*2 = 3*a*^2 # No, a*^2+k*^2+hk gives a_star^2 (from base_design_terms for h^2+k^2+2hk/2) # So params[0] is a*^2, params[1] is c*^2 a_star = safe_sqrt(params[0]) c_star = safe_sqrt(params[1]) da_star = 0.0 if a_star <= ANGLE_EPS else 0.5 * sigma[0] / a_star dc_star = 0.0 if c_star <= ANGLE_EPS else 0.5 * sigma[1] / c_star cell.reciprocal_lengths[:] = [a_star, a_star, c_star] cell.reciprocal_length_sigma[:] = [da_star, da_star, dc_star] # Direct lengths for hexagonal: a = 2/(sqrt(3)a*), c=1/c* a = 0.0 if a_star <= ANGLE_EPS else 2.0 / (math.sqrt(3.0) * a_star) c = 0.0 if c_star <= ANGLE_EPS else 1.0 / c_star # Error propagation for 1/x: d(1/x) = dx/x^2 da = 0.0 if a_star <= ANGLE_EPS else da_star * a / a_star # (2/sqrt(3)) * (-1/a*^2) * da* = (2/sqrt(3)a*) * (-1/a*) * da* = a * (-1/a*) * da* dc = 0.0 if c_star <= ANGLE_EPS else dc_star * c / c_star # c * (-1/c*) * dc* cell.direct_lengths[:] = [a, a, c] cell.direct_length_sigma[:] = [da, da, dc] # Angles for hexagonal cell.direct_cosines[2] = -0.5 # cos(120 deg) cell.reciprocal_cosines[2] = 0.5 # cos(60 deg) cell.direct_angles_deg[2] = 120.0 cell.reciprocal_angles_deg[2] = 60.0 # Direct volume for hexagonal V = sqrt(3)/2 * a^2 * c cell.direct_volume = math.sqrt(3.0) * 0.5 * a * a * c cell.direct_volume_sigma = SigmaPropagator([a, a, c], [da, da, dc]).evaluate(lambda v: math.sqrt(3.0) * 0.5 * v[0] * v[1] * v[2])[1] if cell.direct_volume > 0.0: cell.reciprocal_volume = 1.0 / cell.direct_volume cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume) return cell raise LatticeCalculationError(f"unsupported lattice system: {lattice_system}") # ---------- reporting ---------- def write_job_report(output: TextIO, job: JobData, fit: LeastSquaresResult, cell: CellConstants) -> None: """ 計算結果をフォーマットされたレポートとして出力ストリームに書き込む。 詳細説明: この関数は、単一の計算ジョブに対する詳細なレポートを生成します。 レポートには、ジョブのタイトル、設定、最小二乗法で得られたパラメータと誤差、 各観測データに対する計算値と残差、そして最終的に導出された逆格子および実格子の 定数とその標準偏差が含まれます。 :param output: TextIO: レポートを書き込むための出力ストリーム。 :param job: JobData: 計算に使用されたジョブデータ。 :param fit: LeastSquaresResult: 最小二乗法の結果。 :param cell: CellConstants: 導出された格子定数。 :returns: None """ output.write(f"\n {job.title}\n") cfg = job.config output.write(" HLS LE(1) LE(2) LE(3) LE(4) LE(5) IW IO IP\n") output.write( f"{cfg.lattice_system:7d}{cfg.correction_flags[0]:7d}{cfg.correction_flags[1]:7d}{cfg.correction_flags[2]:7d}" f"{cfg.correction_flags[3]:7d}{cfg.correction_flags[4]:7d}{cfg.weight_mode:7d}{cfg.position_mode:7d}{cfg.profile_mode:7d}\n" ) output.write(f"Wave Length ={job.reference_wavelength:f}\n") output.write(" P(I) SIGP(I)\n") for value, sigma in zip(fit.parameters, fit.parameter_sigma): output.write(f" {value:10.7f} {sigma:10.7f}\n") output.write(" H K L Dobs. Dcal. BO BC Resid. d(2Q)\n") for obs, fitted, resid in zip(fit.kept_observations, fit.fitted_y, fit.residuals): # Calculate d_obs and d_cal from transformed positions # transformed_position = 4 * sin^2(theta) / wavelength^2 = 1/d^2 # So d = 1 / sqrt(transformed_position) d_obs = 1.0 / math.sqrt(max(obs.transformed_position, ANGLE_EPS)) d_cal = 1.0 / math.sqrt(max(fitted, ANGLE_EPS)) # Calculate calculated 2theta from d_cal # 2d sin(theta) = wavelength => sin(theta) = wavelength / (2d) theta_calc = clamp_unit(obs.wavelength / (2.0 * d_cal)) if d_cal > 0 else 1.0 two_theta_calc = 2.0 * math.degrees(math.asin(theta_calc)) # Print observation details, fitted values, residuals, and 2theta difference output.write( f" ({obs.serial_index:2d}){obs.h:5d}{obs.k:5d}{obs.l:5d}{d_obs:8.4f}{d_cal:8.4f}" f"{obs.raw_position:8.3f}{two_theta_calc:8.3f}{resid:10.5f}{(obs.raw_position - two_theta_calc):10.5f}\n" ) output.write("Reciprocal cell constant\n") output.write(" a* b* c*\n") for j in range(3): output.write(f"{cell.reciprocal_lengths[j]:10.7f}({cell.reciprocal_length_sigma[j]:10.7f}) ") output.write("\n alpha* beta* gamma*\n") for j in range(3): output.write(f"{cell.reciprocal_angles_deg[j]:10.6f}({cell.reciprocal_angle_sigma_deg[j]:10.6f}) ") output.write("\n cos a* cos b* cos c*\n") for j in range(3): output.write(f"{cell.reciprocal_cosines[j]:10.6f}({cell.reciprocal_cosine_sigma[j]:10.6f}) ") output.write(f"\n V*={cell.reciprocal_volume:.10g}({cell.reciprocal_volume_sigma:.10g})\n") output.write("\nDirect cell constant\n") output.write(" a b c\n") for j in range(3): output.write(f"{cell.direct_lengths[j]:8.5f}({cell.direct_length_sigma[j]:8.5f}) ") output.write("\n alpha beta gamma\n") for j in range(3): output.write(f"{cell.direct_angles_deg[j]:8.4f}({cell.direct_angle_sigma_deg[j]:8.4f}) ") output.write("\n cos a cos b cos c\n") for j in range(3): output.write(f"{cell.direct_cosines[j]:8.4f}({cell.direct_cosine_sigma[j]:8.4f}) ") output.write(f"\n V={cell.direct_volume:.10g}({cell.direct_volume_sigma:.10g})\n") # ---------- CLI ---------- def build_argument_parser() -> argparse.ArgumentParser: """ コマンドライン引数を解析するためのパーサーを構築する。 詳細説明: この関数は、スクリプトが受け入れるコマンドライン引数を定義します。 具体的には、入力ファイルパス、出力ファイルパス、およびバナー出力を抑制するための `--silent` フラグを設定します。 :returns: argparse.ArgumentParser: 構築されたArgumentParserオブジェクト。 """ parser = argparse.ArgumentParser(description="Least-square calculation of lattice constants from reflection data") parser.add_argument("input_file", type=Path, help="input file") parser.add_argument("output_file", type=Path, help="output report file") parser.add_argument("--silent", action="store_true", help="suppress banner output") return parser def main() -> int: """ スクリプトのエントリポイント。 詳細説明: コマンドライン引数を解析し、入力ファイルから格子定数計算ジョブを読み込みます。 各ジョブに対して重み付き最小二乗計算を実行し、得られた結果から格子定数を導出します。 最終的なレポートは指定された出力ファイルに書き込まれます。 :returns: int: 終了コード（成功の場合は0）。 """ parser = build_argument_parser() args = parser.parse_args() if not args.silent: print("Least Square Calculation for Lattice Constant") print("Refactored Python version") print(f"Input : {args.input_file}") print(f"Output: {args.output_file}") with args.input_file.open("r", encoding="utf-8") as fin: jobs = parse_jobs(fin) with args.output_file.open("w", encoding="utf-8") as fout: fout.write("Least Square Calculation for Lattice Constant\n") fout.write("Refactored Python version\n\n") for job in jobs: fit = weighted_least_squares(job.observations, fout) cell = derive_cell_constants(job.config.lattice_system, fit.parameters, fit.parameter_sigma) write_job_report(fout, job, fit, cell) fout.write("\nLattice Constant Calculation was finished!\n") if not args.silent: print("Lattice Constant Calculation was finished!") return 0 if __name__ == "__main__": raise SystemExit(main())