#!/usr/bin/env python3
"""Independent checks for DJAMILAR manuscript, 2026-09-23.

These calculations test stated reduced equations. They do NOT certify the
covariant HD constraints, a PPN solution, a lattice model, or observational fits.
Requirements: Python 3, NumPy, SciPy. Run from any directory.
"""
from pathlib import Path
import json
import math
import numpy as np
import scipy
from scipy.integrate import solve_bvp, quad

OUT = Path(__file__).resolve().parent
results = {}
results['software_versions'] = {'numpy': np.__version__, 'scipy': scipy.__version__}

# SI constants: numerical diagnostics, not a fit of the microscopic gap.
c = 299792458.0
hbar = 1.054571817e-34
G = 6.67430e-11
eV = 1.602176634e-19
pc = 3.085677581491367e16
GM_sun = 1.32712440018e20
rd = 25e-9
planck_energy = math.sqrt(hbar*c**5/G)
results['normalization'] = {
    'aperture_gap_eV_at_25_nm': hbar*c/rd/eV,
    'deep_gap_eV_for_R_AZ_1_and_C_IR_1': planck_energy/math.sqrt(16*math.pi)/eV,
    'hierarchy_for_R_AZ_1_and_C_IR_1': rd*planck_energy/(math.sqrt(16*math.pi)*hbar*c),
    'QCD_diagnostic_R_AZ': 0.0095/1.7**2,
}

# Explicit counterexample to inferring |alpha2| >= f from alpha2+2xi=-A.
f, KB = 2e-6, 2e-8
A = f + KB*(1-f)/2
results['ppn_algebra_counterexample'] = dict(f=f, K_B=KB, A=A,
    alpha1=-8*A, alpha2=0.0, xi=-A/2,
    interpretation='Algebraic counterexample only; not a covariant solution.')

# The quartic polynomial r^4 H4 is harmonic: Delta sum x_i^4 = 12 r^2,
# Delta (r^2)^2 = 20 r^2. Its zero is checked exactly at coefficient level.
results['cubic_harmonic'] = {
    'laplacian_coefficient': 12 - (3/5)*20,
    'axis': 1-3/5, 'face_diagonal': 1/2-3/5,
    'body_diagonal': 1/3-3/5,
}

# Gaussian momentum-transfer normalization and second moment in dimensionless k.
gauss_norm = 4/math.sqrt(math.pi)*quad(lambda z:z*z*math.exp(-z*z),0,np.inf)[0]
gauss_k2 = 4/math.sqrt(math.pi)*quad(lambda z:z**4*math.exp(-z*z),0,np.inf)[0]/3
results['localization'] = {'gaussian_normalization':gauss_norm,
    'one_axis_dimensionless_second_moment':gauss_k2,
    'expected_second_moment':0.5,
    'D_force_over_D_pp':2.0,
    'heating_coefficient_in_hbar2_Gamma_over_M_rd2':0.75}

# Off-shell check of the FLRW memory energy identity on arbitrary instantaneous
# backgrounds. No equation of motion is substituted; the identity must factor.
rng=np.random.default_rng(430919)
# Independent contraction of the linearized Ricci tensor. In Fourier space,
# d_mu d_nu -> -k_mu k_nu, h_00=-2 Phi, h_ii=-2 Psi.
eta=np.array([-1.,1.,1.,1.])
ricci_error=0.
for _ in range(100):
    k=rng.normal(size=4)
    Phi,Psi=rng.normal(size=2)
    h=np.diag([-2*Phi,-2*Psi,-2*Psi,-2*Psi])
    trace=np.sum(eta*np.diag(h))
    box_factor=-np.sum(eta*k*k)
    R00=.5*(-2*np.sum(eta*k*k[0]*h[:,0])-box_factor*h[0,0]+k[0]**2*trace)
    expected=-np.sum(k[1:]**2)*Phi-3*k[0]**2*Psi
    ricci_error=max(ricci_error,abs(R00-expected)/(1+abs(R00)+abs(expected)))
results['linearized_R00_max_scaled_residual']=ricci_error

max_res=0.
for _ in range(100):
    H, Hd, v, vd, Vp, gamma = rng.normal(size=6)
    rho_dot=v*vd+Vp*v+3*gamma*(Hd*v+H*vd)
    rho_plus_p=v*v+3*gamma*H*v-gamma*vd
    lhs=rho_dot+3*H*rho_plus_p
    rhs=v*(vd+3*H*v+Vp+3*gamma*(Hd+3*H*H))
    max_res=max(max_res,abs(lhs-rhs)/(1+abs(lhs)+abs(rhs)))
results['FLRW_off_shell_identity_max_scaled_residual']=max_res

# Exact homogeneous temporal solution K=K2(Q-Q0)^2 and charge a^3 K_Q=C.
K2,Q0,C=1.3,0.7,0.19
charge_errors=[]
energy_errors=[]
for a in np.geomspace(.3,3,50):
    Q=Q0+C/(2*K2*a**3)
    K=K2*(Q-Q0)**2
    KQ=2*K2*(Q-Q0)
    charge_errors.append(abs(a**3*KQ-C))
    energy_errors.append(abs(Q*KQ-K-(Q0*C/a**3+C*C/(4*K2*a**6))))
results['FLRW_temporal_exact_solution']={
    'charge_max_absolute_error':max(charge_errors),
    'energy_max_absolute_error':max(energy_errors)}

# BVP for the NEW explicitly normalized static deep-response diagnostic:
# f''+2 f'/y-2 f/y^2 = f |f| - 1/y^2, y=r/R_s.
# A point source is used. The inner regular branch has
# f=1/2+b*y+y^2/16+O(y^3); the far boundary uses its asymptotic expansion.
# Repeating once on a larger interval checks a concrete boundary-truncation risk.
def solve_profile(ymin,ymax,tol):
    # t=ln(y) avoids subtracting separate singular 1/y^2 terms near the origin.
    mesh=np.linspace(np.log(ymin),np.log(ymax),900)
    yy=np.exp(mesh)
    initial=np.vstack((1/(yy+2),-yy/(yy+2)**2))
    def ode(t,z):
        return np.vstack((z[1],np.exp(2*t)*z[0]*np.abs(z[0])-1-z[1]+2*z[0]))
    far=1/ymax-1/ymax**2-0.5/ymax**3-1.5/ymax**4
    def bc(za,zb):
        return np.array([za[0]-za[1]-0.5+ymin**2/16,zb[0]-far])
    sol=solve_bvp(ode,bc,mesh,initial,tol=tol,max_nodes=30000)
    if not sol.success:
        raise RuntimeError(sol.message)
    # Residual on an independent grid using derivatives of the returned spline.
    grid=np.linspace(np.log(ymin*1.03),np.log(ymax/1.03),1400)
    z=sol.sol(grid)
    dz=sol.sol(grid,1)
    rhs=ode(grid,z)
    err=np.max(np.abs(dz-rhs)/(1+np.abs(rhs)))
    return sol, float(err), float(np.max(sol.rms_residuals))

sol1,e1,r1=solve_profile(1e-4,1e3,2e-7)
sol2,e2,r2=solve_profile(2e-5,2e3,5e-8)
check_y=np.geomspace(1e-3,300,250)
v1=sol1.sol(np.log(check_y))[0]; v2=sol2.sol(np.log(check_y))[0]
sample_y=np.array([1e-3,1e-2,0.1,1.,10.,100.,1000.])
sample_f=sol2.sol(np.log(sample_y))[0]
results['static_HD_BVP']={
    'status':'Numerical witness of a reduced static equation only',
    'first_domain':[1e-4,1e3], 'second_domain':[2e-5,2e3],
    'requested_solver_tolerances':[2e-7,5e-8],
    'nodes':[int(sol1.x.size),int(sol2.x.size)],
    'max_collocation_RMS_residual':[r1,r2],
    'independent_spline_scaled_residual':[e1,e2],
    'max_relative_profile_difference_on_overlap':float(np.max(abs(v1-v2)/v2)),
    'min_f_on_overlap':float(np.min(v2)),
    'max_logarithmic_derivative_on_overlap':float(np.max(sol2.sol(np.log(check_y))[1])),
    'sample':[dict(y=float(y),f=float(fv),g_over_deep=float(y*fv)) for y,fv in zip(sample_y,sample_f)]}

# Finite uniform sphere: R=1, C/L^2=1 and the homogeneous A_match r
# contribution omitted because the radial operator annihilates it.
u=np.geomspace(1e-3,1,100)
gin=.5*(u-u**3/5)
ginp=.5*(1-3*u*u/5)
ginpp=-3*u/5
rin=ginpp+2*ginp/u-2*gin/u**2+u
uout=np.geomspace(1,100,100)
gout=.5*(1-1/(5*uout*uout))
goutp=1/(5*uout**3)
goutpp=-3/(5*uout**4)
rout=goutpp+2*goutp/uout-2*gout/uout**2+1/uout**2
matching=np.array([gin[-1]-gout[0],ginp[-1]-goutp[0],ginpp[-1]-goutpp[0]])
results['finite_uniform_sphere']={
    'interior_operator_max_absolute_residual':float(np.max(abs(rin))),
    'exterior_operator_max_absolute_residual':float(np.max(abs(rout))),
    'boundary_matching_max_absolute_residual':float(np.max(abs(matching))),
    'scope':'Leading HD-dominated equation only; outer matching coefficient remains fixed by the full solution.'}

# Scale separation estimates use L (the static coefficient), NOT the raw
# covariant ell_theta; L^2=ell_theta^2/(2-K_B).
L=10*pc; aI=1.2e-10
rows=[]
for solar_mass in [1,1e8,1e10,1e11]:
    Cgrav=GM_sun*solar_mass
    rM=math.sqrt(Cgrav/aI)
    Rs=L*L/rM
    rows.append({'mass_in_solar_masses':solar_mass,'r_M_pc':rM/pc,
        'R_s_pc':Rs/pc,'inner_acceleration_m_s2':Cgrav/(2*L*L),
        'leading_R_s_over_r_at_1_kpc':Rs/(1e3*pc)})
results['screening_scale_diagnostics']=rows
results['screening_scale_diagnostics_caution']=(
    'R_s/r is an expansion parameter, not a calculated force correction when it is large. '
    'The physical deep-response approximation also requires g_chi << a_I and appropriate environmental matching.')

# Lensing integral for constant metric slip. The full metric outside a finite
# MOND interval still needs matching; this tests the idealized analytic factor.
b=2.1
integral=quad(lambda z:b/(b*b+z*z),-np.inf,np.inf)[0]
results['lensing_integral_pi']={'numerical':integral,'pi':math.pi,
    'absolute_error':abs(integral-math.pi)}

assert abs(gauss_norm-1)<1e-12 and abs(gauss_k2-.5)<1e-12
assert max_res<1e-13
assert ricci_error<1e-13
assert max(charge_errors)<1e-12 and max(energy_errors)<1e-12
assert np.all(v2>0) and np.all(sol2.sol(np.log(check_y))[1]<0)
assert results['static_HD_BVP']['max_relative_profile_difference_on_overlap']<1e-5
assert abs(integral-math.pi)<1e-12
assert np.max(abs(rin))<1e-10 and np.max(abs(rout))<1e-12
assert np.max(abs(matching))<1e-12
results['scope_limits']=[
    'No proof of a complete microscopic BCC/BGG Hamiltonian.',
    'No full covariant HD or memory Dirac analysis.',
    'No derivation of PPN coefficients on a screened finite source.',
    'No observational fit, galaxy sample reconstruction, or cosmological spectrum.',
    'No proof of vacuum-energy sequestering for the compact action.']

# V5.2: independent scalar checks of the new cross-sector relations.
from scipy.optimize import minimize_scalar
rng_new=np.random.default_rng(20260920)
identity_errors=[]
for beta, kb in zip(rng_new.uniform(0,2,100),rng_new.uniform(.001,1.99,100)):
    alpha=-8*(beta+kb/2)/(1+beta)
    cir=(1+beta)/(1-kb/2)
    identity_errors.append(abs(cir-1/(1+alpha/8))/cir)
results['tracking_identity']={'max_relative_error':max(identity_errors),
    'upper_deviation_at_epsilon_1e4':1e-4/(8-1e-4)}
Omega=.26; K2=1.0; CIR=1.0
amin=math.sqrt(3*Omega/(2*K2*CIR))
pressure_rows=[]
for qmax in [.001,.01,.1,.5,1.0]:
    def drift(logq):
        q=math.exp(logq)
        return math.sqrt(3*Omega/(K2*CIR))*(1+2*q)/(4*math.sqrt(q))
    opt=minimize_scalar(drift,bounds=(-25,math.log(qmax)),method='bounded',options={'xatol':1e-13})
    analytic=math.sqrt(3*Omega/(K2*CIR))*((1+2*qmax)/(4*math.sqrt(qmax)) if qmax<.5 else 1/math.sqrt(2))
    pressure_rows.append({'qmax':qmax,'analytic_floor':analytic,'numerical_floor':opt.fun,
        'relative_difference':abs(opt.fun-analytic)/analytic})
H0=70*1000/(1e6*pc); year=365.25*86400
results['pressure_drift']={'unrestricted_floor':amin,'saturating_q':.5,'saturating_w':1/3,
    'homogeneous_floor_per_year_at_H0_70':amin*H0*year,'restricted':pressure_rows}
AU=149597870700.; gp=GM_sun/(9.58*AU)**2; gbar=1.2e-10; delta=1e-14
excess=gbar/math.expm1(1)
lpair=(excess-delta)/(gp-gbar)
results['single_pair_ppn_gate']={'galaxy_excess':excess,'saturn_baryonic_acceleration':gp,
    'L_pair':lpair,'epsilon1_critical':8*lpair,
    'exact_boundary_gp':gbar+8*(excess-delta)/(8*lpair),
    'historical_stress_not_joint_likelihood':True}
results['screening_additions']={'crossover_mass_solar_at_L10pc':aI*L**2/GM_sun,
    'historical_Mars_stress_L_pc':math.sqrt(GM_sun/(2e-15))/pc,
    'deep_force_fraction_at_12a0':1/math.sqrt(12),
    'RAR_central_excess_fraction_at_12a0':1/math.expm1(math.sqrt(12))}
# Fourier convention exp(i k.x - i omega t): compute divergence directly.
kappa_errors=[]
for _ in range(100):
    kval=rng_new.normal(size=3); ksq=kval@kval; om=rng_new.normal(); kap=rng_new.uniform(-2,3)
    phi,psi,alpha=rng_new.normal(size=3)
    da=-ksq*(phi-1j*om*alpha)
    expansion=-ksq*alpha+3j*om*psi
    direct=da-kap*(-1j*om)*expansion
    formula=-ksq*phi-3*kap*om**2*psi+1j*om*ksq*(1-kap)*alpha
    kappa_errors.append(abs(direct-formula))
results['memory_kappa_linearization_max_error']=max(kappa_errors)
# Independent finite differences of the deep constitutive flux F(v)=|v|v.
v=np.array([.4,-.7,1.1]); norm=np.linalg.norm(v); h=1e-5
jac=np.column_stack([(np.linalg.norm(v+h*np.eye(3)[j])*(v+h*np.eye(3)[j])-
    np.linalg.norm(v-h*np.eye(3)[j])*(v-h*np.eye(3)[j]))/(2*h) for j in range(3)])
exact=norm*np.eye(3)+np.outer(v,v)/norm
results['external_field_jacobian']={'finite_difference_max_error':float(np.max(abs(jac-exact))),
    'eigenvalues':np.linalg.eigvalsh(exact).tolist(),'expected':[norm,norm,2*norm]}
# Poisson generating function and the operational visibility threshold.
mu=2.; dsep=.8; rtest=.7; nu=math.exp(-dsep*dsep/(4*rtest*rtest))
psum=sum(math.exp(-mu)*mu**n/math.factorial(n)*nu**n for n in range(100))
closed=math.exp(-mu*(1-nu)); eta=.4
threshold=2*rtest*math.sqrt(-math.log(1-(-math.log(eta))/mu))
recovered=math.exp(-mu*(1-math.exp(-threshold**2/(4*rtest**2))))
results['quantum_channel_checks']={'poisson_sum_error':abs(psum-closed),
    'visibility_threshold_recovery_error':abs(recovered-eta)}
assert max(identity_errors)<1e-12
assert max(x['relative_difference'] for x in pressure_rows)<1e-6
assert abs(results['single_pair_ppn_gate']['exact_boundary_gp']-gp)<1e-18
assert max(kappa_errors)<1e-12
assert np.max(abs(jac-exact))<1e-8
assert abs(psum-closed)<1e-14 and abs(recovered-eta)<1e-14
results['scope_limits'] += [
    'No microscopic existence proof for Djamilars or their common quantum dynamics.',
    'No laboratory transfer prediction, calculated shift-breaking source or measured drift.',
    'No validation of the original hard size threshold, inter-universal force or black-hole saturation law.']

target=OUT/'DJAMILAR_verification.json'
target.write_text(json.dumps(results,indent=2,ensure_ascii=False)+'\n')
print(json.dumps(results,indent=2,ensure_ascii=False))