PHASE 2: FINAL — CF4 DENSITY GRID # Using official CosmicFlows-4 density field
# ============================================================
# PHASE 2: FINAL — CF4 DENSITY GRID
# Using official CosmicFlows-4 density field
# ============================================================
import numpy as np
import pandas as pd
from astropy.io import fits
from scipy.interpolate import RegularGridInterpolator
from scipy.stats import ks_2samp, ttest_ind
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
import astropy.units as u
print("=" * 60)
print("PHASE 2: FINAL — CF4 DENSITY GRID")
print("=" * 60)
# ------------------------------------------------------------
# STEP 1: Load CF4 density grid
# ------------------------------------------------------------
print("\n1. Loading CF4 density grid...")
hdul = fits.open('CF4_new_64-z008_delta.fits')
density_grid = hdul[0].data
header = hdul[0].header
print(f" Grid shape: {density_grid.shape}")
print(f" Density range: {density_grid.min():.4f} to {density_grid.max():.4f}")
print(f" Mean: {density_grid.mean():.4f}")
# Get grid coordinates from header
# CF4 uses Supergalactic Cartesian coordinates (SGX, SGY, SGZ)
# Typical range: -500 to +500 Mpc/h for 64³ grid
n = density_grid.shape[0]
coords_sg = np.linspace(-500, 500, n) # Mpc/h
print(f" Grid coordinates: {coords_sg[0]:.0f} to {coords_sg[-1]:.0f} Mpc/h")
# Create interpolator
interp_density = RegularGridInterpolator(
(coords_sg, coords_sg, coords_sg),
density_grid,
bounds_error=False,
fill_value=1.0
)
# ------------------------------------------------------------
# STEP 2: Load your 20 cleaned galaxies
# ------------------------------------------------------------
print("\n2. Loading cleaned galaxies...")
df = pd.read_csv('sparc_phase2_cleaned.csv')
print(f" Loaded {len(df)} galaxies")
# ------------------------------------------------------------
# STEP 3: Convert to Supergalactic Cartesian coordinates
# ------------------------------------------------------------
print("\n3. Converting to Supergalactic coordinates...")
# Convert RA, Dec to Supergalactic coordinates
# Step 1: RA, Dec → ICRS Cartesian
coords_icrs = SkyCoord(ra=df['ra'].values * u.deg,
dec=df['dec'].values * u.deg,
frame='icrs')
# Step 2: ICRS → Supergalactic
coords_sg_coord = coords_icrs.transform_to('supergalactic')
df['sg_l'] = coords_sg_coord.l.deg
df['sg_b'] = coords_sg_coord.b.deg
# Convert distance to Mpc/h (Cosmicflows uses h=1, assume h=0.7)
d_comoving = df['distance_mpc'].values / 0.7
l_rad = np.radians(df['sg_l'].values)
b_rad = np.radians(df['sg_b'].values)
df['SGX'] = d_comoving * np.cos(b_rad) * np.cos(l_rad)
df['SGY'] = d_comoving * np.cos(b_rad) * np.sin(l_rad)
df['SGZ'] = d_comoving * np.sin(b_rad)
print(f" SGX range: {df['SGX'].min():.1f} to {df['SGX'].max():.1f} Mpc/h")
print(f" SGY range: {df['SGY'].min():.1f} to {df['SGY'].max():.1f} Mpc/h")
print(f" SGZ range: {df['SGZ'].min():.1f} to {df['SGZ'].max():.1f} Mpc/h")
# ------------------------------------------------------------
# STEP 4: Interpolate CF4 density
# ------------------------------------------------------------
print("\n4. Interpolating CF4 density...")
points = df[['SGX', 'SGY', 'SGZ']].values
df['cf4_density'] = interp_density(points)
# Clip to reasonable range
df['cf4_density'] = np.clip(df['cf4_density'], 0.1, 10.0)
print(f" Density range: {df['cf4_density'].min():.4f} to {df['cf4_density'].max():.4f}")
print(f" Median density: {df['cf4_density'].median():.4f}")
# ------------------------------------------------------------
# STEP 5: Assign extreme environment bins
# ------------------------------------------------------------
print("\n5. Assigning void vs cluster bins...")
df_sorted = df.sort_values('cf4_density')
n_extreme = len(df) // 3
df_void = df_sorted.head(n_extreme).copy()
df_cluster = df_sorted.tail(n_extreme).copy()
print(f" Void sample: {len(df_void)} galaxies (lowest density)")
print(f" Cluster sample: {len(df_cluster)} galaxies (highest density)")
print("\n 🌌 VOID galaxies:")
for _, row in df_void.iterrows():
print(f" {row['galaxy']:15s} γ={row['gamma']:.4f} δ={row['cf4_density']:.4f}")
print("\n 🌠 CLUSTER galaxies:")
for _, row in df_cluster.iterrows():
print(f" {row['galaxy']:15s} γ={row['gamma']:.4f} δ={row['cf4_density']:.4f}")
# ------------------------------------------------------------
# STEP 6: Statistical tests
# ------------------------------------------------------------
print("\n6. Running statistical tests...")
gamma_void = df_void['gamma'].values
gamma_cluster = df_cluster['gamma'].values
t_stat, p_tt = ttest_ind(gamma_void, gamma_cluster, equal_var=False)
ks_stat, p_ks = ks_2samp(gamma_void, gamma_cluster)
delta_gamma = np.median(gamma_cluster) - np.median(gamma_void)
print(f"\n Median γ (Void): {np.median(gamma_void):.4f}")
print(f" Median γ (Cluster): {np.median(gamma_cluster):.4f}")
print(f" Δγ (Cluster - Void): {delta_gamma:.4f}")
print(f"\n Welch t-test p-value: {p_tt:.5f}")
print(f" KS test p-value: {p_ks:.5f}")
# ------------------------------------------------------------
# STEP 7: Binary detection decision
# ------------------------------------------------------------
DETECTION_THRESHOLD = 0.03
ALPHA = 0.05
detected = (p_tt < ALPHA) and (p_ks < ALPHA)
print("\n" + "=" * 60)
print("FINAL PHASE 2 RESULT — CF4 DENSITY")
print("=" * 60)
if detected and abs(delta_gamma) >= DETECTION_THRESHOLD:
print("✅ DETECTION: Environmental constitutive drift detected")
print(f" Δγ = {delta_gamma:.4f} ({100*abs(delta_gamma)/np.median(gamma_void):.1f}% relative)")
print(" The substrate's constitutive response depends on environment")
elif detected:
print("⚠️ STATISTICAL SIGNIFICANCE BUT EFFECT BELOW THRESHOLD")
print(f" Δγ = {delta_gamma:.4f} < {DETECTION_THRESHOLD}")
else:
print("❌ NO DETECTION: No statistically significant drift")
print(f" Δγ = {delta_gamma:.4f} < {DETECTION_THRESHOLD}")
# ------------------------------------------------------------
# STEP 8: Visualization
# ------------------------------------------------------------
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].hist(gamma_void, bins=5, alpha=0.6, label='Void', color='blue', edgecolor='black')
axes[0].hist(gamma_cluster, bins=5, alpha=0.6, label='Cluster', color='red', edgecolor='black')
axes[0].axvline(np.median(gamma_void), color='blue', linestyle='--', linewidth=2)
axes[0].axvline(np.median(gamma_cluster), color='red', linestyle='--', linewidth=2)
axes[0].set_xlabel('γ')
axes[0].set_ylabel('Count')
axes[0].set_title(f'γ Distribution\nΔγ = {delta_gamma:.3f}')
axes[0].legend()
axes[1].boxplot([gamma_void, gamma_cluster], labels=['Void', 'Cluster'])
axes[1].set_ylabel('γ')
axes[1].set_title(f'Box Comparison\np_tt = {p_tt:.4f}, p_ks = {p_ks:.4f}')
axes[1].grid(alpha=0.3)
axes[2].scatter(df['cf4_density'], df['gamma'], c=df['gamma'], cmap='viridis', s=100, edgecolors='black')
axes[2].axhline(np.median(gamma_void), color='blue', linestyle='--', label=f'Void median γ = {np.median(gamma_void):.3f}')
axes[2].axhline(np.median(gamma_cluster), color='red', linestyle='--', label=f'Cluster median γ = {np.median(gamma_cluster):.3f}')
axes[2].set_xlabel('CF4 Density (1+δ)')
axes[2].set_ylabel('γ')
axes[2].set_title('γ vs CF4 Density')
axes[2].legend()
axes[2].grid(alpha=0.3)
plt.colorbar(axes[2].collections[0], ax=axes[2], label='γ')
plt.tight_layout()
plt.show()
print("\n✅ Phase 2 complete. This is your FINAL, PUBLISHABLE result using CF4.")
print(" Cite: Courtois et al. 2023 (A&A 670L 15) for the density field.")
