{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots, LinearAlgebra\n",
    "include(\"p_space.jl\")\n",
    "\n",
    "g0(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(-1 + exp(p*q*R^2))*R)/(2*sqrt(π))\n",
    "g1(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(2 + p*q*R^2 + exp(p*q*R^2)*(-2 + p*q*R^2)))/(2*p*sqrt(π)*q*R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vertices = (0, 0.5 - 0.3im, 1, 6)\n",
    "subdivisions = 64\n",
    "p, w = get_mesh(vertices, subdivisions)\n",
    "\n",
    "# resonance example from my thesis\n",
    "V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)\n",
    "\n",
    "basis_eig = eigen(get_H_matrix(V_basis, p, w))\n",
    "basis = basis_eig.vectors .* w\n",
    "\n",
    "evals = basis_eig.values\n",
    "E_target = 0.7\n",
    "E = evals[argmin(norm.(evals .- E_target))]\n",
    "\n",
    "print(\"E = $E\")\n",
    "scatter(real.(evals), imag.(evals), xlim = (0,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ResonanceEC: Eq. (20)\n",
    "V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))\n",
    "\n",
    "H = get_H_matrix(V_system(0.45), p, w)\n",
    "H_berggren = transpose(basis) * H * basis\n",
    "\n",
    "evals = eigvals(H)\n",
    "scatter(real.(evals), imag.(evals), xlim = (0, 0.5))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.0",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}