Elliptical parameters of rotated ellipse

Stream: new members

[ZhiKai Pong]

A general ellipse (centered at origin) can be characterized by two angles which gives useful information about the orientation and shape as shown here

Given a parametrization, I would like to show the results at the end:

import Mathlib

open Real

variable {a b τ δx δy : ℝ}

/- Parametrization. -/
local notation "x" => a * cos (τ + δx)
local notation "y" => b * cos (τ + δy)

/- Analytic equation of ellipse. -/
theorem parametrizationisEllipse (ha : a ≠ 0) (hb : b ≠ 0) :
    (x / a) ^ 2 + (y / b) ^ 2 - 2 * (x / a) * (y / b) * cos (δy - δx) = sin (δy - δx) ^ 2 := by
  have h1 : (a * cos (τ + δx)) / a * sin δy - (b * cos (τ + δy)) / b * sin δx =
      cos τ * sin (δy - δx) := by
    field_simp
    rw [cos_add, cos_add, sin_sub]
    ring
  have h2 : (a * cos (τ + δx)) / a * cos δy - (b * cos (τ + δy)) / b * cos δx =
      sin τ * sin (δy - δx) := by
    field_simp
    rw [cos_add, cos_add, sin_sub]
    ring
  trans (a * cos (τ + δx) / a) ^ 2 * (sin δy ^ 2 + cos δy ^ 2) +
      (b * cos (τ + δy) / b) ^ 2 * (sin δx ^ 2 + cos δx ^ 2) -
      2 * (a * cos (τ + δx) / a) * (b * cos (τ + δy) / b) * cos (δy - δx)
  · simp
  trans (cos τ * sin (δy - δx)) ^ 2 + (sin τ * sin (δy - δx)) ^ 2
  · rw [← h1, ← h2]
    rw [cos_sub]
    ring
  trans (cos τ ^ 2 + sin τ ^ 2) * sin (δy - δx) ^ 2
  · ring
  simp

variable {χ ψ : ℝ}

/- want to show the following: -/

lemma ellipticity_angle : sin (2 * χ) = sin (2 * arctan (b/a)) * sin (δy - δx) := by
  sorry

lemma orientation_angle : tan (2 * ψ) = a * b / (a ^ 2 - b ^ 2) * cos (δy - δx) := by
  sorry

The proof that I know goes something like this:

variable {χ ψ : ℝ}

variable {a' b' τ δ0 : ℝ}

/- Rotate such that the major and minor axis aligns with x' and y'. -/

local notation "x'" => x * cos ψ + y * sin ψ
local notation "y'" => -x * sin ψ + y * cos ψ

/- Assume this (need to know that rotated ellipse is still an ellipse?) -/

lemma x'_parametrization : x' = a' * cos (τ + δ0) := by
  sorry

lemma y'_parametrization : y' = b' * cos (τ + δ0) := by
  sorry

/- Equate the two definitions of `x'` and `y'` then simplify. -/

but I'm aware this is quite hand-wavy, in particular I have no idea how to justify the second step formally. I would like to try to state this result in lean but I'm not quite sure how to start, and if there's any examples or more generally how to state and prove things in analytic geometry I'd appreciate any advice and recommendations, thanks a lot :)

(I'm using local notation here because if I use def or abbrev I have to carry the parameters around which will look like x a τ δx making the expressions more messy than it already is, advice on how to best state those variables is also much appreciated)

[Kenny Lau]

how do you know that those angles give you the correct result?

[ZhiKai Pong]

Kenny Lau said:

how do you know that those angles give you the correct result?

I'm not really sure how to answer this other than "the maths work out".

/- x' y' are coordinates of rotated axes. -/
x' := x * cos ψ + y * sin ψ
y' := -x * sin ψ + y * cos ψ

/- a' and b' are the major and minor axes of the ellipse. -/
tan χ = ∓ (b' / a')

the angles are given by those relationships, then it's just applying trigonometric identities and simplifying. The question is more about how I should set the problem up and define things as currently I'm understanding this through the geometry where I can just "get" the angle from the diagram and I don't know how to state this in lean.

[Kenny Lau]

An algebraic geometer would probably say that this is a quadratic with a given point so any line passing through the point must intersect the curve at one more place and this parametrises the conic section

[Kenny Lau]

is it possible to reduce the question to first a standard orientation?

[ZhiKai Pong]

Kenny Lau said:

An algebraic geometer would probably say that this is a quadratic with a given point so any line passing through the point must intersect the curve at one more place and this parametrises the conic section

I don't think mathlib has anything on conic sections :(

Kenny Lau said:

is it possible to reduce the question to first a standard orientation?

The orientation angles are what I want in this case, the goal here is to show that those two angles characterizes all ellipses (orientation and ellipticity) up to a "normalization" with a ^ 2 + b ^ 2 = 1. (this is useful for the polarization ellipse)

[ZhiKai Pong]

sorry @Kenny Lau the title for the other thread was too generic but I'm unable to edit the title anymore without breaking it up, might be better if I just ask follow up questions here

[ZhiKai Pong]

posting the last message and full code here again:

I think I'm getting myself confused again so @Kenny Lau I hope you don't mind me asking for further advice.

Essentially I'm having trouble with the sorrys at the end:

code

import Mathlib

structure GeneralEllipse where
  A : ℝ
  B : ℝ
  C : ℝ
  F : ℝ
  h : B ^ 2 - 4 * A * C < 0 ∨ F = 0

namespace GeneralEllipse

open Real

structure EllipseParameters where
  E₀x : ℝ
  E₀y : ℝ
  δx : ℝ
  δy : ℝ

structure Point (E : GeneralEllipse) where
  x : ℝ
  y : ℝ
  eq : E.A * x ^ 2 + E.B * x * y + E.C * y ^ 2 + E.F = 0

namespace EllipseParameters

@[simps]
noncomputable def ellipse (p : EllipseParameters) : GeneralEllipse where
  A := p.E₀y ^ 2
  B := - 2 * p.E₀x * p.E₀y * cos (p.δy - p.δx)
  C := p.E₀x ^ 2
  F := - p.E₀x ^ 2 * p.E₀y ^ 2 * sin (p.δy - p.δx) ^ 2
  h := by
    have : 0 ≤ (p.E₀x * p.E₀y * sin (p.δy - p.δx)) ^ 2 := sq_nonneg _
    obtain h | h := eq_or_lt_of_le this
    · right
      calc
        _ = - (p.E₀x * p.E₀y * sin (p.δy - p.δx)) ^ 2 := by ring
        _ = 0 := by rw [← h]; ring
    · left
      calc
        _ = - 4 * p.E₀y ^ 2 * p.E₀x ^ 2 * (1 - cos (p.δy - p.δx) ^ 2) := by ring
        _ = - 4 * (p.E₀x * p.E₀y * sin (p.δy - p.δx)) ^ 2 := by rw [← sin_sq]; ring
        _ < 0 := by linarith

@[simps]
noncomputable def eval (p : EllipseParameters) (τ : ℝ) : p.ellipse.Point where
  x := p.E₀x * cos (τ + p.δx)
  y := p.E₀y * cos (τ + p.δy)
  eq := by
    simp only [EllipseParameters.ellipse_A, EllipseParameters.ellipse_B,
      EllipseParameters.ellipse_C, EllipseParameters.ellipse_F]
    have h1 : cos (τ + p.δx) * sin p.δy - cos (τ + p.δy) * sin p.δx =
        cos τ * sin (p.δy - p.δx) := by
      rw [cos_add, cos_add, sin_sub]
      ring
    have h2 : cos (τ + p.δx) * cos p.δy - cos (τ + p.δy) * cos p.δx =
        sin τ * sin (p.δy - p.δx) := by
      rw [cos_add, cos_add, sin_sub]
      ring
    trans p.E₀y ^ 2 * (p.E₀x * cos (τ + p.δx)) ^ 2 * (sin p.δy ^ 2 + cos p.δy ^ 2) +
      p.E₀x ^ 2 * (p.E₀y * cos (τ + p.δy)) ^ 2 * (sin p.δx ^ 2 + cos p.δx ^ 2) -
      2 * p.E₀x ^ 2 * p.E₀y ^ 2 * cos (τ + p.δx) * cos (τ + p.δy) * cos (p.δy - p.δx)
      - p.E₀x ^ 2 * p.E₀y ^ 2 * sin (p.δy - p.δx) ^ 2
    · simp only [neg_mul, sin_sq_add_cos_sq, mul_one]
      ring
    trans p.E₀x ^ 2 * p.E₀y ^ 2 * ((cos τ * sin (p.δy - p.δx)) ^ 2 + (sin τ * sin (p.δy - p.δx)) ^ 2
      - sin (p.δy - p.δx) ^ 2)
    · rw [← h1, ← h2]
      rw [cos_sub]
      ring
    trans p.E₀x ^ 2 * p.E₀y ^ 2 * ((cos τ ^ 2 + sin τ ^ 2) * sin (p.δy - p.δx) ^ 2 - sin (p.δy - p.δx) ^ 2)
    · ring
    simp

end EllipseParameters

class IsStandard (E : GeneralEllipse) : Prop where
  hB : E.B = 0

lemma phase_diff (p : EllipseParameters) (hp : p.ellipse.IsStandard) :
    Angle.coe (p.δy - p.δx) = - (π / 2 : ℝ) := by
    -- should the goal be something like (p.E₀x = 0 ∨ p.E₀y = 0) ∨ Angle.coe (p.δy - p.δx) = - (π / 2 : ℝ)?
  have hB := p.ellipse_B
  rw [hp.hB] at hB
  simp at hB
  rcases hB with (hx | hy) | hcos
  · sorry -- E₀x = 0
  · sorry -- E₀y = 0
  · rw [← cos_pi_div_two, ← cos_neg (π/2)] at hcos
    rw [Angle.cos_sin_inj hcos _]
    · simp
    · sorry -- need more info for the sine case?
    /- edit: this might be clearer
    have hcos' : Angle.cos (Angle.coe (p.δy - p.δx)) = 0 := by
      rw [Angle.cos_coe]
      exact hcos
    rw [Angle.cos_eq_zero_iff] at hcos'
    sorry
    -/

noncomputable def standardparametrization (s : EllipseParameters) (hs : s.ellipse.IsStandard) (τ δ₀ : ℝ) : s.ellipse.Point where
  x := s.E₀x * cos (τ + δ₀)
  y := s.E₀y * sin (τ + δ₀)
  -- still not sure how to introduce `δ₀` without explicitly calculating it
  eq := by
    sorry

I think what I want is what I'm calling the standardparametrization, but I'm still not quite sure how to introduce δ₀ without calculating it explicitly. It seems to me that using ∃ δ₀ will still require me to show it.

In addition to this, to prove that the y-coordinate of the parametrization is a sine I'll need to prove phase_diff as you suggested. However the degenerate case is once again getting in the way and I don't know how to justify the phase diff since it doesn't matter in the E₀x = 0 or E₀y = 0 cases.

There are probably more issues down the line, e.g. I'm not sure how to define the rotation on Point directly (i.e. given (p : EllipseParameters), rotate p.ellipse.Point by an angle ψ)? Generally, I'm still not totally convinced that the current implementation is the best, but that could very well be that I'm not familiar enough with the capabilities of defining and accessing the fields of definitions and structures