import PhysicsLogic.Papers.KaplanSavageWise1998

open PhysicsLogic.Papers.KaplanSavageWise1998 Filter
open PhysicsLogic.Framework.EffectiveRangeExpansion
open scoped Topology

set_option backward.isDefEq.respectTransparency false

-- An algebraic unitary-limit EFT model, not measured nuclear data or a
-- finite-range potential realization. Units: M=4*pi, r0=1, Lambda=2.
private noncomputable def unitaryFlow : PionlessRGFlow where
  mass := 4 * Real.pi
  inverseScatteringLength := 0
  effectiveRange := 1
  c0 := fun mu => -mu⁻¹
  c2 := fun mu => 1 / (2 * mu ^ 2)

private theorem observable :
    HasFirstTwoCoefficients (fun q : Real => q ^ 2 / 2) 0 (1 / 2) := by
  have positive : ∀ᶠ q in approach, 0 < q := self_mem_nhdsWithin
  apply tendsto_const_nhds.congr'
  filter_upwards [positive] with q hq
  field_simp [ne_of_gt hq]
  ring

private theorem couplings (mu : Real) (nonpole : -mu + (0 : Real) ≠ 0) :
    HasFirstTwoCoefficients (fun q => -1 / (q ^ 2 / 2 + mu))
      (unitaryFlow.mass / (4 * Real.pi) * unitaryFlow.c0 mu)
      (unitaryFlow.mass / (4 * Real.pi) * unitaryFlow.c2 mu) := by
  apply (ContactMatching.canonical_contact_coefficients_iff_values
    (fun q => q ^ 2 / 2) mu 0 (1 / 2) _ _
    (by simpa using observable) (by simpa using nonpole)).mpr
  have hn : mu ≠ 0 := by simpa using nonpole
  constructor <;> simp [unitaryFlow, Real.pi_ne_zero]
  <;> field_simp

private theorem local_equations (p : Real) (positive : 0 < p) :
    let g0 := unitaryFlow.mass / (4 * Real.pi) * unitaryFlow.c0 p
    let g2 := unitaryFlow.mass / (4 * Real.pi) * unitaryFlow.c2 p
    Complex.I / (p : Complex) + (g0 : Complex) *
        ((p : Complex) + Complex.I * (p : Complex)) * (Complex.I / (p : Complex)) =
      -(g0 : Complex) ∧
    (1 / 2 : Complex) + (g0 : Complex) * ((p : Complex) + Complex.I * (p : Complex)) *
        (1 / 2 : Complex) = -((g2 * p ^ 2 : Real) : Complex) *
          (1 + ((p : Complex) + Complex.I * (p : Complex)) * (Complex.I / (p : Complex))) := by
  have hn := ne_of_gt positive
  have hc : (p : Complex) ≠ 0 := Complex.ofReal_ne_zero.mpr hn
  dsimp [unitaryFlow]
  norm_num [Real.pi_ne_zero]
  constructor <;> push_cast <;> field_simp [hn, hc] <;> ring_nf
  <;> simp [Complex.I_sq]

-- This consumes the complete first-order contract, including actual running,
-- both optical relations and every finite-loop error, not only a smallness
-- inequality. The leading raw series fails while the derivative correction
-- becomes relatively small as p decreases. Its unknown physical remainder
-- remains outside the contract.
example (p : Real) (positive : 0 < p) (small : p < 1) :
    (forall mu : Real, mu ≠ 0 ->
      HasDerivAt (fun x : Real => -x⁻¹) (mu⁻¹ ^ 2) mu ∧
      HasDerivAt (fun x : Real => 1 / (2 * x ^ 2)) (-1 / mu ^ 3) mu) ∧
    ‖(1 / 2 : Complex)‖ = p / 2 * ‖Complex.I / (p : Complex)‖ ∧
    ‖(1 / 2 : Complex)‖ < ‖Complex.I / (p : Complex)‖ ∧
    Complex.normSq (1 + 2 * Complex.I * (p : Complex) * (Complex.I / (p : Complex))) = 1 ∧
    (1 / 2 : Complex).im = 2 * p * (star (Complex.I / (p : Complex)) * (1 / 2 : Complex)).re ∧
    (forall t : Real,
      (Complex.I / (p : Complex) + (t : Complex) * (1 / 2 : Complex)).im -
        p * Complex.normSq (Complex.I / (p : Complex) + (t : Complex) * (1 / 2 : Complex)) =
          -p * t ^ 2 / 4) ∧
    (forall n : Nat,
      let loopSum := (Finset.range (n + 1)).sum (fun j =>
        -((-p⁻¹ : Real) : Complex) * (-((-p⁻¹ : Real) : Complex) *
          ((p : Complex) + Complex.I * (p : Complex))) ^ j)
      ‖loopSum - Complex.I / (p : Complex)‖ =
        ‖(1 : Complex) + Complex.I‖ ^ (n + 1) / p ∧
      ‖Complex.I / (p : Complex)‖ <= ‖loopSum - Complex.I / (p : Complex)‖) := by
  have hn := ne_of_gt positive
  have complete := (PhysicsLogic.Papers.KaplanSavageWise1998.PionlessRGFlow.observable_matching_controls_the_complete_first_order_argument :
    forall (data : PionlessRGFlow), _)
    unitaryFlow (fun q => q ^ 2 / 2) p 2 1 (by dsimp [unitaryFlow]; positivity)
    (by simpa [unitaryFlow, HasFirstTwoCoefficients, approach] using observable)
    couplings positive (by norm_num) (by norm_num [unitaryFlow])
  have point := complete.2 p (Complex.I / (p : Complex)) (1 / 2)
    (by simpa [unitaryFlow] using hn) (local_equations p positive)
  rcases point with ⟨_, _, _, strict, unitary, linearized, defect, errors, _⟩
  have leadingNorm : ‖Complex.I / (p : Complex)‖ = 1 / p := by
    simp [norm_div, Complex.norm_real, Real.norm_eq_abs, abs_of_pos positive]
  have ratio : |unitaryFlow.mass / (4 * Real.pi) * unitaryFlow.c0 p| *
      ‖(p : Complex) + Complex.I * (p : Complex)‖ = ‖(1 : Complex) + Complex.I‖ := by
    have factor : (p : Complex) + Complex.I * (p : Complex) =
        (p : Complex) * (1 + Complex.I) := by ring
    rw [factor, norm_mul]
    simp [unitaryFlow, Real.pi_ne_zero, abs_inv, abs_of_pos positive,
      Complex.norm_real, Real.norm_eq_abs]
    field_simp
  have ratioAboveOne : 1 <= ‖(1 : Complex) + Complex.I‖ := by
    simpa using Complex.abs_im_le_norm ((1 : Complex) + Complex.I)
  refine ⟨?_, ?_, strict (by linarith), unitary, linearized, ?_, ?_⟩
  · intro mu hm
    have running := complete.1 mu (by simpa [unitaryFlow] using hm)
    constructor
    · simpa [unitaryFlow, Real.pi_ne_zero] using running.2.2.1
    · convert running.2.2.2 using 1 <;> simp [unitaryFlow, Real.pi_ne_zero]
      field_simp
      ring
  · rw [leadingNorm]
    norm_num
    field_simp
  · intro t
    convert defect t using 1 <;> norm_num [Complex.normSq_apply] <;> ring
  · intro n
    have error := errors n
    rw [ratio] at error
    refine ⟨?_, ?_⟩
    · simpa [unitaryFlow, Real.pi_ne_zero, leadingNorm, div_eq_mul_inv,
        abs_of_pos positive] using error.1
    · simpa [unitaryFlow, Real.pi_ne_zero] using error.2 ratioAboveOne

-- The same ERE/contact relation permits a zero effective range. Therefore
-- an upper-order assignment cannot supply a nonzero C2 leading coefficient.
example (mu : Real) (nonpole : mu ≠ 0) :
    HasFirstTwoCoefficients (fun q : Real => -1 / (0 + mu)) (-1 / mu) 0 ∧
      ¬(exists lower : Real, 0 < lower ∧ forall p : Real, 0 < p -> lower <= |p ^ 2 * 0|) := by
  constructor
  · exact (ContactMatching.canonical_contact_coefficients_iff_values
      (fun _ => 0) mu 0 0 (-1 / mu) 0
      (by simpa [HasFirstTwoCoefficients] using
        (tendsto_const_nhds : Tendsto (fun _ : Real => 0) approach (nhds (0 : Real))))
      (by simpa using nonpole)).mpr
        (by constructor <;> simp [div_eq_mul_inv])
  · rintro ⟨lower, positive, bound⟩
    have impossible := bound 1 (by norm_num)
    norm_num at impossible
    linarith

#print axioms PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.canonical_contact_coefficients_iff_values
#print axioms PhysicsLogic.Papers.KaplanSavageWise1998.PionlessRGFlow.observable_matching_controls_the_complete_first_order_argument