#+PROPERTY: header-args :session gaussian_farts :kernel python :pandoc yes :async yes

#+begin_src jupyter-python
  import scipy
  import numpy as np
  import matplotlib.pyplot as plt


  from scipy.special import expi
#+end_src

#+RESULTS:

* Inverse Laplace Transform of the Propagator
We need the laplace transform of our spectral density.
#+begin_src jupyter-python
  η = 1
  ω_c = 1
  Ω = 1

  def α(t):
      return η / np.pi * (ω_c / (1 + 1j * ω_c * t)) ** 2


  def α_l(z):
      return (
          (1j * ω_c * η)
          / np.pi
          ,* (1 + 1j * z / ω_c * np.exp(-1j * z / ω_c) * expi(1j * z / ω_c))
      )
  def α_l_num(z):
      return quadpy.quad(lambda t: α(t) * np.exp(-t*z), 0, 1000)
#+end_src

#+RESULTS:

Now let's try calculating \(K_{11}\).
First we define its laplace transform
#+begin_src jupyter-python
  def K_l(z):
      return np.array([[z, Ω * np.ones_like(z)], [-(Ω + np.imag(α_l(z))), z]]) / (z ** 2 + Ω ** 2 + Ω * np.imag(α_l(z)))
#+end_src

#+RESULTS:

Plotting it may help to give us some intuition.
#+begin_src jupyter-python
  res = np.linspace(.1, 10, 100)
  ims = np.linspace(-10, 10, 100)
  zs = np.meshgrid(res, 1j* ims, indexing="ij")
  zs = (zs[0] + zs[1])
  image = plt.imshow(K_l(zs)[0,1].imag)
  plt.colorbar(image, ax=plt.gca(), location='right', anchor=(0, 0.3), shrink=0.7)
#+end_src

#+RESULTS:
:RESULTS:
: <matplotlib.colorbar.Colorbar at 0x7fa8d2f60c70>
[[file:./.ob-jupyter/ed423eb6eb8ae253acdfbf34ebc5c0318e43b0ca.svg]]
:END:


Let's do the inverse transform.
#+begin_src jupyter-python
  import quadpy

  def integrand(s, t):
      return 1 / (2j * np.pi) * K_l(s) * np.exp(s * t)

  def γ(t):
      return 1j * t + 10 + 1.1 * t**2/np.sqrt(t**2 + 1)

  def γ_dot(t):
      return 1.1 * (2*t/np.sqrt(t**2 + 1) - (t/np.sqrt(t**2 + 1))**3) + 1j

  def K(t, lim):
      return quadpy.quad(
          lambda s: integrand(γ(s), t) * γ_dot(s),
          -lim,
          lim,
          epsabs=1e-3,
          epsrel=1e-3,
      )[0]


  t = np.linspace(0, 1, 100)
  #plt.plot(t, [K11(τ, 200, .1).real for τ in t])
  #plt.plot(t, [K11(τ, 200, 10).real for τ in t])
  #plt.plot(t, [K11(τ, 200, 20).real for τ in t])
  #plt.plot(t, [K11(τ, 400).real for τ in t])
  from scipy import interpolate
  K_i = interpolate.interp1d(t, [K(τ, 300) for τ in t], axis=0)
#+end_src

#+RESULTS:


#+begin_src jupyter-python
  t = np.linspace(0, 1, 1000)
  Ks = K_i(t).real
  for i, j in np.meshgrid([0, 1], [0, 1]):
      plt.plot(t, Ks[:, i, j])
#+end_src

#+RESULTS:
[[file:./.ob-jupyter/f30e618e0e7ac37663f69ca4fbebc83317e671ed.svg]]