Pushed rendered site to repo

Author: Soybean44

.gitignore

@@ -1,2 +1 @@
 /.quarto/
-/_site/

_site/about.html

@@ -0,0 +1,58 @@
+[
+  {
+    "objectID": "about.html",
+    "href": "about.html",
+    "title": "About",
+    "section": "",
+    "text": "About this site"
+  },
+  {
+    "objectID": "experiments/voltage/index.html",
+    "href": "experiments/voltage/index.html",
+    "title": "Poisson’s Equation (Voltage)",
+    "section": "",
+    "text": "Import required libraries\n\nimport matplotlib.pyplot as plt\nfrom matplotlib.animation import FuncAnimation\nimport numpy as np\nfrom IPython.display import HTML\n\nSet up the data required\n\nsize = 50\nvoltages = np.zeros((size, size))\nV0 = 50\ncond = -0.5\nvoltages[10, 10] = V0\n\nThe partial differential equation we are solving is \\[\n  \\nabla V = -\\frac{1}{\\epsilon_0}\\rho\n\\]\nthe right hand side is cond. Then we can solve this with the boundaries of a 50x50 having a voltage of 0. To solve this we can use the following formula for \\(u_{i,j}\\) \\[\nu_{i,j} = \\frac{1}{4}(u_{i,j+1}+u_{i+1,j}+u_{i,j-1}+u_{i-1,j}-cond)\n\\]\n\nfor _ in range(500):\n    for i in range(size):\n        for j in range(size):\n            u1 = voltages[i, j + 1] if j + 1 < size else 0\n            u2 = voltages[i, j - 1] if j - 1 > -1 else 0\n            u3 = voltages[i + 1, j] if i + 1 < size else 0\n            u4 = voltages[i - 1, j] if i - 1 > -1 else 0\n            voltages[i, j] = (u1 + u2 + u3 + u4 - cond) / 4\n\nThen we get the gradient of this to find the electric field\n\npos = np.array(range(size))\n\ndy, dx = np.gradient(-voltages, pos, pos)\n\nThen we graph it\n\nskip = 5  # Number of points to skip\nplt.figure()\nplt.imshow(\n    np.abs(voltages),\n    extent=(0, size, 0, size),\n    origin=\"lower\",\n    cmap=\"viridis\",\n)\nplt.colorbar(label=\"Voltage\")\nplt.quiver(\n    pos[::skip],\n    pos[::skip],\n    dx[::skip, ::skip],\n    dy[::skip, ::skip],\n    color=\"r\",\n    headlength=3,\n)\nplt.title(\"Scalar Field using contour\")\nplt.xlabel(\"X-axis\")\nplt.ylabel(\"Y-axis\")\n\nplt.show()"
+  },
+  {
+    "objectID": "experiments/wave_equation/index.html",
+    "href": "experiments/wave_equation/index.html",
+    "title": "Wave Equation",
+    "section": "",
+    "text": "Import required libraries\n\nimport matplotlib.pyplot as plt\nfrom matplotlib.animation import FuncAnimation\nimport numpy as np\nfrom IPython.display import HTML\n\nSet up the values for the wave equation\n\nL = 2*np.pi # Length of the string\nNx = 100 # Number of spatial points\nNt = 100 # Number of time steps\ndx = L/(Nx-1) # Spatial step size\nv = 5 # Wave speed\ndt = np.sqrt(1)*dx/v # Time step size\n\nprint(f\"dt: {dt}s\")\n# Initialize the spring\nx = np.linspace(0, L, Nx)\ny = np.zeros(Nx)\ny_old = np.zeros(Nx)\n\n# initial disturbance\ny_old = np.sin(x)\n# Apply fixed boundary conditions\ny_old[0] = y_old[Nx-1] =  0\ny[0] = y[Nx-1] =  0\n\n#  make sure the C calculated comes out to 1\nC = (v*dt/dx)**2\nprint(C)\n\ndt: 0.012693303650867852s\n1.0\n\n\nThe partial differential equation we are solving is \\[\n\\frac{\\partial^2 y}{\\partial^2 t} = \\frac{T}{\\mu} \\frac{\\partial^2 y}{\\partial^2 x}\n\\] Then to solve the differential equation we can use the formula \\[\ny_{i+1,j}= 2y_{i,j} -y_{i-1,j} + C(y_{i,j+1} - 2y_{i,j} + y_{i,j-j})\n\\] where the i index represents time and the j index represents position. C is defined by \\[\nC=(T/\\mu)(\\Delta t)^2/(\\Delta x)^2\n\\] for the first timestep, we dont have a previous time, so instead we use our initial velocity condition \\(\\frac{\\partial y}{\\partial t} = 0\\) at \\(t=0\\) which in turn we can use the forward difference formula and the previous formula for the differential equation to gain a new formula \\[\ny_{i,j}= y_{i-1,j} + \\frac{1}{2}C(y_{i-1,j+1} - 2y_{i-1,j} + y_{i-1,j+1j})\n\\] One important thing to note is that for the numerical approximation to be stable, \\(C<=1\\) where the theoretical solution is when equivalence holds\n\n# compute initial timestep assuming initial velocity is 0\ny[1:-1] = y_old[1:-1] + 0.5 * C * (y_old[2:] - 2 * y_old[1:-1] + y_old[:-2])\ndef update_string(y_old, y):\n  y_new = np.zeros(Nx)\n  # Use np slices to compute everything in parallel, \n  y_new[1:-1] = 2*y[1:-1] - y_old[1:-1] + C*(y[2:] - 2*y[1:-1] + y[:-2])\n  # Update arrays for next iteration\n  y_old[:] = y[:]\n  y[:] = y_new[:]\n\nFinally we can plot the data and also view it as an animation\n\n# Set up the plot\nfig, ax = plt.subplots()\nline, = ax.plot(x, y)\nax.set_ylim(-1.5, 1.5)\nax.set_xlim(0, L)\nax.set_title('Vibrating String')\nax.set_xlabel('Position (m)')\nax.set_ylabel('Displacement (m)')\n# Animation function\nskip = 2\ndef animate(frame):\n    global y_old, y, y_new, skip\n    for _ in range(skip):\n      update_string(y_old, y)\n    line.set_ydata(y)\n    return line,\n# Create animation\nani = FuncAnimation(fig, animate, frames=Nt, blit=True, interval=dt*1000*skip)\nani.save('../../videos/standing_wave.mp4', writer='ffmpeg', fps=1/(dt*skip))\nHTML(ani.to_jshtml())\n\n\n\n\n\n\n\n\n  \n  \n    \n    \n      \n          \n      \n        \n      \n          \n      \n          \n      \n          \n      \n          \n      \n          \n      \n          \n      \n          \n    \n    \n      \n      Once\n      \n      Loop\n      \n      Reflect"
+  },
+  {
+    "objectID": "experiments/schrodingers/index.html",
+    "href": "experiments/schrodingers/index.html",
+    "title": "Schrödinger’s Equation",
+    "section": "",
+    "text": "Start with imports and setup\nimport scipy.sparse as sparse\nfrom skimage import measure\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nimport plotly.graph_objects as go\nimport numpy as np\nimport torch\nimport time\nfrom datetime import timedelta\nfrom IPython.display import Image\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')"
+  },
+  {
+    "objectID": "experiments/schrodingers/index.html#theory",
+    "href": "experiments/schrodingers/index.html#theory",
+    "title": "Schrödinger’s Equation",
+    "section": "Theory",
+    "text": "Theory\nSchrödinger’s Equation is the equation governing motion at the quantum scale. We have a wavefunction \\(\\psi(x)\\) which satisfies the following differential equation \\[\n-\\frac{\\hbar}{2m}\\nabla^2 \\psi(x) + V(x) \\psi(x) = E \\psi(x)\n\\]\nWe can then see that this differential equation actually forms an eigenvector problem or more specifically, \\[\n\\hat{H} \\ket{\\psi} = (\\hat{T} + \\hat{V})\\ket{\\psi} = E \\ket{\\psi}\n\\]\nwe can then find the central difference of \\(\\frac{d^2 \\psi}{dx^2}\\), being \\[\n\\frac{d^2 \\psi}{dx^2} = \\frac{\\psi_{x-1}-2\\psi_x+\\psi_{x+1}}{(\\Delta x)^2}\n\\]\nThe next step is to recognize our boundary conditions, which for us is \\(\\psi(0)=\\psi(L)=0\\). Then we can represent \\(\\hat{H}\\) as a discrete matrix instead of an arbitrary linear operator, given our finite difference as \\[\n\\hat{H} = -\\frac{\\hbar}{2m(\\Delta x)^2}D+V\n\\] Where \\(D\\) and \\(V\\) are matrices representing the Laplacian and Potential, represented as such \\[\nD = \\begin{pmatrix}\n  -2 & 1 & 0 \\\\\n  1 & -2 & 1 \\\\\n  0 & 1 & -2 \\\\\n\\end{pmatrix}, V= \\begin{pmatrix}\n  V_1 & 0 & 0 \\\\\n  0 & V_2 & 0 \\\\\n  0 & 0 & V_3 \\\\\n\\end{pmatrix}\n\\] Given we have three datapoints each diagonal entry in \\(V\\) is the voltage at that given position. These matrices are then “extended” to fit within the number of datapoints we have, making sure to keep the matrix hemertian.\nHowever what if we want to extend this to higher dimensions. Then we need to take what is known as a kroneker sum of the matrix \\(D\\) and use it in its place. It is defined as so \\[\nD \\oplus D = D \\otimes I + I \\otimes D\n\\] where \\(\\otimes\\) is the kroneker product and is defined as \\[\nD \\otimes I =\n\\begin{pmatrix}\nD & 0 & 0\\\\\n0 & D & 0\\\\\n0 & 0 & D\\\\\n\\end{pmatrix}\n\\] then for 2 dimensions we just use one kroneker sum and for 3 dimensions we can use 3 as so \\[\nD \\oplus D \\oplus D\n\\] which happens to work out. Then we can flatten our grid of sample points and make it a simple n dimensional vector our final equation then becomes \\[\n\\left[\\frac{\\hbar}{2m}(D\\oplus D\\oplus D)+V\\right]\\vec{\\psi} = E \\vec{\\psi}\n\\] which we can simplify using natural units, setting \\(\\hbar = 1\\) and as well making the entries in \\(V\\) return \\(m\\delta^2V(x,y,z)\\) as well as the energy eigenvalues being now \\(m\\delta^2E\\), where \\(\\delta\\) is a generalized finite differential for all directions in space. This becomes the following equation then \\[\n\\left[\\frac{1}{2}(D\\oplus D\\oplus D)+V\\right]\\vec{\\psi} = m\\delta^2E \\vec{\\psi}\n\\] Noting that the matrix \\(V\\) itself is returning a different value now of \\(m\\delta^2V(x)\\)"
+  },
+  {
+    "objectID": "experiments/schrodingers/index.html#calculations",
+    "href": "experiments/schrodingers/index.html#calculations",
+    "title": "Schrödinger’s Equation",
+    "section": "Calculations",
+    "text": "Calculations\nWith this we can finally start on the calculations. Start by defining a Potential (Noting this is really potenetial energy and not electric potential). The potential is then the potential energy of an electron in an atom of hydrogen which is \\[\nV(r) = -\\frac{e^2}{4 \\pi \\epsilon_0 r} = -\\frac{1}{ma_0 r},\\, a_0 = \\frac{4\\pi\\epsilon_0\\hbar^2}{e^2m}\\approx 5.29 \\times 10^-11m\n\\] This allows us to define the radius in terms of the bohr radius \\(a_0\\). We can then rearange the equation to solve for \\(m\\delta^2V(x)\\) and write our function as \\[\nm\\delta^2V(x)= -\\frac{\\delta^2}{a_0 r} = -\\frac{(\\delta/a_0)^2}{(r/a_0)}\n\\] Changing the units to a_0 by dividing by \\(a_0/a_0\\)\nStart by defining values for our grid and delta\n\nN = 120\n# Note that both the coordinates and delta are measured in units of a_0\nX,Y,Z = np.mgrid[-25:25:N*1j,-25:25:N*1j,-25:25:N*1j] # create a grid from -25 to 25 a_0\ndelta = np.diff(X[:,0,0])[0] # get the step size\n\nNow we write a function to get the potential as follows\n\ndef get_potential(x,y,z):\n  # since the units are in a_0 we can treat the a_0 in the equation defined above as 1\n  # As well we add 1e-10 to the value of r to avoid any infinities at r=0\n  return -delta**2/np.sqrt(x**2+y**2+z**2 + 1e-10)\n\nV0 = get_potential(X,Y,Z) # This creates a 3d grid of potentials\nV0.shape # shape should be n*n*n we really want this to be a diagonal matrix in M_{n^3}(R)\n\n(120, 120, 120)\n\n\nWe want to get our proper \\(V\\) matrix which we can do along with our \\(D\\) matrix\n\ndiag = np.ones([N]) # create a vector of ones in R^n\ndiags = np.array([diag, -2*diag, diag]) # create a n*3 matrix representing the diagonals in D\nD = sparse.spdiags(diags, np.array([-1,0,1]), N, N) # Create D from the diags\nT = -1/2 * sparse.kronsum(sparse.kronsum(D,D), D) # Create T from D\nV = sparse.diags(V0.reshape(N**3), (0)) # Create a sparse matrix from V0\nH = T+V # Finally create the hamiltonian matrix\nH.shape # Print out the shape of H which should be (N^3,N^3)\n\n(1728000, 1728000)\n\n\nConvert H so that it is on the gpu. This makes finding the eigenvalues easier\n\nH = H.tocoo()\nH = torch.sparse_coo_tensor(indices=torch.tensor([H.row, H.col]), values=torch.tensor(H.data), size=H.shape).to(device)\n\n/tmp/ipykernel_7383/2973072149.py:2: UserWarning:\n\nCreating a tensor from a list of numpy.ndarrays is extremely slow. Please consider converting the list to a single numpy.ndarray with numpy.array() before converting to a tensor. (Triggered internally at /build/python-pytorch/src/pytorch/torch/csrc/utils/tensor_new.cpp:254.)\n\n\n\nFinally we can compute the eigenvalues/eigenvectors as follows\n\nnum_states = 20 # number of states to calculate\nstart = time.perf_counter()\neigenvalues, eigenvectors = torch.lobpcg(H, k=num_states, largest=False)\nend = time.perf_counter()\nelapsed = timedelta(seconds = end-start)\nprint(\"Time to calculate:\", elapsed)\n\nTime to calculate: 0:49:01.725120\n\n\nHere is a function to get a an eigenvector given an energy level\n\ndef get_e(n):\n    return eigenvectors.T[n].reshape((N,N,N)).cpu().numpy()\n\nWe can compute an object given a certain energy level and plot it\n\neigenstate = 6 # Eigenstate to display\nverts, faces, _, _ = measure.marching_cubes(get_e(eigenstate)**2, 1e-6, spacing=(0.1, 0.1, 0.1))\nintensity = np.linalg.norm(verts, axis=1)\n\nfig = go.Figure(data=[go.Mesh3d(x=verts[:, 0], y=verts[:, 1], z=verts[:, 2], \n                                i=faces[:, 0], j=faces[:, 1], k=faces[:, 2],\n                                intensity=intensity,\n                                colorscale='Agsunset',\n                                opacity=0.5)])\n\nfig = fig.update_layout(scene=dict(xaxis=dict(visible=False),\n                              yaxis=dict(visible=False),\n                              zaxis=dict(visible=False),\n                              bgcolor='rgb(0, 0, 0)'),\n                  margin=dict(l=0, r=0, b=0, t=0))\n\n\nfig.show()\n\n        \n        \n        \n\n\n                            \n                                            \n\n\nplot the eigenvalues by first converting the units to joules\n\nhbar = 1.055e-34\na0 = 5.29e-11\nm = 9.11e-31 # Mass of an electron\nJ_to_eV = 6.242e18 # convert from Joules to eV\nconversion = hbar**2 / m / delta**2 / a0**2 * J_to_eV \nconverted_eigenvalues = eigenvalues.cpu() * conversion\n\nplt.figure()\nplt.scatter(list(range(num_states)), converted_eigenvalues, label=\"Enegy Eigenstates\")\nplt.xlim(0,20)\nplt.title(\"Eigenstates of a Hydrogen atom (Orbitals)\")\nplt.xlabel(\"Energy Level\")\nplt.ylabel(\"Energy Eigenvalue (eV)\")\nplt.legend()\nplt.show()"
+  },
+  {
+    "objectID": "index.html",
+    "href": "index.html",
+    "title": "Computational Physics",
+    "section": "",
+    "text": "This is a Quarto website.\nTo learn more about Quarto websites visit https://quarto.org/docs/websites.\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nCpp_test\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nPoisson's Equation (Voltage)\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nSchrödinger's Equation\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nWave Equation\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\nNo matching items"
+  },
+  {
+    "objectID": "experiments/cpp_test/index.html",
+    "href": "experiments/cpp_test/index.html",
+    "title": "Cpp_test",
+    "section": "",
+    "text": "This should test the render cpp script.\nextern void render_rect(int x, int y, int w, int h);\n\nint main() {\n  render_rect(0,0,50,50);\n  return 0;\n}\nNow we display a canvas with a square"
+  }
+]