defining variables for ellipse transformations

Stream: new members

[ZhiKai Pong]

import Mathlib

open Real

/- this doesn't work -/
section

variable {x y τ ω δx δy : ℝ}

variable {ψ : ℝ}

local notation "α" => arctan (y/x)

lemma auxiliary_angle : tan α = y/x := by
  apply Real.tan_arctan

noncomputable def ψ_def := ψ = 2⁻¹ * arctan (tan (2 * α) * cos (δy - δx))

example :
    tan (2 * ψ) = tan (2 * α) * cos (δy - δx) := by
  rw [ψ_def] --failed to rewrite using equation theorems for 'ψ_def'.
  simp

end


/- this doesn't work either -/
section

variable {x y τ ω δx δy : ℝ}

local notation "α" => arctan (y/x)

noncomputable def ψ' := 2⁻¹ * arctan (tan (2 * α) * cos (δy - δx))

example :
    tan (2 * ψ') = tan (2 * arctan (y/x)) * cos (δy - δx) := by --don't know how to synthesize implicit argument
  rw [ψ']
  simp

end


/- the following works -/

noncomputable def ψ (x y δx δy : ℝ) := 2⁻¹ * arctan (tan (2 * arctan (y/x)) * cos (δy - δx))

--is there a way to not carry the `x y δx δy` around
example :
    tan (2 * ψ x y δx δy) = tan (2 * arctan (y/x)) * cos (δy - δx) := by
  rw [ψ]
  simp

I have some expression that depends on x y δx δy, and I would like to represent this expression by a variable ψ. What I actually want is the implicit relation tan (2 * ψ) = tan (2 * α) * cos (δy - δx), but I'm not sure show to state this and if there's a way to do this please let me know.

The code above has some attempts: the third one works, but I don't want to carry around the ψ x y δx δy since those variables are used in all the derivations and definitions in this section. Is there a simple way to let lean understand these without me having to spell everything out? the first two attempts doesn't quite work and I don't really understand why

[ZhiKai Pong]

import Mathlib

open Real

section

variable {x y τ ω δx δy : ℝ}

local notation "α" => arctan (y/x)

local notation  "ψ" => 2⁻¹ * arctan (tan (2 * α) * cos (δy - δx))

example :
    tan (2 * ψ) = tan (2 * α) * cos (δy - δx) := by
  simp

end

I can do this, but this doesn't seem that idiomatic as I want ψ to be an object with meaning rather than just placeholder for an expression.

[Kenny Lau]

@ZhiKai Pong

import Mathlib

open Real

variable (x y δx δy ψ : ℝ) (hψ : ψ = 2⁻¹ * arctan (tan (2 * arctan (y/x)) * cos (δy - δx)))

include hψ in
lemma foo : tan (2 * ψ) = tan (2 * arctan (y/x)) * cos (δy - δx) := by
  rw [hψ]
  simp

[ZhiKai Pong]

Kenny Lau said:

ZhiKai Pong

import Mathlib

open Real

variable (x y δx δy ψ : ℝ) (hψ : ψ = 2⁻¹ * arctan (tan (2 * arctan (y/x)) * cos (δy - δx)))

include hψ in
lemma foo : tan (2 * ψ) = tan (2 * arctan (y/x)) * cos (δy - δx) := by
  rw [hψ]
  simp

thanks, what's the difference between your example and the following? I want hf to be a standalone declaration if possible

import Mathlib

open Real

variable (x y δx δy ψ : ℝ)

noncomputable def hψ := ψ = 2⁻¹ * arctan (tan (2 * arctan (y/x)) * cos (δy - δx))

lemma foo : tan (2 * ψ) = tan (2 * arctan (y/x)) * cos (δy - δx) := by
  rw [hψ] --failed to rewrite using equation theorems for 'hψ'.
  simp

[Kenny Lau]

@ZhiKai Pong why does hψ have to be a def?

[Kenny Lau]

also the difference is that you have defined a Prop called hψ

[Kenny Lau]

you can see that if you do #check hψ

[Kenny Lau]

do you run into any problem if you use the approach I suggested?

[Michael Rothgang]

You can make $\psi$ a definition instead:

import Mathlib

open Real
variable (x y δx δy : ℝ)
noncomputable def ψ : ℝ := 2⁻¹ * arctan (tan (2 * arctan (y/x)) * cos (δy - δx))

lemma foo : tan (2 * ψ x y δx δy) = tan (2 * arctan (y/x)) * cos (δy - δx) := by rw [ψ]; simp

[Michael Rothgang]

Note that $\psi$ depends on x, y, δx and δy.

[Kenny Lau]

@Michael Rothgang that's the same thing as the original message

[ZhiKai Pong]

Kenny Lau said:

ZhiKai Pong why does hψ have to be a def?

Beacuse ψ represents a physical angle see other thread and I want it to have a standalone docstring etc. (Proving the geometry seems out of reach for me at the moment so I've decide to just hardcode it for now)

Kenny Lau said:

also the difference is that you have defined a Prop called hψ

is (hψ : ψ = 2⁻¹ * arctan (tan (2 * arctan (y/x)) * cos (δy - δx))) not a Prop then? my understanding is that a proof of ψ = 2⁻¹ * arctan (tan (2 * arctan (y/x)) * cos (δy - δx))) is an object that has this type, but I'm still not quite clear what actually is a Prop

[ZhiKai Pong]

ok I refreshed my memory with Ch3 of #tpil and I think I can understand that a proof is different from a Prop. so a Prop on its own doesn't carry information, I have to prove this Prop is True for it to be an actual proof that can be used in rw?

[Kenny Lau]

correct

[Kenny Lau]

ZhiKai Pong said:

Kenny Lau said:

ZhiKai Pong why does hψ have to be a def?

Beacuse ψ represents a physical angle see other thread and I want it to have a standalone docstring etc. (Proving the geometry seems out of reach for me at the moment so I've decide to just hardcode it for now)

I know what the context is, my question is about what is the best way to translate a given piece of maths into Lean, which is not always as trivial as you think

[Kenny Lau]

my question is more about the implementation choices

[ZhiKai Pong]

that's also what I'm asking because I don't have a good idea :(

[Kenny Lau]

right, so my point is that, why don't you try my approach first, and see if you bump into any problem later on

[ZhiKai Pong]

yea I think I'll go with your recommendation of using hf for now, and write a short paragraph explaining the implementation for the whole section. Many thanks!

[Kenny Lau]

and then when you have to use it, you can always do _ rfl

[Kenny Lau]

(i.e. let ψ be _, and let hψ be rfl)

[Kenny Lau]

ZhiKai Pong said:

that's also what I'm asking because I don't have a good idea :(

well sometimes for things like this you kinda want to also start at the end and look back

[ZhiKai Pong]

Kenny Lau said:

(i.e. let ψ be _, and let hψ be rfl)

sorry I don't quite get this, do you mind giving an example

[Kenny Lau]

it'd be more helpful actually the other way round, what is your end result?

[Kenny Lau]

It's worth actually giving it more a think of what the goal is and how you want to translate it into Lean

[Kenny Lau]

the literal translation isn't always the best

[ZhiKai Pong]

actually I think all the expressions is just rewriting using foo, so as long as I have foo that I should be good (as in the algebra simplifies to some sin^2 + cos^2 = 1 and everything cancels)

[Kenny Lau]

image.png

[Kenny Lau]

so this is the page you want to translate

[ZhiKai Pong]

right, so the whole thing that I want to do is this:

structure StokesParameters_monochromatic (E₀x E₀y δx δy : ℝ) where
  s0 : ℝ
  s1 : ℝ
  s2 : ℝ
  s3 : ℝ
  hs0 : s0 = E₀x ^ 2 + E₀y ^ 2
  hs1 : s1 = E₀x ^ 2 - E₀y ^ 2
  hs2 : s2 = 2 * E₀x * E₀y * cos (δy - δx)
  hs3 : s3 = 2 * E₀x * E₀y * sin (δy - δx)

theorem stokesparametrizePoincareSphere (s : StokesParameters_monochromatic E₀x E₀y δx δy) :
    s.s1 = s.s0 * cos (2 * χ) * cos (2 * ψ) ∧
    s.s2 = s.s0 * cos (2 * χ) * sin (2 * ψ) ∧
    s.s3 = s.s0 * sin (2 * χ) := by

I want to show that the stokes parameters (I'm not familiar with structure and will update that later) parameterizes a sphere using those angles.

[ZhiKai Pong]

Kenny Lau said:

image.png

In physics the angle is defined from the geometry, but it seems difficult to formalize, so I'm thinking maybe I ignore that for now, hardcode χ and ψ and just verify the algebra.

[Kenny Lau]

@ZhiKai Pong can you explain to me what e0x e0y dx dy s0 s1 s2 s3 are?

[Michael Rothgang]

Here's one way how you could write the above.

import Mathlib

open Real

structure StokesParameters where
  s0 : ℝ
  s1 : ℝ
  s2 : ℝ
  s3 : ℝ

@[simps]
noncomputable def StokesParameters_monochromatic (E₀x E₀y δx δy : ℝ) : StokesParameters where
  s0 := E₀x ^ 2 + E₀y ^ 2
  s1 := E₀x ^ 2 - E₀y ^ 2
  s2 := 2 * E₀x * E₀y * cos (δy - δx)
  s3 := 2 * E₀x * E₀y * sin (δy - δx)

-- need to pass in hypotheses what χ and ψ are
theorem stokesparametrizePoincareSphere (χ ψ : ℝ) (E₀x E₀y δx δy : ℝ) :
    letI s := StokesParameters_monochromatic E₀x E₀y δx δy;
    s.s1 = s.s0 * cos (2 * χ) * cos (2 * ψ) ∧
    s.s2 = s.s0 * cos (2 * χ) * sin (2 * ψ) ∧
    s.s3 = s.s0 * sin (2 * χ) := sorry

[Michael Rothgang]

This doesn't include what $\chi$ and $\psi$ are (as you didn't mention that above).

[Michael Rothgang]

Note the letI inside the main statement: this constructs s as desired, but also tells Lean to replace any mention of the expression s by its definition. (Otherwise, your desired theorem would have type let s := ...; s.s1 = s.s0 * cos (2 * χ) * cos (2 * ψ) ...`.)

[Kenny Lau]

aha, you can do the same with χ and ψ then! (wait, i'm still not sure if that's the best approach)

[Michael Rothgang]

Notably, I'd change the stokes parameters from their particular definition: StokesParameters_monochromatic tells you the values of your parameters, for particular E₀x E₀y δx δy.

[Michael Rothgang]

Kenny Lau said:

aha, you can do the same with χ and ψ then! (wait, i'm still not sure if that's the best approach)

Agreed, with both parts.

[Kenny Lau]

i think i need to know the mathematical significances of everything before i can tell you what the best approach is

[Kenny Lau]

what i'm suggesting is basically, as a middlestep between translating maths to lean, the middlstep is to make the maths more "rigorous" first

[Kenny Lau]

instead of saying "just look at this picture", that's not rigorous mathematics

[Kenny Lau]

you need to figure out mathematically what the relation between all the symbols are first

[Michael Rothgang]

(deleted)

[Michael Rothgang]

Without knowing the maths in your domain too well, the standard way I'd formalise this is the following.

import Mathlib

open Real

noncomputable section
-- TODO: this needs a better name and documentation
structure AngleSomething where
  (x y δx δy : ℝ)

def AngleSomething.χ (A : AngleSomething) := 2⁻¹ * arctan (tan (2 * arctan (A.y / A.x)) * cos (A.δy - A.δx))

def AngleSomething.ψ (A : AngleSomething) : ℝ := sorry -- TODO: fill in a definition

structure StokesParameters where
  s0 : ℝ
  s1 : ℝ
  s2 : ℝ
  s3 : ℝ

@[simps]
def StokesParameters_monochromatic (E₀x E₀y : ℝ) (A : AngleSomething) : StokesParameters where
  s0 := E₀x ^ 2 + E₀y ^ 2
  s1 := E₀x ^ 2 - E₀y ^ 2
  s2 := 2 * E₀x * E₀y * cos (A.δy - A.δx)
  s3 := 2 * E₀x * E₀y * sin (A.δy - A.δx)
-- need to pass in hypotheses what χ and ψ are
theorem stokesparametrizePoincareSphere (A : AngleSomething) (E₀x E₀y : ℝ) :
    letI s := StokesParameters_monochromatic E₀x E₀y A;
    s.s1 = s.s0 * cos (2 * A.χ) * cos (2 * A.ψ) ∧
    s.s2 = s.s0 * cos (2 * A.χ) * sin (2 * A.ψ) ∧
    s.s3 = s.s0 * sin (2 * A.χ) := sorry

[Michael Rothgang]

Explanation

Basically, you can combine all your parameters x, y, \deltax, \deltay into some structure (which I've named AngleSomething --- obviously, you should give it a more descriptive name). Then, you could make \psi a def depending on an instance of AngleSomething. Mathematically, this dependency makes sense: \psi is determined by a choice of these parameters. You can make this fairly short to write (A.\psi in my example) using "dot notation".
I don't know if including E₀x E₀y in your AngleSomething structure makes sense, or you'd e.g. prefer a structure containing an AngleSomething and these two parameters, or something else. As I said, don't know very well what you're doing mathematically.

[ZhiKai Pong]

ok let me attempt to tell the full story:

consider a time-harmonic electromagnetic planewave propagating in the z-direction. since EM waves are transverse, the general form is

local notation "Ex" => E₀x * cos (τ + δx)
local notation "Ey" => E₀y * cos (τ + δy)

the locus of these parametric equations give an ellipse known as the polarization ellipse:

theorem polarizationEllipse (hx : E₀x ≠ 0) (hy : E₀y ≠ 0) :
    (Ex / E₀x)^2 + (Ey / E₀y)^2 - 2 * (Ex / E₀x) * (Ey / E₀y) * cos (δy - δx) =
    sin (δy - δx) ^ 2

here's the step that I tried to ask about yesterday, but I've come to the conclusion that it's too handwavy and difficult to justify:
image.png

I think going from (18) to (19) requires knowledge of conic sections etc. (I don't even know that an ellipse after rotation is still an ellipse so that's a massive jump in logic)

The gist is that consider rotating the ellipse to standard orientation, then after suitable simplification the expressions χ and ψ can be obtained
image.png
where tan α = Ey/Ex and δ = δy - δx.

If we define the stokes parameters as

  s0 := E₀x ^ 2 + E₀y ^ 2
  s1 := E₀x ^ 2 - E₀y ^ 2
  s2 := 2 * E₀x * E₀y * cos (A.δy - A.δx)
  s3 := 2 * E₀x * E₀y * sin (A.δy - A.δx)

then we have the relationships

    s.s1 = s.s0 * cos (2 * A.χ) * cos (2 * A.ψ) ∧
    s.s2 = s.s0 * cos (2 * A.χ) * sin (2 * A.ψ) ∧
    s.s3 = s.s0 * sin (2 * A.χ)

full code

import Mathlib

open Real

variable {E₀x E₀y τ δx δy : ℝ}

local notation "Ex" => E₀x * cos (τ + δx)
local notation "Ey" => E₀y * cos (τ + δy)

/-- The locus of the electric field vector of an monochromatic time-harmonic
  electromagnetic plane wave is an ellipse. -/
theorem polarizationEllipse (hx : E₀x ≠ 0) (hy : E₀y ≠ 0) :
    (Ex / E₀x)^2 + (Ey / E₀y)^2 - 2 * (Ex / E₀x) * (Ey / E₀y) * cos (δy - δx) =
    sin (δy - δx) ^ 2 := by
  have h1 : (E₀x * cos (τ + δx)) / E₀x * sin δy - (E₀y * cos (τ + δy)) / E₀y * sin δx =
      cos τ * sin (δy - δx) := by
    field_simp
    rw [cos_add, cos_add, sin_sub]
    ring
  have h2 : (E₀x * cos (τ + δx)) / E₀x * cos δy - (E₀y * cos (τ + δy)) / E₀y * cos δx =
      sin τ * sin (δy - δx) := by
    field_simp
    rw [cos_add, cos_add, sin_sub]
    ring
  trans (E₀x * cos (τ + δx) / E₀x) ^ 2 * (sin δy ^ 2 + cos δy ^ 2) +
      (E₀y * cos (τ + δy) / E₀y) ^ 2 * (sin δx ^ 2 + cos δx ^ 2) -
      2 * (E₀x * cos (τ + δx) / E₀x) * (E₀y * cos (τ + δy) / E₀y) * 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

noncomputable section

structure monochrom_wave_parameters where
  (E₀x E₀y δx δy : ℝ)

def monochrom_wave_parameters.χ (A : monochrom_wave_parameters) := 2⁻¹ * arctan (tan (2 * arctan (A.E₀y / A.E₀x)) * cos (A.δy - A.δx))

def monochrom_wave_parameters.ψ (A : monochrom_wave_parameters) := 2⁻¹ * arcsin (sin (2 * arctan (A.E₀y / A.E₀x)) * sin (A.δy - A.δx))

structure StokesParameters where
  s0 : ℝ
  s1 : ℝ
  s2 : ℝ
  s3 : ℝ

@[simps]
def StokesParameters_monochromatic (E₀x E₀y : ℝ) (A : monochrom_wave_parameters) : StokesParameters where
  s0 := E₀x ^ 2 + E₀y ^ 2
  s1 := E₀x ^ 2 - E₀y ^ 2
  s2 := 2 * E₀x * E₀y * cos (A.δy - A.δx)
  s3 := 2 * E₀x * E₀y * sin (A.δy - A.δx)

-- need to pass in hypotheses what χ and ψ are
theorem stokesparametrizePoincareSphere (A : monochrom_wave_parameters) (E₀x E₀y : ℝ) :
    letI s := StokesParameters_monochromatic E₀x E₀y A;
    s.s1 = s.s0 * cos (2 * A.χ) * cos (2 * A.ψ) ∧
    s.s2 = s.s0 * cos (2 * A.χ) * sin (2 * A.ψ) ∧
    s.s3 = s.s0 * sin (2 * A.χ) := by
  sorry

[ZhiKai Pong]

so what's troubling me is that ideally I would want to derive the expressions for χ and ψfrom the geometry, but it seems to complicated for me at the moment. If I hard code them as I'm asking in this thread, I'll still be able to show that they parameterize the sphere (just that we don't know where those expressions come from). Physically this result says that every possible "state of polarization" (all possible forms of ellipses) can be represented by a unique point on the sphere.

[Kenny Lau]

I agree that we don't want the geometry, but what exactly is the process of de-geometrisation?

[Kenny Lau]

ZhiKai Pong said:

I think going from (18) to (19) requires knowledge of conic sections etc. (I don't even know that an ellipse after rotation is still an ellipse so that's a massive jump in logic)

no, I believe this is all just some explicit trigonometry

[ZhiKai Pong]

I agree that everything that comes afterwards is just trigonometry.

image.png
image.png

It's precisely the step going form (18) to (19) I don't know how to state in lean.

[ZhiKai Pong]

Eq. (13) is this

Ex = E₀x * cos (τ + δx)
Ey = E₀y * cos (τ + δy)

[ZhiKai Pong]

so equation (19) is claiming that we can "just write down the parametric from of a ellipse in standard orientation" which I believe is not justified, together with the rotation in (18) which is also geometrical

[Michael Rothgang]

Kenny Lau said:

ZhiKai Pong said:

I think going from (18) to (19) requires knowledge of conic sections etc. (I don't even know that an ellipse after rotation is still an ellipse so that's a massive jump in logic)

no, I believe this is all just some explicit trigonometry

To elaborate: you can write down a rotation as an explicit matrix (in terms of sine and cosine), so this particular fact boils down to some trigonometry computations.

[Kenny Lau]

@ZhiKai Pong what are a1 and a2?

[ZhiKai Pong]

Michael Rothgang said:

To elaborate: you can write down a rotation as an explicit matrix (in terms of sine and cosine), so this particular fact boils down to some trigonometry computations.

I understand this, but saying that the parameterization step and the rotation commute(?) is something I don't know how to make rigorous.

[ZhiKai Pong]

Kenny Lau said:

ZhiKai Pong what are a1 and a2?

sorry for the confusion of notation here, a1 a2 are just E₀x and E₀y, which are the dimensions of the original ellipse, and a and b are the dimensions of the rotated ellipse (major and minor axes).

[ZhiKai Pong]

from another textbook
image.png

[Kenny Lau]

@ZhiKai Pong consider the following translation:

import Mathlib

noncomputable section

open Real

-- Ax^2+By^2+Cxy+Dx+Ey=F
structure Ellipse : Type where
  A : ℝ
  B : ℝ
  C : ℝ
  D : ℝ
  E : ℝ
  F : ℝ
  -- maybe insert some positivity condition

namespace Ellipse

/-- A point on the ellipse -/
structure Point (E : Ellipse) : Type where
  x : ℝ
  y : ℝ
  hxy : E.A * x^2 + E.B * y^2 + E.C * x * y + E.D * x + E.E * y = 1

/-- The information needed to parametrise an ellipse, namely a, b, δ₀ -/
structure Parameters : Type where
  a : ℝ
  ha : 0 < a
  b : ℝ
  hb : 0 < b
  δ₀ : ℝ

namespace Parameters

/-- The ellipse parametrised by a set of parameters. -/
def ellipse (p : Parameters) : Ellipse :=
  sorry

/-- The point parametrised -/
def eval (p : Parameters) (τ : ℝ) : p.ellipse.Point where
  x := p.a * cos (τ + p.δ₀)
  y := p.b * cos (τ + p.δ₀)
  hxy := sorry

def χ (p : Parameters) : ℝ := arctan (p.b / p.a)

end Parameters

class IsStandard (E : Ellipse) : Prop where
  ha : 0 < E.A
  hb : 0 < E.B
  hc : E.C = 0
  hd : E.D = 0
  he : E.E = 0
  hf : 0 < E.F

/-- Rotate the ellipse parametrised by (a, b, 0) by angle ψ. -/
def rotate (a b ψ : ℝ) : Parameters :=
  sorry

-- I don't know what these are
def a₁ : ℝ := sorry
def a₂ : ℝ := sorry
def α : ℝ := arctan (a₂ / a₁)
def δ : ℝ := sorry

-- (31)
theorem main_result (a b ψ : ℝ) (P₁) (hp₁ : P₁ = rotate a b ψ) :
    a ^ 2 + b ^ 2 = a₁ ^ 2 + a₂ ^ 2 ∧
    tan (2 * ψ) = tan (2 * α) * cos δ ∧
    sin (2 * P₁.χ) = sin (2 * α) * sin δ :=
  sorry

end Ellipse

[Kenny Lau]

and from this I hope you understand me more when I say "you really need to think about how to translate a given piece of maths into Lean, and it is not trivial"

[Kenny Lau]

like, you really need to think about what everything means, and what are the relations between the different symbols, as I tried to say above

[Kenny Lau]

the first question you should ask yourself is, what is an ellipse, and what do you mean by parametrising an ellipse?

[Kenny Lau]

and for this question I hope you will agree that Ellipse and Ellipse.Parameters above are the correct translations

[ZhiKai Pong]

@Kenny Lau many thanks for this, I think this is a good starting point for me to work on!

This sort of exposition of definitions is what I'm currently trying to grasp, appreciate it a lot :)

[ZhiKai Pong]

I have a follow-up question (sorry didn't find the time to work on this until now):

full code

import Mathlib

open Real

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

namespace GeneralEllipse

structure EllipseParameters where
  E₀x : ℝ
  E₀y : ℝ
  δx : ℝ
  δy : ℝ
  non_zero : ¬((E₀y = 0 ∧ E₀x = 0))

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

@[simps]
noncomputable def EllipseParameters.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
    by_cases h : sin (p.δy - p.δx) = 0
    · right
      simp [h]
    · have hxy : ¬((p.E₀y = 0 ∨ p.E₀x = 0)) := by
        have h_non_zero := p.non_zero
        push_neg
        constructor
        · by_contra
          sorry -- if p.E₀x = 0, E.A * x ^ 2 + E.B * x * y + E.C * y ^ 2 + E.F = 0 implies E.A * x ^ 2 = 0, therefore p.E₀y = 0
        · by_contra
          sorry
      left
      push_neg at h
      calc
        _ = - 4 * p.E₀y ^ 2 * p.E₀x ^ 2 * (1 - cos (p.δy - p.δx) ^ 2) := by ring
        _ = - (2 * p.E₀y * p.E₀x * sin (p.δy - p.δx)) ^ 2 := by rw [← sin_sq]; ring
        _ < 0 := by
          apply neg_pos.mp
          simp only [neg_neg, sq_pos_iff, ne_eq, mul_eq_zero, OfNat.ofNat_ne_zero, false_or, hxy]
          exact h

noncomputable def EllipseParameters.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

The ellipse I'm defining wants to take into account the cases where it degenerates into a straight line, therefore I have the condition h : B ^ 2 - 4 * A * C < 0 ∨ F = 0. While proving the coefficients defined by

  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

satisfies the condition, I need to prove that sin (p.δy - p.δx) = 0 implies ¬((p.E₀y = 0 ∨ p.E₀x = 0)). However to prove this I think I'll need to use the equation form E.A * x ^ 2 + E.B * x * y + E.C * y ^ 2 + E.F = 0, but I'm not sure how to change the definition to make this work while keeping the structure you suggested (which is useful beacuse I do want to show EllipseParameters.eval give the correct form of parameterization in the next def). Any suggestions on how do I modify this to fill in the missing sorrys?

[Kenny Lau]

@ZhiKai Pong well, good to see that you're adapting well to my suggestions. I had a look at your code, and it seems like you artificially introduced those goals there.

The maths proof goes like, if sin (p.δy - p.δx) = 0, F is 0, so you're good; otherwise, then cos (p.δy - p.δx) is not +-1, so the discriminant is negative!

[ZhiKai Pong]

Maybe I wrote this too hastily, please allow me to clarify: currently the logic goes like "an ellipse is defined by 4 real numbers A B C F satisfying some condition" but makes no reference to the actual algebraic equation, which is then defined in GeneralEllipse.Point. This feels a bit strange to me, but I can't quite pinpoint how to restructure x y and A B C F to allow me to use the equation.

[Kenny Lau]

@ZhiKai Pong I'm saying that Line 36 is your artificially introduced goal, which is not mathematically sound

[Kenny Lau]

your definitions before that are mathematically correct

[Kenny Lau]

in other words, all your definitions are fine, all your theorems are valid, but you took the wrong turn in a proof, that's why you're stuck

[ZhiKai Pong]

Kenny Lau said:

ZhiKai Pong I'm saying that Line 36 is your artificially introduced goal, which is not mathematically sound

I got that goal while working through the calc block, which makes sense to me because the final step of the goal -(2 * p.E₀y * p.E₀x * sin (p.δy - p.δx)) ^ 2 < 0 requires me to prove neither E₀x nor E₀y zero.

[ZhiKai Pong]

I think usually when we talk about ellipses we require the major and minor axis to both be non-zero so this is avoided, but now I'm allowing the degenerate cases where one of them can be zero, just not both at the same time.

[ZhiKai Pong]

Kenny Lau said:

The maths proof goes like, if sin (p.δy - p.δx) = 0, F is 0, so you're good; otherwise, then cos (p.δy - p.δx) is not +-1, so the discriminant is negative!

maybe I'm misunderstanding this part - do you mind elaborating how does cos not equal to +-1 imply the discriminant is negative?

[ZhiKai Pong]

I think it's true if I just had B ^ 2 - 4 * A * C ≤ 0, but I think the definition I had is the correct one in the sense that B ^ 2 - 4 * A * C ≤ 0 allows for parabolas

[Kenny Lau]

@ZhiKai Pong have you tried writing out b^2-4ac by hand?

[ZhiKai Pong]

the calc block is mimicking what I'll do by hand

[Kenny Lau]

ok, fair enough, I missed that, my apologies

[Kenny Lau]

so as you can see, you ended up with - (2 * p.E₀y * p.E₀x * sin (p.δy - p.δx)) ^ 2

[Kenny Lau]

can you tell me mathematically when this is negative?

[ZhiKai Pong]

when p.E₀y p.E₀x sin (p.δy - p.δx) are all non-zero. The last one I have from the assumption, the other two is the reason I had the introduced goal hxy.

[Kenny Lau]

wait, i just realised

[Kenny Lau]

isn't b^2-4ac just 4 times f?

[Kenny Lau]

so h is just equivalent to f <= 0

[Kenny Lau]

that simplifies your proof

[Kenny Lau]

so, start over from L33, just case on the whole thing being 0

[ZhiKai Pong]

aha this works!! (any advice on how to write this in a better way?)

noncomputable def EllipseParameters.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
    by_cases h : p.E₀x ^ 2 * p.E₀y ^ 2 * sin (p.δy - p.δx) ^ 2 = 0
    · right
      simp [h]
    · left
      push_neg at h
      calc
        _ = - 4 * p.E₀y ^ 2 * p.E₀x ^ 2 * (1 - cos (p.δy - p.δx) ^ 2) := by ring
        _ = - (2 * p.E₀y * p.E₀x * sin (p.δy - p.δx)) ^ 2 := by rw [← sin_sq]; ring
        _ < 0 := by
          apply neg_pos.mp
          repeat rw [mul_ne_zero_iff] at h
          repeat rw [pow_two] at h
          repeat rw [mul_self_ne_zero] at h
          simp only [neg_neg, sq_pos_iff, ne_eq, mul_eq_zero, OfNat.ofNat_ne_zero, false_or, not_or]
          conv =>
            enter [1]
            rw [And.comm]
          exact h

Thanks a lot!

[Kenny Lau]

Congratulations!

[Kenny Lau]

@ZhiKai Pong

@[simps]
noncomputable def EllipseParameters.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

[ZhiKai Pong]

this is very nice, thanks again!

[ZhiKai Pong]

Unfortunately I'm stuck on the next step - I've attempted to define the rotation, but it seems like I'm missing something.

full code

import Mathlib

open Real

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

namespace GeneralEllipse

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

noncomputable def χ (p : EllipseParameters) : ℝ := Complex.arg (p.E₀x + Complex.I * p.E₀y)

@[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

/-- Rotate the ellipse parametrised by `p` by angle ψ to the standard orientation. -/
/- this probably needs modification but not sure how -/
@[simps]
noncomputable def ellipse_rotate (p : EllipseParameters) (δ₀ ψ : ℝ) : EllipseParameters where
  E₀x := p.E₀x * cos ψ  + p.E₀y * sin ψ
  E₀y := - p.E₀x * sin ψ + p.E₀y * cos ψ
  δx := δ₀
  δy := δ₀ + π/2

lemma rotation_is_standard {p s : EllipseParameters} (δ₀ ψ : ℝ)
    (hs : s = ellipse_rotate p δ₀ ψ) :
    IsStandard s.ellipse where
  hB := by simp [hs]

/- Attempt at stating this : -/
lemma evalx_eq {p : EllipseParameters} {δ₀ ψ τ : ℝ} :
    ((ellipse_rotate p δ₀ ψ).eval τ).x = (p.eval τ).x * cos ψ + (p.eval τ).y * sin ψ := by
  sorry

/- I want this result: -/
lemma s_eval_Ex {p s : EllipseParameters} (δ₀ ψ τ : ℝ) (hs : s = ellipse_rotate p δ₀ ψ) :
    s.E₀x * cos (τ + δ₀) = p.E₀x * cos (τ + p.δx) * cos ψ + p.E₀y * cos (τ + p.δy) * sin ψ := by
  trans (s.eval τ).x
  · simp [hs]
  rw [hs, evalx_eq]
  rfl

I've defined ellipse_rotation by stating how the ellipse parameters transform. However, there is another rotation here, which is how the Point actually transform. I want the result at the end, which follows from the two rotations coinciding, and evalx_eq is an attempt at stating this result but I'm not quite sure how to state/prove this properly. I've also tried defining the ellipse_rotation directly on the Point instead, but that didn't get me anywhere either :/ If you have time any help is much appreciated

I think on paper I'll define

  δx := δ₀
  δy := δ₀ + π/2

to be a property of IsStandard, but I'm not sure where everything should go or if current definitions need some tweaking.

[Kenny Lau]

@ZhiKai Pong since you're doing away with D and E, i.e. your ellipse is centered at the origin, then rotation (of Point) should just be rotation.

Just a note, I feel like you're trying to replace doing maths on pen-and-paper by Lean, which might not be the best approach; you should definitely first verify on paper (or using wolframalpha, mathematica, pari/gp, desmos, etc.) all the formulas you want to verify, before typing it out in Lean

Let's start over mathematically and make sure everything is neat. Tell me your result to each question as a function of t, depending on additional constants (i.e. parameters)

1. What is the standard parametrisation of the standard unit circle? Make sure that when t=0, the result you get is (1,0), and that it goes counter-clockwise and comes back when t=2π. (You can use Desmos to confirm.)
2. Now can you change the parametrisation to make it into a standard ellipse? How many extra parameters do you need at this stage? Still, make sure that you start at the positive x-axis when t=0.
3. What is the equation of the ellipse you formed by your previous parametrisation in 2.?
4. Now can you rotate your parametrised point counter-clockwise by a degree of ψ?
5. Now can you introduce a phase difference δ0 so that you don't start on an axis (of the ellipse)?

See, at this point I run the risk of deviating from your source paper. I apologise if this is the case. I tried to reread the screenshots you sent me, but I got lost because your screenshots did not contain a definition of δ1 and δ2. Can you try to read the paper again to find those definitions? Then I might have to change the questions above.

[ZhiKai Pong]

Well from my perspective I'd like to think I understand the equations quite well, it's just that I don't know how to translate it into lean. Maybe what you're observing is that there are gaps in my thinking that I usually wouldn't realize if I'm working on paper and I'm trying to flesh them out while working in lean, but I'm not quite used to this yet (as you can probably tell I'm not a mathematician).

Anyway here's my attempt at your questions:

1. standard parametrization of unit circle

x := cos (t)
y := sin (t) -- or cos(t - π/2)

2.a Parametrization of standard ellipse

x := a * cos (t)
y := b * sin (t) -- or b * cos (t - π/2)

3.a Equation of standard ellipse in (2.a)

b ^ 2 * x ^ 2 + a ^ 2 * y ^ 2 - a ^ 2 * b ^ 2 = 0
-- more idiomatic form would be
(x/a)^2 + (y/b)^2 = 1
-- but I want to allow for the cases where a or b is 0

4.a Rotation of parametrised point counter-clockwise by ψ

-- multiply expression in (2.a) by rotation matrix
x' := x * cos ψ - y * sin ψ = a * cos t * cos ψ - b * sin t * sin ψ
y' := x * sin ψ + y * cos ψ = a * cos t * sin ψ + b * cos t * cos ψ

5.a Introducing a phase difference δ₀ for a standard ellipse

-- the expression in (4.a) should be changed to
x := a * cos (t + δ₀)
y := b * sin (t + δ₀) -- or cos(t + δ₀- π/2)

But this is only half of the story, to follow the textbook there is an alternative formulation

2.b Parametrization of general ellipse

x' := a₁ * cos (t)
y' := b₁ * cos (t + δ) -- where δ = δ₂ - δ₁

3.b Equation of general ellipse in (2.b)

b₁ ^ 2 * x' ^ 2 - 2 * a₁ * b₁ * cos δ * x' * y' + a₁ ^ 2 * y' ^ 2 - a₁ ^ 2 * b₁ ^ 2 * sin δ ^ 2 = 0

5.b Introducing a phase difference δ₁ for a general ellipse

x' := a₁ * cos (t + δ₁)
y' := b₁ * cos (t + δ₂) -- δ₂ = δ + δ₁

4.b Rotating point of general ellipse clockwise by angle ψ

x := x' * cos ψ + y' * sin ψ = a * cos (t + δ₁) * cos ψ - b * cos (t + δ₂) * sin ψ
y := - x' * sin ψ + y' * cos ψ = a * cos (t + δ₁) * sin ψ + b * cos (t + δ₂) * cos ψ

6. Rotate the general ellipse to standard orientation (how to state this in lean?), then equate 4.b and 5.a

x = a * cos (t + δ₀) = x' * cos ψ + y' * sin ψ = a * cos (t + δ₁) * cos ψ - b * cos (t + δ₂) * sin ψ
y = b * sin (t + δ₀) = - x' * sin ψ + y' * cos ψ = a * cos (t + δ₁) * sin ψ + b * cos (t + δ₂) * cos ψ
-- which is the expression in the textbook

Apologies for the confusion, now the notation should follow the textbook more closely.

So to me the difficulty is I don't know how to state the formulation via (a) is discussing the same thing as the formulation via (b) and making the equality in 6.

On a more technical level, I agree the rotation should be defined on Point (which is consistent with the exposition above), but I don't know how to access the fields since currently Point is defined from the EllipseParameters. I don't have a clear picture of hierachies in mind on which comes first (parameters, ellipse i.e. equation coefficients, or the parametrization?) and I think that's what is getting me stuck.

(I'm aware this is long and appreciate your time in helping out!)

[ZhiKai Pong]

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

[Notification Bot]

3 messages were moved from this topic to #new members > Check consistency between two statements by Kenny Lau.

[ZhiKai Pong]

(or if you have the permissions would you mind editing the title for this thread? something like defining variables for ellipse transformations, thanks)

[Kenny Lau]

@Kevin Buzzard I don't have permission, could you please help me move this thread?

[Kenny Lau]

@ZhiKai Pong Yes, you should add the degenerate conditions to phase_diff. Note that the difference between the angles can also be pi/2 -- in which case the orientation is flipped, but the ellipse itself is still in the standard orientation (axis parallel to axis).

The way I would define standard parametrisation is through a Prop:

structure EllipseParameters.IsStandard (p : EllipseParameters) : Prop where
  standard : p.E₀x = 0 ∨ p.E₀y = 0 ∨ cos (p.δy - p.δx) = 0

Note that since this is a Prop defined on EllipseParameters, there is then no need to produce any explicit or implicit δ₀. (But in any case, δ₀ is just p.δx, by comparing the equation for x.)

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

lemma isStandard_def (E : GeneralEllipse) :
    E.IsStandard ↔ E.B = 0 :=
  ⟨fun ⟨h⟩ ↦ h, fun h ↦ ⟨h⟩⟩

class EllipseParameters.IsStandard (p : EllipseParameters) : Prop where
  standard : p.E₀x = 0 ∨ p.E₀y = 0 ∨ cos (p.δy - p.δx) = 0

lemma EllipseParameters.isStandard_def (p : EllipseParameters) :
    p.IsStandard ↔ p.E₀x = 0 ∨ p.E₀y = 0 ∨ cos (p.δy - p.δx) = 0 :=
  ⟨fun ⟨h⟩ ↦ h, fun h ↦ ⟨h⟩⟩

lemma EllipseParameters.isStandard_iff (p : EllipseParameters) :
    p.ellipse.IsStandard ↔ p.IsStandard := by
  rw [isStandard_def, GeneralEllipse.isStandard_def, p.ellipse_B]
  simp [or_assoc]

With these correct definitions you'll notice that each proof comes for free :upside_down:

[Kenny Lau]

Note that the connection between cosine being 0 and the value of the angle is for free at  Real.Angle.cos_eq_zero_iff

[Kenny Lau]

If you want a challenge, you can try to prove the following theorem:

lemma EllipseParameters.isStandard_iff' (p : EllipseParameters) :
    p.IsStandard ↔ ∃ E₀x E₀y δ₀ : ℝ, ∀ t,
      (p.eval t).x = E₀x * cos (t + δ₀) ∧ (p.eval t).y = E₀y * sin (t + δ₀) :=
  sorry

[Kenny Lau]

this might be the connection you wanted

[Kenny Lau]

(you'll even discover that the theorem is false! finding what's wrong is left to the reader as an exercise :melt: )

[Kenny Lau]

ZhiKai Pong said:

I'm not sure how to define the rotation on Point directly (i.e. given (p : EllipseParameters), rotate p.ellipse.Point by an angle ψ)?

For this, you define rotation on the whole R^2, and also define how to form the correct new parameters, and show that "everything works out".

[Kenny Lau]

def rotateR2 (p : ℝ × ℝ) (ψ : ℝ) : ℝ × ℝ :=
  (sorry, sorry)

def Point.toR2 {E : GeneralEllipse} (P : E.Point) : ℝ × ℝ :=
  (P.x, P.y)

def EllipseParameters.rotate (p : EllipseParameters) (ψ : ℝ) : EllipseParameters where
  E₀x := sorry
  E₀y := sorry
  δx := sorry
  δy := sorry

-- This might be harder. What's the equation of an ellipse after rotation?
def rotate (E : GeneralEllipse) (ψ : ℝ) : GeneralEllipse where
  A := sorry
  B := sorry
  C := sorry
  F := E.F -- F shouldn't be changed.
  h := sorry

def Point.rotate {E : GeneralEllipse} (P : E.Point) (ψ : ℝ) : (E.rotate ψ).Point where
  x := (rotateR2 P.toR2 ψ).1
  y := (rotateR2 P.toR2 ψ).2
  eq := sorry

-- First correctness theorem:
-- The new ellipse formed by the new parameters, is the old ellipse rotated by ψ.
theorem everything_works_out? (p : EllipseParameters) (ψ : ℝ) :
    p.ellipse.rotate ψ = (p.rotate ψ).ellipse :=
  _

-- Second correctness theorem:
-- Even the parametrised point (at t=t) is correctly cotated.
-- Unfortunately we're getting into dependent type hell here, I'll figure something out afterwards...
theorem everything_works_out! (p : EllipseParameters) (ψ t : ℝ) :
    (p.rotate ψ).eval t = cast (congrArg Point <| everything_works_out? p ψ) ((p.eval t).rotate ψ) :=
  sorry

[Kenny Lau]

@ZhiKai Pong I've provided the skeleton for you. Try to understand my approach (and critique/ask if you want to!) and then you can try to have fun with the sorrys.

[ZhiKai Pong]

thanks so much again, I'll give it a go :)

[ZhiKai Pong]

long code

import Mathlib

open Real

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

namespace GeneralEllipse

structure EllipseParameters where
  a : ℝ
  b : ℝ
  δ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

noncomputable def χ (p : EllipseParameters) : ℝ := Complex.arg (p.a + Complex.I * p.b)

@[simps]
noncomputable def ellipse (p : EllipseParameters) : GeneralEllipse where
  A := p.b ^ 2
  B := - 2 * p.a * p.b * cos (p.δy - p.δx)
  C := p.a ^ 2
  F := - p.a ^ 2 * p.b ^ 2 * sin (p.δy - p.δx) ^ 2
  h := by
    have : 0 ≤ (p.a * p.b * sin (p.δy - p.δx)) ^ 2 := sq_nonneg _
    obtain h | h := eq_or_lt_of_le this
    · right
      calc
        _ = - (p.a * p.b * sin (p.δy - p.δx)) ^ 2 := by ring
        _ = 0 := by rw [← h]; ring
    · left
      calc
        _ = - 4 * p.b ^ 2 * p.a ^ 2 * (1 - cos (p.δy - p.δx) ^ 2) := by ring
        _ = - 4 * (p.a * p.b * 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.a * cos (τ + p.δx)
  y := p.b * 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.b ^ 2 * (p.a * cos (τ + p.δx)) ^ 2 * (sin p.δy ^ 2 + cos p.δy ^ 2) +
      p.a ^ 2 * (p.b * cos (τ + p.δy)) ^ 2 * (sin p.δx ^ 2 + cos p.δx ^ 2) -
      2 * p.a ^ 2 * p.b ^ 2 * cos (τ + p.δx) * cos (τ + p.δy) * cos (p.δy - p.δx)
      - p.a ^ 2 * p.b ^ 2 * sin (p.δy - p.δx) ^ 2
    · simp only [neg_mul, sin_sq_add_cos_sq, mul_one]
      ring
    trans p.a ^ 2 * p.b ^ 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.a ^ 2 * p.b ^ 2 * ((cos τ ^ 2 + sin τ ^ 2) * sin (p.δy - p.δx) ^ 2 -
      sin (p.δy - p.δx) ^ 2)
    · ring
    simp

structure IsStandard (p : EllipseParameters) : Prop where
  standard : p.a = 0 ∨ p.b = 0 ∨ Angle.cos (↑(p.δy - p.δx)) = 0

lemma isStandard_def (p : EllipseParameters) :
    p.IsStandard ↔ p.a = 0 ∨ p.b = 0 ∨ Angle.cos (↑(p.δy - p.δx)) = 0 :=
  ⟨fun ⟨h⟩ ↦ h, fun h ↦ ⟨h⟩⟩

end EllipseParameters

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

lemma isStandard_def (E : GeneralEllipse) :
    E.IsStandard ↔ E.B = 0 :=
  ⟨fun ⟨h⟩ ↦ h, fun h ↦ ⟨h⟩⟩

lemma EllipseParameters.isStandard_iff (p : EllipseParameters) :
    p.ellipse.IsStandard ↔ p.IsStandard := by
  rw [isStandard_def, GeneralEllipse.isStandard_def, p.ellipse_B, Angle.cos_coe]
  simp [or_assoc]

lemma EllipseParameters.isStandard_parametrization (p : EllipseParameters) :
    p.IsStandard → ∃ a b δ₀ : ℝ, ∀ t,
    (p.eval t).x = a * cos (t + δ₀) ∧
    ((p.eval t).y = b * sin (t + δ₀) ∨ (p.eval t).y = - b * sin (t + δ₀)) := by
  intro h
  simp
  rw [isStandard_def] at h
  rcases h with hx | hy | hδ
  · use 0, p.b, p.δy + π/2
    intro t
    simp
    constructor
    · left
      exact hx
    · left
      simp [← add_assoc, sin_add_pi_div_two]
  · use p.a, 0, p.δx
    intro t
    simp
    left
    exact hy
  · rw [Angle.cos_eq_zero_iff] at hδ
    use p.a, p.b, p.δx
    intro t
    rcases hδ with hpos | hneg
    · simp
      have h : (p.δy : Angle) = ↑p.δx + ↑(π / 2) := by
        rw [← sub_left_inj (a := ↑p.δx), ← Angle.coe_sub, hpos]
        simp
      right
      simp [← Angle.cos_coe, ← Angle.sin_coe, ← add_assoc, h, Angle.cos_add_pi_div_two]
    · rw [neg_div] at hneg
      simp
      have h : (p.δy : Angle) = ↑p.δx - ↑(π / 2) := by
        rw [← sub_left_inj (a := ↑p.δx), ← Angle.coe_sub, hneg]
        simp
      left
      left
      simp [← Angle.cos_coe, ← Angle.sin_coe, ← add_sub_assoc, h, Angle.cos_sub_pi_div_two]

@[simps]
noncomputable def rotate (E : GeneralEllipse) (ψ : ℝ) : GeneralEllipse where
  A := E.A * cos ψ ^ 2 + E.B * sin ψ * cos ψ + E.C * sin ψ ^ 2
  B := 2 * (E.C - E.A) * sin ψ * cos ψ + E.B * (cos ψ ^ 2 - sin ψ ^ 2)
  C := E.A * sin ψ ^ 2 - E.B * sin ψ * cos ψ + E.C * cos ψ ^ 2
  F := E.F
  h := by
    obtain h1 | h2 := E.h
    · left
      ring_nf
      calc
        _ = E.B ^ 2 * (sin ψ ^ 2 + cos ψ ^ 2) ^ 2 -
            4 * E.A * E.C * (sin ψ ^ 2 + cos ψ ^ 2) ^ 2 := by ring
        _ = E.B ^ 2 - 4 * E.A * E.C := by simp
        _ < 0 := by exact h1
    · right
      exact h2

lemma rotate_to_standard (p : EllipseParameters) (ψ : ℝ)
    (h : IsStandard (p.ellipse.rotate ψ)) :
    (p.a ^ 2 - p.b ^ 2) * sin (2 * ψ) =
    2 * p.a * p.b * cos (p.δy - p.δx) * cos (2 * ψ) := by
  rw [isStandard_def] at h
  simp at h
  rw [sin_two_mul, cos_two_mul']
  linarith

-- I think everything up to this point is correct

-- Note: require explicit calculation of rotated δx?
noncomputable def EllipseParameters.rotate (p : EllipseParameters) (ψ : ℝ) :
    EllipseParameters where
  a := p.a * cos ψ - p.b * sin ψ
  b := p.a * sin ψ + p.b * cos ψ
  δx := sorry
  δy := sorry

@[simp]
noncomputable def rotateR2 (p : ℝ × ℝ) (ψ : ℝ) : ℝ × ℝ :=
  (p.1 * cos ψ + p.2 * sin ψ, - p.1 * sin ψ + p.2 * cos ψ)

@[simp]
def Point.toR2 {E : GeneralEllipse} (P : E.Point) : ℝ × ℝ :=
  (P.x, P.y)

noncomputable def Point.rotate {E : GeneralEllipse} (P : E.Point) (ψ : ℝ) :
    (E.rotate ψ).Point where
  x := (rotateR2 P.toR2 ψ).1
  y := (rotateR2 P.toR2 ψ).2
  eq := by
    simp
    calc
      _ = (P.x ^ 2 * E.A + P.x * P.y * E.B + P.y ^ 2 * E.C) *
          (sin ψ ^ 2 + cos ψ ^ 2) ^ 2 + E.F := by ring
      _ = (P.x ^ 2 * E.A + P.x * P.y * E.B + P.y ^ 2 * E.C) + E.F := by simp
      _ = 0 := by rw [← P.eq]; ring

-- I think I want this, can this be defined purely through Point.rotate?
lemma parametrization_eq (p s : EllipseParameters) (ψ : ℝ)
    (hs : IsStandard s.ellipse) (hr : s.ellipse = p.ellipse.rotate ψ) :
    (s.eval t).x = ((p.eval t).rotate ψ).x := by
  sorry

-- First correctness theorem:
-- The new ellipse formed by the new parameters, is the old ellipse rotated by ψ.
-- Note: follows from above, require explicit calculation of rotated δx
theorem everything_works_out? (p : EllipseParameters) (ψ : ℝ) :
    p.ellipse.rotate ψ = (p.rotate ψ).ellipse := by
  sorry

-- Second correctness theorem:
-- Even the parametrised point (at t=t) is correctly cotated.
-- Unfortunately we're getting into dependent type hell here, I'll figure something out afterwards...
theorem everything_works_out! (p : EllipseParameters) (ψ t : ℝ) :
    (p.rotate ψ).eval t = cast (congrArg Point <| everything_works_out? p ψ) ((p.eval t).rotate ψ) :=
  sorry

[ZhiKai Pong]

@Kenny Lau I've tried to fill in as much as I can, and I'm happy with everything before EllipseParameters.rotate (Proofs can probably be shortened, and I've changed the dimensions of the ellipse back to a and b instead of E₀x and E₀y for clarity).

For EllipseParameters.rotateI'm still not sure because it seems that I'll still have to give explicit transformations for δx and δy but I don't really think that's tractable. The textbook way of doing this is via what I named as parametrization_eq, which I'm not sure whether it's doable prely through Point.rotate. on the other hand, I got my first wanted result rotate_to_standard without using the parameter or point rotation, I'll try and see whether I can also get the others to work.

[ZhiKai Pong]

I think this is what I want at the end as the goal of this project:

lemma parametrization_eq (p s : EllipseParameters) (ψ : ℝ)
    (hs : IsStandard s.ellipse) (hr : s.ellipse = p.ellipse.rotate ψ) :
    (s.eval t).x = ((p.eval t).rotate ψ).x := by
  -- need `everything_works_out!`
  sorry

lemma aux1 {p s : EllipseParameters} {ψ : ℝ}
    (hs : IsStandard s.ellipse) (hr : s.ellipse = p.ellipse.rotate ψ) :
    s.a ^ 2 + s.b ^ 2 = p.a ^ 2 + p.b ^ 2 := by
  sorry

lemma aux2 {p s : EllipseParameters} {ψ : ℝ}
    (hs : IsStandard s.ellipse) (hr : s.ellipse = p.ellipse.rotate ψ) :
    s.a * s.b = p.a * p.b * sin (p.δy - p.δx) := by
  sorry

lemma _root_.Real.two_mul_tan_sq_div_one_add_tan_sq {x : ℝ} (hx : cos x ≠ 0) :
    2 * tan x / (1 + tan x ^ 2) = sin (2 * x) := by
  rw [tan_eq_sin_div_cos, sin_two_mul]
  field_simp
  rw [pow_two]
  ring

lemma ellipticity_result (p s : EllipseParameters) (ψ χ : ℝ)
    (hs : IsStandard s.ellipse) (hr : s.ellipse = p.ellipse.rotate ψ)
    (hχ : χ = Complex.arg (s.a + Complex.I * s.b)) :
    sin (2 * χ) = 2 * (p.a * p.b * sin (p.δy - p.δx)) / (p.a ^ 2 + p.b ^ 2) := by
  rw [← aux1 hs hr, ← aux2 hs hr]
  trans 2 * tan χ / (1 + tan χ ^ 2)
  rw [two_mul_tan_sq_div_one_add_tan_sq]
  sorry -- figure out how to deal with s.a = 0 and cos χ = 0 case
  have htan : tan χ = s.b/s.a := by rw [hχ, Complex.tan_arg]; simp
  rw [hχ]
  field_simp
  sorry

[ZhiKai Pong]

if the idea here doesn't work, another way for me to retrieve the major and minor axes from E:GeneralEllipse instead of relying on rotation is just to hard code
image.png
which may still give what I want purely through algebraic manipulations.

[ZhiKai Pong]

Aha I think I got it without having to define point rotation (well it is already implicitly defined in the Cartesian Equation, rotation on the parametrized form is too difficult)

lemma AC_invariant (p : EllipseParameters) (ψ : ℝ) :
    (p.ellipse.rotate ψ).A + (p.ellipse.rotate ψ).C = p.ellipse.A + p.ellipse.C := by
  simp
  calc
    _ = p.a ^ 2 * (sin ψ ^ 2 + cos ψ ^ 2) + p.b ^ 2 * (sin ψ ^ 2 + cos ψ ^ 2) := by ring
    _ = p.b ^ 2 + p.a ^ 2 := by simp [add_comm]

lemma F_invariant (p : EllipseParameters) (ψ : ℝ) :
    4 * (p.ellipse.rotate ψ).A * (p.ellipse.rotate ψ).C - (p.ellipse.rotate ψ).B ^ 2 =
    4 * p.ellipse.A * p.ellipse.C - p.ellipse.B ^ 2 := by
  simp
  calc
    _ = 4 * p.a ^ 2 * p.b ^ 2 * (sin ψ ^ 2 + cos ψ ^ 2) ^ 2 -
        4 * p.a ^ 2 * p.b ^ 2 * (sin ψ ^ 2 + cos ψ ^ 2) ^ 2 * cos (p.δy - p.δx) ^ 2 := by ring
    _ = 4 * p.b ^ 2 * p.a ^ 2 - (2 * p.a * p.b * cos (p.δy - p.δx)) ^ 2 := by simp; ring

lemma aux1 {p s : EllipseParameters} {ψ : ℝ}
    (hr : s.ellipse = p.ellipse.rotate ψ) :
    s.a ^ 2 + s.b ^ 2 = p.a ^ 2 + p.b ^ 2 := by
  have hAC : s.ellipse.A + s.ellipse.C = p.ellipse.A + p.ellipse.C := by
    rw [hr]
    rw [AC_invariant p ψ]
  simp at hAC
  linarith

lemma aux2 {p s : EllipseParameters} {ψ : ℝ}
    (hs : IsStandard s.ellipse) (hr : s.ellipse = p.ellipse.rotate ψ) :
    s.a * s.b = p.a * p.b * sin (p.δy - p.δx) ∨ s.a * s.b = - p.a * p.b * sin (p.δy - p.δx) := by
  have hF : 4 * s.ellipse.A * s.ellipse.C - s.ellipse.B ^ 2 =
      4 * p.ellipse.A * p.ellipse.C - p.ellipse.B ^ 2 := by
    rw [hr]
    rw [F_invariant p ψ]
  simp at hF
  have hsp : s.IsStandard := s.isStandard_iff.mp hs
  rw [EllipseParameters.isStandard_def] at hsp
  have hs0 : 2 * s.a * s.b * cos (s.δy - s.δx) = 0 := by
    rcases hsp with hx | hy | hδ
    · simp; left; left; exact hx
    · simp; left; right; exact hy
    · simp; right; exact hδ
  rw [hs0]at hF
  simp at hF
  rw [← mul_one (4 * p.b ^ 2 * p.a ^ 2), ← sin_sq_add_cos_sq (p.δy - p.δx)] at hF
  ring_nf at hF
  have h : (s.a * s.b) ^ 2  = (p.a * p.b * sin (p.δy - p.δx)) ^ 2 := by
    linarith
  rw [sq_eq_sq_iff_eq_or_eq_neg] at h
  simp [h]

hopefully I'll be able to derive what I need from this. Many thanks again, but let me know if you have any additional comments/ideas