import Mathlib.Analysis.Complex.Basic
import Mathlib.Analysis.Complex.TaylorSeries
import Mathlib.Analysis.Calculus.LineDeriv.Basic
import Mathlib.Analysis.Calculus.ContDiff.Defs
import Mathlib.Analysis.Calculus.FDeriv.Basic
import Mathlib.Analysis.Calculus.FDeriv.Symmetric
import Mathlib.Analysis.SpecialFunctions.Complex.LogDeriv
import Mathlib.Data.Complex.Module
import Mathlib.Data.Complex.Order
import Mathlib.Data.Complex.Exponential
import Mathlib.Data.Fin.Tuple.Basic
import Mathlib.Analysis.RCLike.Basic
import Mathlib.Topology.Algebra.InfiniteSum.Module
import Mathlib.Topology.Instances.RealVectorSpace
import Nevanlinna.cauchyRiemann
import Nevanlinna.laplace
import Nevanlinna.partialDeriv

variable {F : Type*} [NormedAddCommGroup F] [NormedSpace ℝ F]


def Harmonic (f : ℂ → F) : Prop :=
  (ContDiff ℝ 2 f) ∧ (∀ z, Complex.laplace f z = 0)


theorem holomorphic_is_harmonic {f : ℂ → ℂ} (h : Differentiable ℂ f) :
  Harmonic f := by

  -- f is real C²
  have f_is_real_C2 : ContDiff ℝ 2 f :=
    ContDiff.restrict_scalars ℝ (Differentiable.contDiff h)

  have fI_is_real_differentiable : Differentiable ℝ (partialDeriv ℝ 1 f) := by
    exact (partialDeriv_contDiff ℝ f_is_real_C2 1).differentiable (Submonoid.oneLE.proof_2 ℕ∞)

  constructor
  · -- f is two times real continuously differentiable
    exact f_is_real_C2

  · -- Laplace of f is zero
    unfold Complex.laplace
    rw [CauchyRiemann₄ h]

    -- This lemma says that partial derivatives commute with complex scalar
    -- multiplication. This is a consequence of partialDeriv_compContLin once we
    -- note that complex scalar multiplication is continuous ℝ-linear.
    have : ∀ v, ∀ s : ℂ, ∀ g : ℂ → ℂ, Differentiable ℝ g → partialDeriv ℝ v (s • g) = s • (partialDeriv ℝ v g) := by
      intro v s g hg

      -- Present scalar multiplication as a continuous ℝ-linear map. This is
      -- horrible, there must be better ways to do that.
      let sMuls : ℂ →L[ℝ] ℂ :=
      {
        toFun := fun x ↦ s * x
        map_add' := by
          intro x y
          ring
        map_smul' := by
          intro m x
          simp
          ring
      }

      -- Bring the goal into a form that is recognized by
      -- partialDeriv_compContLin.
      have : s • g = sMuls ∘ g := by rfl
      rw [this]

      rw [partialDeriv_compContLin ℝ hg]
      rfl

    rw [this]
    rw [partialDeriv_comm f_is_real_C2 Complex.I 1]
    rw [CauchyRiemann₄ h]
    rw [this]
    rw [← smul_assoc]
    simp

    -- Subgoals coming from the application of 'this'
    -- Differentiable ℝ (Real.partialDeriv 1 f)
    exact fI_is_real_differentiable
    -- Differentiable ℝ (Real.partialDeriv 1 f)
    exact fI_is_real_differentiable


theorem re_of_holomorphic_is_harmonic {f : ℂ → ℂ}  (h : Differentiable ℂ f) :
  Harmonic (Complex.reCLM ∘ f) := by

  constructor
  · -- Continuous differentiability
    apply ContDiff.comp
    exact ContinuousLinearMap.contDiff Complex.reCLM
    exact ContDiff.restrict_scalars ℝ (Differentiable.contDiff h)
  · rw [laplace_compContLin]
    simp
    intro z
    rw [(holomorphic_is_harmonic h).right z]
    simp
    exact ContDiff.restrict_scalars ℝ (Differentiable.contDiff h)


theorem im_of_holomorphic_is_harmonic {f : ℂ → ℂ}  (h : Differentiable ℂ f) :
  Harmonic (Complex.imCLM ∘ f) := by

  constructor
  · -- Continuous differentiability
    apply ContDiff.comp
    exact ContinuousLinearMap.contDiff Complex.imCLM
    exact ContDiff.restrict_scalars ℝ (Differentiable.contDiff h)
  · rw [laplace_compContLin]
    simp
    intro z
    rw [(holomorphic_is_harmonic h).right z]
    simp
    exact ContDiff.restrict_scalars ℝ (Differentiable.contDiff h)


theorem logabs_of_holomorphic_is_harmonic
  {f : ℂ → ℂ}
  (h₁ : Differentiable ℂ f)
  (h₂ : ∀ z, f z ≠ 0)
  (h₃ : ∀ z, f z ∈ Complex.slitPlane) :
  Harmonic (fun z ↦ Real.log ‖f z‖) := by

  /- We start with a number of lemmas on regularity of all the functions involved -/

  -- The norm square is real C²
  have normSq_is_real_C2 : ContDiff ℝ 2 Complex.normSq := by
    unfold Complex.normSq
    simp
    conv =>
      arg 3
      intro x
      rw [← Complex.reCLM_apply, ← Complex.imCLM_apply]
    apply ContDiff.add
    apply ContDiff.mul
    apply ContinuousLinearMap.contDiff Complex.reCLM
    apply ContinuousLinearMap.contDiff Complex.reCLM
    apply ContDiff.mul
    apply ContinuousLinearMap.contDiff Complex.imCLM
    apply ContinuousLinearMap.contDiff Complex.imCLM

  -- f is real C²
  have f_is_real_C2 : ContDiff ℝ 2 f :=
    ContDiff.restrict_scalars ℝ (Differentiable.contDiff h₁)

  -- Complex.log ∘ f is real C²
  have log_f_is_holomorphic : Differentiable ℂ (Complex.log ∘ f) := by
    intro z
    apply DifferentiableAt.comp
    exact Complex.differentiableAt_log (h₃ z)
    exact h₁ z

  -- Real.log |f|² is real C²
  have t₄ : ContDiff ℝ 2 (Real.log ∘ ⇑Complex.normSq ∘ f) := by
    rw [contDiff_iff_contDiffAt]
    intro z
    apply ContDiffAt.comp
    apply Real.contDiffAt_log.mpr
    simp
    exact h₂ z
    apply ContDiff.comp_contDiffAt z normSq_is_real_C2
    exact ContDiff.contDiffAt f_is_real_C2

  have t₂ : Complex.log ∘ ⇑(starRingEnd ℂ) ∘ f = Complex.conjCLE ∘ Complex.log ∘ f := by
    funext z
    unfold Function.comp
    rw [Complex.log_conj]
    rfl
    exact Complex.slitPlane_arg_ne_pi (h₃ z)

  constructor
  · -- logabs f is real C²
    have : (fun z ↦ Real.log ‖f z‖) = (2 : ℝ)⁻¹ • (Real.log ∘ Complex.normSq ∘ f) := by
      funext z
      simp
      unfold Complex.abs
      simp
      rw [Real.log_sqrt]
      rw [div_eq_inv_mul (Real.log (Complex.normSq (f z))) 2]
      exact Complex.normSq_nonneg (f z)
    rw [this]

    have : (2 : ℝ)⁻¹ • (Real.log ∘ Complex.normSq ∘ f) = (fun z ↦ (2 : ℝ)⁻¹ • ((Real.log ∘ ⇑Complex.normSq ∘ f) z)) := by
      exact rfl
    rw [this]
    apply ContDiff.const_smul
    exact t₄

  · -- Laplace vanishes
    have : (fun z ↦ Real.log ‖f z‖) = (2 : ℝ)⁻¹ • (Real.log ∘ Complex.normSq ∘ f) := by
      funext z
      simp
      unfold Complex.abs
      simp
      rw [Real.log_sqrt]
      rw [div_eq_inv_mul (Real.log (Complex.normSq (f z))) 2]
      exact Complex.normSq_nonneg (f z)
    rw [this]
    rw [laplace_smul]
    simp

    have : ∀ (z : ℂ), Complex.laplace (Real.log ∘ ⇑Complex.normSq ∘ f) z = 0 ↔ Complex.laplace (Complex.ofRealCLM ∘ Real.log ∘ ⇑Complex.normSq ∘ f) z = 0 := by
      intro z
      rw [laplace_compContLin]
      simp
      -- ContDiff ℝ 2 (Real.log ∘ ⇑Complex.normSq ∘ f)
      exact t₄
    conv =>
      intro z
      rw [this z]

    have : Complex.ofRealCLM ∘ Real.log ∘ Complex.normSq ∘ f = Complex.log ∘ Complex.ofRealCLM ∘ Complex.normSq ∘ f := by
      unfold Function.comp
      funext z
      apply Complex.ofReal_log
      exact Complex.normSq_nonneg (f z)
    rw [this]

    have : Complex.ofRealCLM ∘ ⇑Complex.normSq ∘ f = ((starRingEnd ℂ) ∘ f) * f := by
      funext z
      simp
      exact Complex.normSq_eq_conj_mul_self
    rw [this]

    have : Complex.log ∘ (⇑(starRingEnd ℂ) ∘ f * f) = Complex.log ∘ ⇑(starRingEnd ℂ) ∘ f + Complex.log ∘ f := by
      unfold Function.comp
      funext z
      simp
      rw [Complex.log_mul_eq_add_log_iff]

      have : Complex.arg ((starRingEnd ℂ) (f z)) = - Complex.arg (f z) := by
        rw [Complex.arg_conj]
        have : ¬ Complex.arg (f z) = Real.pi := by
          exact Complex.slitPlane_arg_ne_pi (h₃ z)
        simp
        tauto
      rw [this]
      simp
      constructor
      · exact Real.pi_pos
      · exact Real.pi_nonneg
      exact (AddEquivClass.map_ne_zero_iff starRingAut).mpr (h₂ z)
      exact h₂ z
    rw [this]
    rw [laplace_add]

    rw [t₂, laplace_compCLE]
    intro z
    simp
    rw [(holomorphic_is_harmonic log_f_is_holomorphic).2 z]
    simp

    -- ContDiff ℝ 2 (Complex.log ∘ f)
    exact ContDiff.restrict_scalars ℝ (Differentiable.contDiff log_f_is_holomorphic)

    -- ContDiff ℝ 2 (Complex.log ∘ ⇑(starRingEnd ℂ) ∘ f)
    rw [t₂]
    apply ContDiff.comp
    exact ContinuousLinearEquiv.contDiff Complex.conjCLE
    exact ContDiff.restrict_scalars ℝ (Differentiable.contDiff log_f_is_holomorphic)

    -- ContDiff ℝ 2 (Complex.log ∘ f)
    exact ContDiff.restrict_scalars ℝ (Differentiable.contDiff log_f_is_holomorphic)

    -- ContDiff ℝ 2 (Real.log ∘ ⇑Complex.normSq ∘ f)
    exact t₄