1+ import numpy as np
2+ import pandas as pd
3+ import matplotlib .pyplot as plt
4+ from scipy .integrate import quad
5+ from scipy .interpolate import interp1d
6+
7+ # =========================
8+ # 前回の最適化結果(勝利パラメータ)
9+ # =========================
10+ res_csgt_x = [0.010 , 0.100 , 1.200 , - 19.34 , 72.99 , 0.35 , - 1.10 ] # A, sigma, zp, M, H0, Om, w_off
11+ res_lcdm_x = [- 19.34 , 72.99 , 0.35 ] # M, H0, Om
12+
13+ # =========================
14+ # グラフ用計算関数
15+ # =========================
16+ def get_plot_data (params , is_csgt = True ):
17+ if is_csgt :
18+ A , sigma , zp , M , H0 , Om , w_off = params
19+ else :
20+ M , H0 , Om = params
21+ A , sigma , zp , w_off = 0.0 , 1.0 , 0.7 , - 1.0
22+
23+ Ode = 1.0 - Om
24+ z_range = np .linspace (0.001 , 2.3 , 200 )
25+
26+ # w(z) の計算
27+ w_z = w_off + A * np .exp (- (z_range - zp )** 2 / (2 * sigma ** 2 ))
28+
29+ # H(z) の計算
30+ H_z = []
31+ for z in z_range :
32+ w_int , _ = quad (lambda zp_val : ((1.0 + w_off ) + A * np .exp (- (zp_val - zp )** 2 / (2 * sigma ** 2 ))) / (1.0 + zp_val ), 0 , z )
33+ Ez = np .sqrt (Om * (1 + z )** 3 + Ode * np .exp (3.0 * w_int ))
34+ H_z .append (H0 * Ez )
35+
36+ return z_range , w_z , np .array (H_z )
37+
38+ # データの取得
39+ z_axis , w_csgt , H_csgt = get_plot_data (res_csgt_x , True )
40+ _ , w_lcdm , H_lcdm = get_plot_data (res_lcdm_x , False )
41+
42+ # =========================
43+ # プロット作成
44+ # =========================
45+ fig , ax = plt .subplots (1 , 2 , figsize = (15 , 6 ))
46+
47+ # 左図:状態方程式 w(z)
48+ ax [0 ].plot (z_axis , w_csgt , 'r-' , lw = 2 , label = 'CSGT (A=0.01, z_p=1.2)' )
49+ ax [0 ].plot (z_axis , w_lcdm , 'b--' , lw = 2 , label = 'ΛCDM (w=-1)' )
50+ ax [0 ].axhline (- 1 , color = 'gray' , linestyle = ':' , alpha = 0.5 )
51+ ax [0 ].set_xlabel ('Redshift z' , fontsize = 12 )
52+ ax [0 ].set_ylabel ('Equation of State w(z)' , fontsize = 12 )
53+ ax [0 ].set_title ('Dark Energy Evolution' , fontsize = 14 )
54+ ax [0 ].legend ()
55+ ax [0 ].grid (alpha = 0.3 )
56+
57+ # 右図:膨張率 H(z) の比較(偏差)
58+ # 標準モデルからのズレをパーセントで表示すると分かりやすいわ
59+ H_diff = (H_csgt - H_lcdm ) / H_lcdm * 100
60+ ax [1 ].plot (z_axis , H_diff , 'g-' , lw = 2 , label = '(H_CSGT - H_LCDM) / H_LCDM [%]' )
61+ ax [1 ].axhline (0 , color = 'blue' , linestyle = '--' , alpha = 0.5 )
62+ ax [1 ].set_xlabel ('Redshift z' , fontsize = 12 )
63+ ax [1 ].set_ylabel ('Difference in H(z) (%)' , fontsize = 12 )
64+ ax [1 ].set_title ('Expansion Rate Deviation from ΛCDM' , fontsize = 14 )
65+ ax [1 ].legend ()
66+ ax [1 ].grid (alpha = 0.3 )
67+
68+ plt .tight_layout ()
69+ plt .show ()
70+
71+ print ("Graphs generated successfully." )
0 commit comments