From fa2dc35365c2fc6ceee6f09c2349cd66cec92d25 Mon Sep 17 00:00:00 2001
From: ysyapa <ysyapa@ncsu.edu>
Date: Thu, 18 Jan 2024 23:12:46 +0000
Subject: [PATCH] Basic resonance solution

---
 p_space.jl         | 27 ++++++++++++++++++++
 p_space_test.ipynb | 61 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 88 insertions(+)
 create mode 100644 p_space.jl
 create mode 100644 p_space_test.ipynb

diff --git a/p_space.jl b/p_space.jl
new file mode 100644
index 0000000..a763c53
--- /dev/null
+++ b/p_space.jl
@@ -0,0 +1,27 @@
+using FastGaussQuadrature
+
+function gausslegendre_shifted(a, b, n)
+    scale = (b - a) / 2
+    shift = (a + b) / 2
+    p, w = gausslegendre(n)
+    p = p .* scale .+ shift
+    w = w .* scale
+    return (p, w)
+end
+
+function get_mesh(vertices, subdivisions)
+    p = Vector{ComplexF64}()
+    w = Vector{ComplexF64}()
+    for (a, b) in zip(vertices, vertices[2:end])
+        p_new, w_new = gausslegendre_shifted(a, b, subdivisions)
+        append!(p, p_new)
+        append!(w, w_new)
+    end
+    return (p, w)
+end
+
+get_V_matrix(V_pq, p, w) = V_pq.(p, transpose(p)) .* transpose(w)
+
+get_K_matrix(p, μ) = collect(Diagonal(p.*p ./ (2*μ)))
+
+get_H_matrix(V_pq, p, w, μ=0.5) = get_K_matrix(p, μ) + get_V_matrix(V_pq, p, w)
diff --git a/p_space_test.ipynb b/p_space_test.ipynb
new file mode 100644
index 0000000..c1c79b1
--- /dev/null
+++ b/p_space_test.ipynb
@@ -0,0 +1,61 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using Plots, LinearAlgebra\n",
+    "include(\"p_space.jl\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vertices = (0, 1 - 0.5im, 2, 6)\n",
+    "subdivisions = 64\n",
+    "p, w = get_mesh(vertices, subdivisions)\n",
+    "\n",
+    "scatter(real.(p), imag.(p))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "g1(α, p, q) = sqrt(α/(4*π))*((1/α+2/(p*q))*exp(-(p+q)^2/(4*α))+(1/α-2/(p*q))*exp(-(p-q)^2/(4*α)))\n",
+    "V_pq(p, q) = -10 * g1(1, p, q)\n",
+    "\n",
+    "H = get_H_matrix(V_pq, p, w)\n",
+    "evals = eigen(H).values\n",
+    "\n",
+    "E_target = 0.258 - 0.164im\n",
+    "E = evals[argmin(norm.(evals .- E_target))]\n",
+    "\n",
+    "print(\"E = $E\")\n",
+    "scatter(real.(evals), imag.(evals), xlim = (0,1))"
+   ]
+  }
+ ],
+ "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
+}