"""
このモジュールは、宇宙論的シミュレーションデータの処理に役立つ関数群を提供します。
主に以下の機能を含みます：
- 密度のラジアルプロファイル作成
- 周期境界条件下での座標変換
- comoving座標から物理座標(kpc)への変換
定数:
--------
h : float
ハッブル定数の無次元パラメータ（H0 = 100 h km/s/Mpc単位）。
rho_crit : float
宇宙の臨界密度 (単位: M_sun / Mpc^3)。
redshift : np.ndarray
外部ファイルから読み込んだ赤方偏移データ。
scale_factor : np.ndarray
外部ファイルから読み込んだスケール因子データ。
"""
import numpy as np
h = 0.67742
"""ハッブル定数の無次元化パラメータ。"""
rho_crit = 1.27e11
"""宇宙の臨界密度 [M_sun / Mpc^3]。"""
redshift = np.loadtxt("/home/nishihama/erosita/241119/redshift.txt", unpack=True)[2]
"""赤方偏移データを外部ファイルから読み込み。"""
scale_factor = np.loadtxt("/home/nishihama/erosita/241119/redshift.txt", unpack=True)[3]
"""スケール因子データを外部ファイルから読み込み。"""
def create_radial_profile_bin(density, x, xmin=-3, xmax=1, n_bins=30):
"""
ログスケールで分割した半径ビンごとに密度の統計量（中央値、25パーセンタイル、75パーセンタイル）を計算する。
Parameters
----------
density : np.ndarray
各点の密度データ。
x : np.ndarray
半径データ（logスケールでビン分けする対象）。
xmin : float, optional
ビンの最小値（log10スケール）, デフォルトは-3。
xmax : float, optional
ビンの最大値（log10スケール）, デフォルトは1。
n_bins : int, optional
ビンの数, デフォルトは30。
Returns
-------
x_common : np.ndarray
各ビンの中心値。
fidu_d_median : np.ndarray
各ビン内の密度の中央値。
fidu_25tile : np.ndarray
各ビン内の密度の25パーセンタイル値。
fidu_75tile : np.ndarray
各ビン内の密度の75パーセンタイル値。
"""
x_common_egde = np.logspace(xmin, xmax, n_bins+1)
x_common = np.logspace(xmin, xmax, n_bins)
fidu_indices = np.digitize(x, x_common_egde, right=True)
fidu_d_grouped = [density[fidu_indices == j] for j in range(n_bins+1)]
fidu_d_grouped = fidu_d_grouped[1:]
fidu_d_median = [np.nanmedian(i) for i in fidu_d_grouped]
fidu_25tile = [np.nanpercentile(i, 25) for i in fidu_d_grouped]
fidu_75tile = [np.nanpercentile(i, 75) for i in fidu_d_grouped]
return x_common, np.array(fidu_d_median), np.array(fidu_25tile), np.array(fidu_75tile)
def conv(x, y, z, cx, cy, cz, L):
"""
周期境界条件を考慮した座標変換を行う。
Parameters
----------
x, y, z : np.ndarray
各軸の座標配列。
cx, cy, cz : float
変換の基準となる中心座標。
L : float
ボックスサイズ（周期境界条件のスケール）。
Returns
-------
x_new, y_new, z_new : np.ndarray
変換後の座標。
"""
x = x - cx - L * np.round((x - cx) / L)
y = y - cy - L * np.round((y - cy) / L)
z = z - cz - L * np.round((z - cz) / L)
return x, y, z
def conv2(coordinates, c, L):
"""
周期境界条件を考慮して3次元座標配列を変換する。
Parameters
----------
coordinates : np.ndarray
形状 (N, 3) の座標配列。
c : tuple of float
中心座標 (cx, cy, cz)。
L : float
ボックスサイズ。
Returns
-------
np.ndarray
変換後の座標配列 (N, 3)。
"""
cx, cy, cz = c
x, y, z = coordinates[:, 0], coordinates[:, 1], coordinates[:, 2]
x = x - cx - L * np.round((x - cx) / L)
y = y - cy - L * np.round((y - cy) / L)
z = z - cz - L * np.round((z - cz) / L)
return np.column_stack((x, y, z))
def to_kpc(x, a):
"""
comoving座標を物理スケール（kpc単位）に変換する。
Parameters
----------
x : float or np.ndarray
comoving座標の値。
a : float
スケール因子。
Returns
-------
float or np.ndarray
物理スケールでの座標値（kpc単位）。
"""
return x * 1000 / h * a