""" 最小二乗法を用いた格子定数計算スクリプト。 このスクリプトは、X線回折データなどの回折データから結晶の格子定数を導出します。 様々な格子系に対応し、観測データの前処理、重み付け、最小二乗フィッティング、 そして最終的な格子定数（逆格子および実格子）の導出を行います。 外れ値の自動除去機能も備えています。 詳細説明: 入力ファイルは特定のフォーマットに従い、複数の計算ジョブを連続して処理できます。 各ジョブでは、格子系、補正フラグ、重み付けモード、回折位置の変換モードなどが 設定され、それに従って観測データが処理されます。 重み付き最小二乗法によりパラメータが決定され、その結果から逆格子定数と実格子定数、 およびそれらの標準偏差が計算されます。計算結果は整形されたレポートファイルとして出力されます。 関連リンク: :doc:`lsq_latt2_usage` """ from __future__ import annotations import argparse import math import sys from dataclasses import dataclass, field from pathlib import Path from typing import Callable, Iterable, TextIO import matplotlib.pyplot as plt 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: """ 最小二乗法の結果を保持するデータクラス。 """ 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 rejected_observations: list[Observation] = field(default_factory=list) @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: """ 観測データに対して重み付き最小二乗法を実行する。 外れ値は逐次除去するが、詳細な報告は最終レポート側でまとめて行う。 """ if not observations: raise LatticeCalculationError("no observations were provided") current = observations[:] rejected: list[Observation] = [] 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") 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 if dof == 0: sigma_obs = np.zeros(len(current), dtype=float) sigma_params = np.zeros(n_params, dtype=float) covariance = covariance_unscaled.copy() return LeastSquaresResult( params, sigma_params, fitted, residuals, sigma_obs, current, covariance, rejected, ) weighted_ss = float(np.sum(w * residuals * residuals)) scale = math.sqrt(weighted_ss / dof) sigma_obs = np.where(w > 0.0, scale / np.sqrt(w), 0.0) outlier_mask = np.abs(residuals) > OUTLIER_SIGMA * sigma_obs if np.any(outlier_mask): removed = [obs for obs, flag in zip(current, outlier_mask) if flag] rejected.extend(removed) current = [obs for obs, flag in zip(current, outlier_mask) if not flag] continue covariance = covariance_unscaled * (scale * scale) sigma_params = np.sqrt(np.maximum(np.diag(covariance), 0.0)) return LeastSquaresResult( params, sigma_params, fitted, residuals, sigma_obs, current, covariance, rejected, ) # ---------- 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 ---------- class TeeStream: """ファイルとコンソールへ同時に書き込む簡易ストリーム。""" def __init__(self, *streams: TextIO): self.streams = streams def write(self, text: str) -> int: for stream in self.streams: stream.write(text) stream.flush() return len(text) def flush(self) -> None: for stream in self.streams: stream.flush() def make_default_output_path(input_path: Path) -> Path: """出力ファイル未指定時は入力ファイルの拡張子を .out に置き換える。""" return input_path.with_suffix('.out') def lattice_system_name(code: int) -> str: names = { 1: "Triclinic / 三斜晶", 2: "Monoclinic(unique b) / 単斜晶(b軸特異)", 3: "Monoclinic(unique c) / 単斜晶(c軸特異)", 4: "Orthorhombic / 斜方晶", 5: "Tetragonal / 正方晶", 6: "Cubic / 立方晶", 7: "Rhombohedral(Trigonal) / 菱面体晶(三方晶)", 8: "Hexagonal / 六方晶", } return names.get(code, f"Unknown({code})") def lattice_parameter_labels(lattice_system: int) -> list[str]: if lattice_system == 1: return [ "a*^2", "b*^2", "c*^2", "2 b* c* cos(alpha*)", "2 c* a* cos(beta*)", "2 a* b* cos(gamma*)", ] if lattice_system == 2: return [ "a*^2", "b*^2", "c*^2", "2 c* a* cos(beta*)", ] if lattice_system == 3: return [ "a*^2", "b*^2", "c*^2", "2 a* b* cos(gamma*)", ] if lattice_system == 4: return ["a*^2", "b*^2", "c*^2"] if lattice_system == 5: return ["a*^2 (= b*^2)", "c*^2"] if lattice_system == 6: return ["a*^2 (= b*^2 = c*^2)"] if lattice_system == 7: return [ "a*^2-like term (see rhombohedral model)", "cos(alpha*)-related cross term", ] if lattice_system == 8: return ["a*^2-like term", "c*^2"] return [f"lattice_param_{i+1}" for i in range(6)] def correction_metadata() -> list[dict[str, str]]: return [ { "name": "C1", "formula": "(pi/2 - theta) * tan(pi/2 - theta)", "en": "Possible interpretation: specimen displacement / zero-shift related term.", "ja": "推定される意味: 試料位置ずれ / ゼロシフト関連の補正項。", }, { "name": "C2", "formula": "4*s^2*(1 - s^2)", "en": "Possible interpretation: Lorentz-polarization related term.", "ja": "推定される意味: ローレンツ・偏光因子関連の補正項。", }, { "name": "C3", "formula": "4*s^2*(1 - s^2)*(1/s + 1/theta)", "en": "Possible interpretation: transparency / beam-divergence related term.", "ja": "推定される意味: 試料透過 / ビーム発散関連の補正項。", }, { "name": "C4", "formula": "sin(2*theta)", "en": "Possible interpretation: sample-height / specimen-displacement / preferred-orientation related term.", "ja": "推定される意味: 試料高さずれ / 試料位置ずれ / 配向性関連の補正項。", }, { "name": "C5", "formula": "sin(2*theta) * cos(theta)", "en": "Possible interpretation: asymmetry / orientation related term.", "ja": "推定される意味: 非対称性 / 配向性関連の補正項。", }, ] def format_value_sigma(value: float, sigma: float, *, width: int = 14, prec: int = 7) -> str: return f"{value:{width}.{prec}f} ± {sigma:{width}.{prec}f}" def write_section_title(output: TextIO, title_en: str, title_ja: str) -> None: output.write(f"\n[{title_en} / {title_ja}]\n") def parameter_blocks(job: JobData, fit: LeastSquaresResult) -> tuple[list[tuple[str, float, float]], list[tuple[str, float, float]]]: lattice_labels = lattice_parameter_labels(job.config.lattice_system) n_lattice = len(lattice_labels) lattice_part = [ (label, float(value), float(sigma)) for label, value, sigma in zip( lattice_labels, fit.parameters[:n_lattice], fit.parameter_sigma[:n_lattice], ) ] corr_meta = correction_metadata() corr_part: list[tuple[str, float, float]] = [] corr_index = n_lattice for flag, meta in zip(job.config.correction_flags, corr_meta): if flag: corr_part.append( ( meta["name"], float(fit.parameters[corr_index]), float(fit.parameter_sigma[corr_index]), ) ) corr_index += 1 return lattice_part, corr_part def calc_fit_statistics(fit: LeastSquaresResult) -> tuple[float, float, int]: if len(fit.kept_observations) == 0: return 0.0, 0.0, 0 w = np.array([obs.weight for obs in fit.kept_observations], dtype=float) y = np.array([obs.transformed_position for obs in fit.kept_observations], dtype=float) resid = fit.residuals n_params = len(fit.parameters) dof = max(len(y) - n_params, 0) weighted_ss = float(np.sum(w * resid * resid)) rms = math.sqrt(float(np.mean(resid * resid))) if len(resid) else 0.0 return weighted_ss, rms, dof def mode_description_position(mode: int) -> str: mapping = { 1: "sin(theta) from angle in degree / 角度(deg)から sin(theta) へ変換", 2: "raw position / 生の位置値をそのまま使用", 3: "zero-shift corrected sin(theta) / ゼロシフト補正付き sin(theta)", 4: "sin(theta/2) / sin(theta/2)", 5: "lambda/(2*d)-like transform / lambda/(2*d) 型変換", } return mapping.get(mode, "geometry-dependent transform / 幾何条件依存の変換") def mode_description_weight(mode: int) -> str: mapping = { 0: "sigma(theta)-based weight / sigma(theta) ベース重み", 1: "sigma(2theta)-based weight / sigma(2theta) ベース重み", 2: "constant weight / 一定重み", 3: "user-supplied weight / ユーザー指定重み", 4: "1/sigma^2 / 1/sigma^2", } return mapping.get(mode, "unknown / 不明") def write_job_report(output: TextIO, job: JobData, fit: LeastSquaresResult, cell: CellConstants) -> None: """ 計算結果を見やすく整理したレポートとして出力する。 """ cfg = job.config lattice_part, corr_part = parameter_blocks(job, fit) weighted_ss, rms_resid, dof = calc_fit_statistics(fit) output.write("\n" + "=" * 88 + "\n") output.write(f"Job: {job.title.strip() or '(untitled job)'}\n") output.write("=" * 88 + "\n") write_section_title(output, "Job summary", "ジョブ概要") output.write(f"Lattice system / 格子系 : {lattice_system_name(cfg.lattice_system)}\n") output.write(f"Reference wavelength / 参照波長 : {job.reference_wavelength:.7f}\n") output.write(f"Weight mode / 重みモード : {cfg.weight_mode} ({mode_description_weight(cfg.weight_mode)})\n") output.write(f"Position mode / 位置変換モード : {cfg.position_mode} ({mode_description_position(cfg.position_mode)})\n") output.write(f"Profile mode / プロファイルモード : {cfg.profile_mode}\n") output.write(f"Input observations / 入力反射数 : {len(job.observations)}\n") output.write(f"Used observations / 採用反射数 : {len(fit.kept_observations)}\n") output.write(f"Rejected observations / 除外反射数 : {len(fit.rejected_observations)}\n") output.write(f"Degrees of freedom / 自由度 : {dof}\n") output.write(f"Weighted SS / 重み付き残差二乗和 : {weighted_ss:.7e}\n") output.write(f"RMS residual in y / y空間RMS残差 : {rms_resid:.7e}\n") write_section_title(output, "Direct lattice constants", "実格子定数") output.write(f"a = {format_value_sigma(cell.direct_lengths[0], cell.direct_length_sigma[0])}\n") output.write(f"b = {format_value_sigma(cell.direct_lengths[1], cell.direct_length_sigma[1])}\n") output.write(f"c = {format_value_sigma(cell.direct_lengths[2], cell.direct_length_sigma[2])}\n") output.write(f"alpha = {format_value_sigma(cell.direct_angles_deg[0], cell.direct_angle_sigma_deg[0])} deg\n") output.write(f"beta = {format_value_sigma(cell.direct_angles_deg[1], cell.direct_angle_sigma_deg[1])} deg\n") output.write(f"gamma = {format_value_sigma(cell.direct_angles_deg[2], cell.direct_angle_sigma_deg[2])} deg\n") output.write(f"cos(a) = {format_value_sigma(cell.direct_cosines[0], cell.direct_cosine_sigma[0])}\n") output.write(f"cos(b) = {format_value_sigma(cell.direct_cosines[1], cell.direct_cosine_sigma[1])}\n") output.write(f"cos(c) = {format_value_sigma(cell.direct_cosines[2], cell.direct_cosine_sigma[2])}\n") output.write(f"V = {format_value_sigma(cell.direct_volume, cell.direct_volume_sigma)}\n") write_section_title(output, "Reciprocal lattice constants", "逆格子定数") output.write(f"a* = {format_value_sigma(cell.reciprocal_lengths[0], cell.reciprocal_length_sigma[0])}\n") output.write(f"b* = {format_value_sigma(cell.reciprocal_lengths[1], cell.reciprocal_length_sigma[1])}\n") output.write(f"c* = {format_value_sigma(cell.reciprocal_lengths[2], cell.reciprocal_length_sigma[2])}\n") output.write(f"alpha* = {format_value_sigma(cell.reciprocal_angles_deg[0], cell.reciprocal_angle_sigma_deg[0])} deg\n") output.write(f"beta* = {format_value_sigma(cell.reciprocal_angles_deg[1], cell.reciprocal_angle_sigma_deg[1])} deg\n") output.write(f"gamma* = {format_value_sigma(cell.reciprocal_angles_deg[2], cell.reciprocal_angle_sigma_deg[2])} deg\n") output.write(f"cos(a*) = {format_value_sigma(cell.reciprocal_cosines[0], cell.reciprocal_cosine_sigma[0])}\n") output.write(f"cos(b*) = {format_value_sigma(cell.reciprocal_cosines[1], cell.reciprocal_cosine_sigma[1])}\n") output.write(f"cos(c*) = {format_value_sigma(cell.reciprocal_cosines[2], cell.reciprocal_cosine_sigma[2])}\n") output.write(f"V* = {format_value_sigma(cell.reciprocal_volume, cell.reciprocal_volume_sigma)}\n") write_section_title(output, "Lattice fit parameters", "格子フィットパラメータ") for name, value, sigma in lattice_part: output.write(f"{name:<36s} = {format_value_sigma(value, sigma)}\n") write_section_title(output, "Correction model", "補正モデル") output.write("s = transformed_position, theta = diffraction angle in radian / θはラジアンでの回折角\n") output.write("The following physical meanings are tentative interpretations inferred from the legacy formulas and comments.\n") output.write("以下の物理的意味は旧コードの式とコメントから推定した暫定解釈です。\n") output.write("The user is responsible for confirming each interpretation for the actual experimental setup.\n") output.write("各補正項の解釈が実験条件に対して妥当かどうかの最終確認は、ユーザーの責任で行ってください。\n") active_corr_names = {name for name, _, _ in corr_part} fitted_corr_map = {name: (value, sigma) for name, value, sigma in corr_part} for flag, meta in zip(cfg.correction_flags, correction_metadata()): state = "ON" if flag else "OFF" output.write(f"\n{meta['name']} : {state}\n") output.write(f" formula / 式 : {meta['formula']}\n") output.write(f" {meta['en']}\n") output.write(f" {meta['ja']}\n") output.write(" This interpretation is inferred from the legacy formula and comments; the user is responsible for confirming it.\n") output.write(" この解釈は旧コードの式とコメントからの推定であり、最終確認はユーザーの責任で行ってください。\n") if meta["name"] in active_corr_names: value, sigma = fitted_corr_map[meta["name"]] output.write(f" fitted coefficient / フィット係数 : {format_value_sigma(value, sigma)}\n") if corr_part: write_section_title(output, "Fitted correction coefficients", "フィットされた補正係数") for name, value, sigma in corr_part: output.write(f"{name:<8s} = {format_value_sigma(value, sigma)}\n") else: write_section_title(output, "Fitted correction coefficients", "フィットされた補正係数") output.write("No correction term is enabled. / 補正項は有効化されていません。\n") write_section_title(output, "Reflection table", "反射一覧") output.write( f"{'No.':>4s} {'h':>4s} {'k':>4s} {'l':>4s} " f"{'d_obs':>12s} {'d_calc':>12s} {'raw_pos':>12s} {'2theta_calc':>12s} " f"{'residual_y':>14s} {'delta_2theta':>14s} {'weight':>12s}\n" ) for obs, fitted, resid, sigma_obs in zip( fit.kept_observations, fit.fitted_y, fit.residuals, fit.observation_sigma, ): d_obs = 1.0 / math.sqrt(max(obs.transformed_position, ANGLE_EPS)) d_calc = 1.0 / math.sqrt(max(fitted, ANGLE_EPS)) theta_calc = clamp_unit(obs.wavelength / (2.0 * d_calc)) if d_calc > 0.0 else 1.0 two_theta_calc = 2.0 * math.degrees(math.asin(theta_calc)) delta_2theta = obs.raw_position - two_theta_calc output.write( f"{obs.serial_index:4d} {obs.h:4d} {obs.k:4d} {obs.l:4d} " f"{d_obs:12.6f} {d_calc:12.6f} {obs.raw_position:12.5f} {two_theta_calc:12.5f} " f"{resid:14.6e} {delta_2theta:14.6e} {obs.weight:12.5e}\n" ) write_section_title(output, "Rejected reflections", "除外された反射") if fit.rejected_observations: output.write(f"Outlier criterion / 外れ値判定基準 : |residual| > {OUTLIER_SIGMA:.1f} * sigma_obs\n") output.write(f"{'No.':>4s} {'h':>4s} {'k':>4s} {'l':>4s} {'raw_pos':>12s} {'note':>18s}\n") for obs in fit.rejected_observations: output.write( f"{obs.serial_index:4d} {obs.h:4d} {obs.k:4d} {obs.l:4d} " f"{obs.raw_position:12.5f} {'rejected as outlier':>18s}\n" ) else: output.write("No reflection was rejected. / 除外された反射はありません。\n") def sanitize_stem(text: str) -> str: return ''.join(ch if ch.isalnum() or ch in ('-', '_') else '_' for ch in text).strip('_') or 'job' def plot_job_result( job: JobData, fit: LeastSquaresResult, output_file: Path, job_index: int, *, show_plot: bool = True, ) -> Path: """各ジョブの Q-1/d プロットを保存する。""" one_over_d_calc = np.sqrt(np.maximum(fit.fitted_y, 0.0)) one_over_d_obs = np.sqrt(np.maximum([obs.transformed_position for obs in fit.kept_observations], 0.0)) sigma_one_over_d = np.zeros_like(one_over_d_calc) valid = one_over_d_calc > ANGLE_EPS sigma_one_over_d[valid] = fit.observation_sigma[valid] / (2.0 * one_over_d_calc[valid]) q_obs = 2.0 * math.pi * one_over_d_obs q_min = 0.0 q_max = max(float(np.max(q_obs)) if q_obs.size else 1.0, 1.0) * 1.05 q_line = np.linspace(q_min, q_max, 400) one_over_d_line = q_line / (2.0 * math.pi) fig, ax_left = plt.subplots(figsize=(8, 5)) ax_left.plot(q_line, one_over_d_line, label='theory: 1/d = Q / 2π') ax_left.errorbar(q_obs, one_over_d_calc, yerr=sigma_one_over_d, fmt='o', capsize=3, label='calculated 1/d ±σ') ax_left.plot(q_obs, one_over_d_obs, 's', label='observed 1/d') ax_left.set_xlabel('Scattering vector Q') ax_left.set_ylabel('1/d') ax_left.grid(True, alpha=0.3) ax_right = ax_left.twinx() ax_right.plot(q_obs, sigma_one_over_d, '^--', label='σ(1/d)') ax_right.set_ylabel('σ(1/d)') title = job.title.strip() or f'Job {job_index}' ax_left.set_title(title) handles_left, labels_left = ax_left.get_legend_handles_labels() handles_right, labels_right = ax_right.get_legend_handles_labels() ax_left.legend(handles_left + handles_right, labels_left + labels_right, loc='best') fig.tight_layout() plot_path = output_file.with_name(f"{output_file.stem}_job{job_index:02d}_{sanitize_stem(title)[:40]}.png") fig.savefig(plot_path, dpi=150) if show_plot: plt.show() plt.close(fig) return plot_path # ---------- CLI ---------- def build_argument_parser() -> argparse.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, nargs='?', help='output report file (.out by default)') parser.add_argument('--silent', action='store_true', help='suppress banner output') parser.add_argument('--no-show', action='store_true', help='save plots but do not open plot windows') return parser def main() -> int: """スクリプトのエントリポイント。""" parser = build_argument_parser() args = parser.parse_args() output_file = args.output_file if args.output_file is not None else make_default_output_path(args.input_file) if not args.silent: print('Least Square Calculation for Lattice Constant') print('Refactored Python version') print(f'Input : {args.input_file}') print(f'Output: {output_file}') with args.input_file.open('r', encoding='utf-8') as fin: jobs = parse_jobs(fin) plot_paths: list[Path] = [] with output_file.open('w', encoding='utf-8') as fout: tee = TeeStream(fout, sys.stdout) tee.write('Least Square Calculation for Lattice Constant\n') tee.write('Refactored Python version\n\n') for idx, job in enumerate(jobs, start=1): fit = weighted_least_squares(job.observations, tee) cell = derive_cell_constants(job.config.lattice_system, fit.parameters, fit.parameter_sigma) write_job_report(tee, job, fit, cell) plot_paths.append(plot_job_result(job, fit, output_file, idx, show_plot=not args.no_show)) tee.write('\nLattice Constant Calculation was finished!\n') if not args.silent: print('Saved plot files:') for path in plot_paths: print(f' {path}') return 0 if __name__ == '__main__': raise SystemExit(main())