import PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching

set_option backward.isDefEq.respectTransparency false

open Filter PhysicsLogic.Framework.EffectiveRangeExpansion
open scoped Topology

example : forall (phase : Real -> Real) (constant quadratic : Real),
    (Tendsto phase (nhdsWithin 0 (Set.Ioi 0)) (𝓝 constant) ∧
      Tendsto (fun q => (phase q - constant) / q ^ 2)
        (nhdsWithin 0 (Set.Ioi 0)) (𝓝 quadratic)) ↔
    Tendsto (fun q => (phase q - constant) / q ^ 2)
      (nhdsWithin 0 (Set.Ioi 0)) (𝓝 quadratic) :=
  PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.first_two_coefficients_iff_quadratic_limit

#print axioms PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.first_two_coefficients_iff_quadratic_limit

example : forall (phase : Real -> Real) (mu g0 g1 : Real),
    (Tendsto (fun q => (-1 / (phase q + mu) - g0) / q ^ 2)
        (nhdsWithin 0 (Set.Ioi 0)) (nhds g1) ∧
      ∀ᶠ q in nhdsWithin 0 (Set.Ioi 0), phase q + mu ≠ 0) ↔
    ∃ contact : Real -> Real,
      Tendsto (fun q => (contact q - g0) / q ^ 2)
        (nhdsWithin 0 (Set.Ioi 0)) (nhds g1) ∧
      ∀ᶠ q in nhdsWithin 0 (Set.Ioi 0), (phase q + mu) * contact q = -1 :=
  PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.inverse_amplitude_contact_coefficients_iff_germ

#print axioms PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.inverse_amplitude_contact_coefficients_iff_germ

example : ∀ (phase contact : Real -> Real) (mu : Real),
    (∀ᶠ p in nhdsWithin 0 (Set.Ioi 0), (phase p + mu) * contact p = -1) ->
    (∀ᶠ p in nhdsWithin 0 (Set.Ioi 0), phase p + mu ≠ 0) ∧
      (∀ᶠ p in nhdsWithin 0 (Set.Ioi 0), contact p = -1 / (phase p + mu)) :=
  PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.inverse_amplitude_identity_fixes_contact_germ

-- Near-threshold matching fixes neither the contact at zero nor its value
-- at a distant prediction momentum. Both witnesses are explicit numbers.
example :
    let contact : Real -> Real := fun p => if 0 < p ∧ p < 1 then 1 else 42
    (∀ᶠ p in approach, ((-1 : Real) + 0) * contact p = -1) ∧
      contact 0 ≠ -1 / ((-1 : Real) + 0) ∧
      contact 2 ≠ -1 / ((-1 : Real) + 0) := by
  dsimp only
  refine ⟨?_, by norm_num, by norm_num⟩
  have positive : ∀ᶠ p in approach, 0 < p := self_mem_nhdsWithin
  have below : ∀ᶠ p in approach, p < 1 :=
    the_approach_really_tends_to_zero.eventually (gt_mem_nhds (by norm_num))
  filter_upwards [positive, below] with p hp hb
  norm_num [hp, hb]

-- A fixed-scale threshold neighborhood need not be uniform as the
-- subtraction scale approaches the contact pole, even for a polynomial K.
example :
    (∀ mu : Real, mu ≠ 1 -> ∀ᶠ p in approach, -1 + p ^ 2 + mu ≠ 0) ∧
      ¬∃ radius : Real, 0 < radius ∧
        ∀ mu : Real, mu ≠ 1 -> ∀ p : Real, 0 < p -> p < radius -> -1 + p ^ 2 + mu ≠ 0 := by
  constructor
  · intro mu nonpole
    have limit : Tendsto (fun p : Real => -1 + p ^ 2 + mu) approach (𝓝 (-1 + mu)) := by
      simpa using ((the_approach_really_tends_to_zero.pow 2).const_add (-1)).add_const mu
    exact limit.eventually (eventually_ne_nhds (by intro h; apply nonpole; linarith))
  · rintro ⟨radius, positive, uniform⟩
    let p := min (radius / 2) 1
    have hp : 0 < p := lt_min (by linarith) (by norm_num)
    have below : p < radius := (min_le_left _ _).trans_lt (by linarith)
    have nonpole : 1 - p ^ 2 ≠ 1 := by nlinarith [sq_pos_of_pos hp]
    exact uniform (1 - p ^ 2) nonpole p hp below (by ring)

example : ∀ (phase remainder : Real -> Real) (inverseLength rangeCoefficient : Real),
    (∀ᶠ p in approach,
      phase p = -inverseLength + rangeCoefficient * p ^ 2 + remainder p) ->
    Tendsto (fun p => remainder p / p ^ 2) approach (𝓝 0) ->
    HasFirstTwoCoefficients phase (-inverseLength) rangeCoefficient :=
  PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.negligible_remainder_identifies_observable_coefficients

example : ∀ (phase contact : Real -> Real) (mu inverseLength rangeCoefficient g0 g1 : Real),
    HasFirstTwoCoefficients phase (-inverseLength) rangeCoefficient ->
    HasFirstTwoCoefficients contact g0 g1 ->
    (∀ᶠ p in approach, (phase p + mu) * contact p = -1) ->
    (inverseLength - mu) * g0 = 1 ∧
      (inverseLength - mu) * g1 = rangeCoefficient * g0 :=
  PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.observable_coefficients_determine_contact_matching

-- A nonzero cubic remainder is enough for coefficient identification, but
-- does not amount to an analytic expansion in p^2 or a quartic error bound.
example : HasFirstTwoCoefficients (fun p : Real => -2 + 3 * p ^ 2 + p ^ 3) (-2) 3 := by
  apply PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.negligible_remainder_identifies_observable_coefficients
    _ (fun p => p ^ 3) 2 3 (Eventually.of_forall (fun _ => rfl))
  have equal : (fun p : Real => p ^ 3 / p ^ 2) =ᶠ[approach] (fun p => p) := by
    have positive : ∀ᶠ p in approach, 0 < p := self_mem_nhdsWithin
    filter_upwards [positive] with p hp
    field_simp [ne_of_gt hp]
  exact the_approach_really_tends_to_zero.congr' equal.symm

-- An unrestricted remainder can instead hide a wrong quadratic coefficient.
example : (∀ p : Real, -2 + 3 * p ^ 2 = -2 + 7 * p ^ 2 + (-4 * p ^ 2)) ∧
    ¬Tendsto (fun p : Real => (-4 * p ^ 2) / p ^ 2) approach (𝓝 0) := by
  refine ⟨fun _ => by ring, ?_⟩
  intro vanishes
  letI : approach.NeBot := the_positive_energy_approach_is_nontrivial
  have equal : (fun p : Real => (-4 * p ^ 2) / p ^ 2) =ᶠ[approach] (fun _ => -4) := by
    have positive : ∀ᶠ p in approach, 0 < p := self_mem_nhdsWithin
    filter_upwards [positive] with p hp
    field_simp [ne_of_gt hp]
  have impossible : (0 : Real) = -4 :=
    tendsto_nhds_unique vanishes (tendsto_const_nhds.congr' equal.symm)
  norm_num at impossible

example : forall (phase : Real -> Real) (mu inverseLength rangeCoefficient g0 g1 : Real),
    Tendsto (fun q => (phase q - (-inverseLength)) / q ^ 2)
      (nhdsWithin 0 (Set.Ioi 0)) (nhds rangeCoefficient) ->
    inverseLength - mu ≠ 0 ->
    (Tendsto (fun q => (-1 / (phase q + mu) - g0) / q ^ 2)
        (nhdsWithin 0 (Set.Ioi 0)) (nhds g1) ↔
      g0 = 1 / (inverseLength - mu) ∧
        g1 = rangeCoefficient / (inverseLength - mu) ^ 2) :=
  PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.canonical_contact_coefficients_iff_values

#print axioms PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.canonical_contact_coefficients_iff_values
#print axioms PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.negligible_remainder_identifies_observable_coefficients
#print axioms PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.observable_coefficients_determine_contact_matching
#print axioms PhysicsLogic.Papers.KaplanSavageWise1998.ContactMatching.inverse_amplitude_identity_fixes_contact_germ