{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FLIM Data analysis\n",
    "\n",
    "## Code blocks for FLIM data phasor analysis\n",
    "\n",
    "First, we try to load some calibration fluorophore of known lifetime (Rhodamine B, Erytrosine B and the mixture) and visulize them in the phasor plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary packeges\n",
    "import numpy as np\n",
    "from sdtfile import SdtFile\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the sdt file\n",
    "file_path = \"./static/sdt/calibration/rhodamineB_in-methanol_ex805_1024x1024_px180_em610_05AOM_05_internal_001.sdt\"\n",
    "# Load the sdt file and print the number of channels\n",
    "flim_file = SdtFile(file_path)\n",
    "print(\"The number of channels is: \", len(flim_file.data))\n",
    "# Load the data from first channel\n",
    "data = flim_file.data[0]\n",
    "print(\"The shape of the data is (x, y, time bins): \", data.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an intensity image of the FLIM data and plot it\n",
    "intensity_image = np.sum(data, axis=2)\n",
    "# Plot the intensity image\n",
    "plt.imshow(intensity_image, cmap='hot')\n",
    "# Add colorbar\n",
    "plt.colorbar()\n",
    "# Add title\n",
    "plt.title(\"Intensity image\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bin the data in x, y to increase the photon count\n",
    "# Bin factor\n",
    "bin_factor = 8\n",
    "# Bin the data according to the bin factor\n",
    "binned_image = intensity_image.reshape((intensity_image.shape[0] // bin_factor, bin_factor,\n",
    "                                  intensity_image.shape[1] // bin_factor, bin_factor)).sum(axis=(1, 3))\n",
    "# Plot the binned image\n",
    "plt.imshow(binned_image, cmap='hot')\n",
    "plt.colorbar()\n",
    "plt.title(\"Binned intensity image\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the flourescence lifetime decay for a pixel and binned data\n",
    "# Create a time array for the x-axis\n",
    "time = np.linspace(0, 12.5, data.shape[2])\n",
    "# Binned data according to the bin factor and keep the third dimension\n",
    "data_binned = data.reshape((data.shape[0] // bin_factor, bin_factor,\n",
    "                            data.shape[1] // bin_factor, bin_factor,\n",
    "                            data.shape[2])).sum(axis=(1, 3))\n",
    "# Select a pixel from the data with all the time bins\n",
    "data_pixel = data[512, 512, :]\n",
    "# Select a pixel from the binned data\n",
    "data_binned_pixel = data_binned[512 // bin_factor, 512 // bin_factor, :]\n",
    "# Plot the results\n",
    "fig, ax = plt.subplots(2, 1)\n",
    "ax[0].scatter(time, data_pixel, s = 2)\n",
    "ax[0].set_title(\"No binning\")\n",
    "ax[0].set_xlabel(\"Time (ns)\"), ax[0].set_ylabel(\"Photon frequency\")\n",
    "ax[1].scatter(time, data_binned_pixel, s = 2)\n",
    "# ax[1].set_yscale(\"log\")\n",
    "ax[1].set_title(\"Binning factor \" + str(bin_factor))\n",
    "ax[1].set_xlabel(\"Time (ns)\"), ax[1].set_ylabel(\"Photon frequency\")\n",
    "fig.tight_layout() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the index of the maximum value in decay data\n",
    "max_index = np.argmax(data_binned_pixel)\n",
    "\n",
    "# Slice the vector to keep only the values on the right side behind the peak\n",
    "sliced_data = data_binned_pixel[max_index:]\n",
    "sliced_time = time[max_index:]\n",
    "sliced_time = sliced_time - sliced_time[0]\n",
    "\n",
    "# Plot the sliced data\n",
    "plt.scatter(sliced_time, sliced_data, s = 2)\n",
    "plt.xlabel(\"Time (ns)\"), plt.ylabel(\"Photon frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Laser frequency input\n",
    "f_laser = 80e6\n",
    "omega = 2 * np.pi * f_laser\n",
    "\n",
    "# Time vector\n",
    "t = sliced_time*1e-9\n",
    "\n",
    "# Calculate the phasor coordinates\n",
    "g = np.sum(sliced_data * np.cos(omega * t)) / np.sum(sliced_data)\n",
    "s = np.sum(sliced_data * np.sin(omega * t)) / np.sum(sliced_data)\n",
    "\n",
    "# Plot the phasor plot\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(g, s, s=35)\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_title('Phasor Plot')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the universal circle \n",
    "center_x = 0.5\n",
    "center_y = 0\n",
    "radius = 0.5\n",
    "x_corr = np.linspace(-radius + center_x, radius + center_x, 100)\n",
    "y_corr = radius * np.sqrt(1 - ((x_corr - center_x) / radius) ** 2) + center_y\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.plot(x_corr, y_corr, 'k')\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_title('Phasor Plot')\n",
    "ax.axes.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the bound (3.4 ns) and free (0.4) NADH points to the phasor plot\n",
    "tau = [3.4, 0.4]\n",
    "g_NADH = []\n",
    "s_NADH = []\n",
    "# Calculate the phasor coordinates for the NADH points\n",
    "for tau_m in tau: \n",
    "    s_now = (2*np.pi*f_laser*tau_m*1e-9)/(1 + (2*np.pi*f_laser)**2*((tau_m*1e-9)**2))\n",
    "    g_now = 1/(1 + (2*np.pi*f_laser)**2*((tau_m*1e-9)**2))\n",
    "    g_NADH.append(g_now)\n",
    "    s_NADH.append(s_now)\n",
    "NADH = np.vstack((g_NADH, s_NADH))\n",
    "# Plot the phasor plot with the NADH points\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(g, s, s = 35, c='b')\n",
    "ax.plot(x_corr, y_corr, 'k')\n",
    "ax.scatter(g_NADH, s_NADH, s=15, c='r')\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_title('Phasor Plot')\n",
    "# Plot text in x, y coordinates\n",
    "ax.text(g_NADH[0]-0.15, s_NADH[0], '3.4 ns', color='k', fontsize=12)\n",
    "ax.text(g_NADH[1]-0.15, s_NADH[1], '0.4 ns', color='k', fontsize=12)\n",
    "ax.axes.set_aspect('equal')\n",
    "# According to the literature, the Rhodamine B has a lifetime of ~ 2.2 ns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the phasor coordinates for all the point in the image\n",
    "g_all = []\n",
    "s_all = []\n",
    "\n",
    "for idx_x in range(data_binned.shape[0]):\n",
    "    for idx_y in range(data_binned.shape[1]):\n",
    "        pixel_decay = data_binned[idx_x, idx_y, :]\n",
    "        # Find the index of the maximum value in decay data\n",
    "        max_index = np.argmax(pixel_decay)   \n",
    "        # Slice the vector to keep only the values on the right side behind the peak\n",
    "        sliced_data = pixel_decay[max_index:]\n",
    "        sliced_time = time[max_index:]\n",
    "        sliced_time = sliced_time - sliced_time[0]\n",
    "        # Time vector\n",
    "        t = sliced_time*1e-9\n",
    "        g = np.sum(sliced_data * np.cos(omega * t)) / np.sum(sliced_data)\n",
    "        s = np.sum(sliced_data * np.sin(omega * t)) / np.sum(sliced_data)\n",
    "        g_all.append(g)\n",
    "        s_all.append(s)\n",
    "# Plot the phasor plot with the NADH points\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(g_all, s_all, s = 2, c='b', alpha=0.1)\n",
    "ax.plot(x_corr, y_corr, 'k')\n",
    "ax.scatter(g_NADH, s_NADH, s=15, c='r')\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_title('Phasor Plot')\n",
    "# Plot text in x, y coordinates\n",
    "ax.text(g_NADH[0]-0.15, s_NADH[0], '3.4 ns', color='k', fontsize=12)\n",
    "ax.text(g_NADH[1]-0.15, s_NADH[1], '0.4 ns', color='k', fontsize=12)\n",
    "ax.axes.set_aspect('equal') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load erytrosine B data and plot the monoexponential decay in phasor\n",
    "# Select the sdt file\n",
    "file_path = \"./static/sdt/calibration/erytrosineB_in-methanol_ex750_1024x1024_px180_em610_10AOM_10_internal_001.sdt\"\n",
    "# Load the sdt file and print the number of channels\n",
    "flim_file = SdtFile(file_path)\n",
    "print(\"The number of channels is: \", len(flim_file.data))\n",
    "# Load the data from first channel\n",
    "data = flim_file.data[0]\n",
    "print(\"The shape of the data is (x, y, time bins): \", data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the flourescence lifetime decay for a pixel and binned data\n",
    "# Create a time array for the x-axis\n",
    "time = np.linspace(0, 12.5, data.shape[2])\n",
    "# Binned data according to the bin factor and keep the third dimension\n",
    "data_binned = data.reshape((data.shape[0] // bin_factor, bin_factor,\n",
    "                            data.shape[1] // bin_factor, bin_factor,\n",
    "                            data.shape[2])).sum(axis=(1, 3))\n",
    "# Calculate the phasor coordinates for all the point in the image\n",
    "g1_all = []\n",
    "s1_all = []\n",
    "\n",
    "for idx_x in range(data_binned.shape[0]):\n",
    "    for idx_y in range(data_binned.shape[1]):\n",
    "        pixel_decay = data_binned[idx_x, idx_y, :]\n",
    "        # Find the index of the maximum value in decay data\n",
    "        max_index = np.argmax(pixel_decay)   \n",
    "        # Slice the vector to keep only the values on the right side behind the peak\n",
    "        sliced_data = pixel_decay[max_index:]\n",
    "        sliced_time = time[max_index:]\n",
    "        sliced_time = sliced_time - sliced_time[0]\n",
    "        # Time vector\n",
    "        t = sliced_time*1e-9\n",
    "        g = np.sum(sliced_data * np.cos(omega * t)) / np.sum(sliced_data)\n",
    "        s = np.sum(sliced_data * np.sin(omega * t)) / np.sum(sliced_data)\n",
    "        g1_all.append(g)\n",
    "        s1_all.append(s)\n",
    "# Plot the phasor plot with the NADH points\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(g1_all, s1_all, s = 2, c='b', alpha=0.1)\n",
    "ax.scatter(g_all, s_all, s = 2, c='g', alpha=0.1)\n",
    "ax.plot(x_corr, y_corr, 'k')\n",
    "ax.scatter(g_NADH, s_NADH, s=15, c='r')\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_title('Phasor Plot')\n",
    "# Plot text in x, y coordinates\n",
    "ax.text(g_NADH[0]-0.15, s_NADH[0], '3.4 ns', color='k', fontsize=12)\n",
    "ax.text(g_NADH[1]-0.15, s_NADH[1], '0.4 ns', color='k', fontsize=12)\n",
    "ax.axes.set_aspect('equal')\n",
    "# According to the literature, the Erythrosin B has a lifetime of ~ 0.48 ns. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mixture of rhodamineB(10%) and erytrosineB(90%) in methanol\n",
    "# Load erytrosine B data and plot the monoexponential decay in phasor\n",
    "# Select the sdt file\n",
    "file_path = \"./static/sdt/calibration/rhodamineB(1)+erytrosineB(9)_in-methanol_ex750_1024x1024_px180_em610_10AOM_10_internal_001.sdt\"\n",
    "# Load the sdt file and print the number of channels\n",
    "flim_file = SdtFile(file_path)\n",
    "print(\"The number of channels is: \", len(flim_file.data))\n",
    "# Load the data from first channel\n",
    "data = flim_file.data[0]\n",
    "print(\"The shape of the data is (x, y, time bins): \", data.shape)\n",
    "# Plot the flourescence lifetime decay for a pixel and binned data\n",
    "# Create a time array for the x-axis\n",
    "time = np.linspace(0, 12.5, data.shape[2])\n",
    "# Binned data according to the bin factor and keep the third dimension\n",
    "data_binned = data.reshape((data.shape[0] // bin_factor, bin_factor,\n",
    "                            data.shape[1] // bin_factor, bin_factor,\n",
    "                            data.shape[2])).sum(axis=(1, 3))\n",
    "# Calculate the phasor coordinates for all the point in the image\n",
    "g2_all = []\n",
    "s2_all = []\n",
    "\n",
    "for idx_x in range(data_binned.shape[0]):\n",
    "    for idx_y in range(data_binned.shape[1]):\n",
    "        pixel_decay = data_binned[idx_x, idx_y, :]\n",
    "        # Find the index of the maximum value in decay data\n",
    "        max_index = np.argmax(pixel_decay)   \n",
    "        # Slice the vector to keep only the values on the right side behind the peak\n",
    "        sliced_data = pixel_decay[max_index:]\n",
    "        sliced_time = time[max_index:]\n",
    "        sliced_time = sliced_time - sliced_time[0]\n",
    "        # Time vector\n",
    "        t = sliced_time*1e-9\n",
    "        g = np.sum(sliced_data * np.cos(omega * t)) / np.sum(sliced_data)\n",
    "        s = np.sum(sliced_data * np.sin(omega * t)) / np.sum(sliced_data)\n",
    "        g2_all.append(g)\n",
    "        s2_all.append(s)\n",
    "# Plot the phasor plot with the NADH points\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(g1_all, s1_all, s = 2, c='b', alpha=0.1)\n",
    "ax.scatter(g_all, s_all, s = 2, c='g', alpha=0.1)\n",
    "ax.scatter(g2_all, s2_all, s = 2, c='c', alpha=0.1)\n",
    "ax.plot(x_corr, y_corr, 'k')\n",
    "ax.scatter(g_NADH, s_NADH, s=15, c='r')\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_title('Phasor Plot')\n",
    "# Plot text in x, y coordinates\n",
    "ax.text(g_NADH[0]-0.15, s_NADH[0], '3.4 ns', color='k', fontsize=12)\n",
    "ax.text(g_NADH[1]-0.15, s_NADH[1], '0.4 ns', color='k', fontsize=12)\n",
    "ax.axes.set_aspect('equal') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bound/free NADH is usually related to the oxidative phosphorillation or glycolysis. The conclusion is made based on the position of the transformed data in the phasor plot in relation to the free and bound components.\n",
    "\n",
    "<div>\n",
    "<img src=\"./static/images/FLIM/metabolic_trajcectory.png\" alt=\"Drawing\" style=\"width: 590px;\"/>\n",
    "</div>\n",
    "\n",
    "<em> Created by Stringari, C. et al. PLOS ONE, 2012. [Phasor Fluorescence Lifetime Microscopy of Free and Protein-Bound NADH Reveals Neural Stem Cell Differentiation Potential](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0048014)</em>\n",
    "\n",
    "Second, we load the image of the brain cells and display the data in phasor plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mixture of rhodamineB(10%) and erytrosineB(90%) in methanol\n",
    "# Load erytrosine B data and plot the monoexponential decay in phasor\n",
    "# Select the sdt file\n",
    "file_path = \"./static/sdt/240906_mouse_001_non_tumor/1024x1024x256_FLIM_NADH_750nm_472-30.sdt\"\n",
    "# Load the sdt file and print the number of channels\n",
    "flim_file = SdtFile(file_path)\n",
    "print(\"The number of channels is: \", len(flim_file.data))\n",
    "# Load the data from first channel\n",
    "data = flim_file.data[0]\n",
    "print(\"The shape of the data is (x, y, time bins): \", data.shape)\n",
    "# Plot the flourescence lifetime decay for a pixel and binned data\n",
    "# Create a time array for the x-axis\n",
    "time = np.linspace(0, 12.5, data.shape[2])\n",
    "# Binned data according to the bin factor and keep the third dimension\n",
    "data_binned = data.reshape((data.shape[0] // bin_factor, bin_factor,\n",
    "                            data.shape[1] // bin_factor, bin_factor,\n",
    "                            data.shape[2])).sum(axis=(1, 3))\n",
    "# create image to display\n",
    "intensity_image = np.sum(data_binned, axis=2)\n",
    "\n",
    "# Calculate the phasor coordinates for all the point in the image\n",
    "g_cells_all = []\n",
    "s_cells_all = []\n",
    "\n",
    "for idx_x in range(data_binned.shape[0]):\n",
    "    for idx_y in range(data_binned.shape[1]):\n",
    "        pixel_decay = data_binned[idx_x, idx_y, :]\n",
    "        # Find the index of the maximum value in decay data\n",
    "        max_index = np.argmax(pixel_decay)   \n",
    "        # Slice the vector to keep only the values on the right side behind the peak\n",
    "        sliced_data = pixel_decay[max_index:]\n",
    "        sliced_time = time[max_index:]\n",
    "        sliced_time = sliced_time - sliced_time[0]\n",
    "        # Time vector\n",
    "        t = sliced_time*1e-9\n",
    "        g = np.sum(sliced_data * np.cos(omega * t)) / np.sum(sliced_data)\n",
    "        s = np.sum(sliced_data * np.sin(omega * t)) / np.sum(sliced_data)\n",
    "        g_cells_all.append(g)\n",
    "        s_cells_all.append(s)\n",
    "# Plot the phasor plot with the NADH points\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.imshow(intensity_image, cmap='gray')\n",
    "ax.axis('off')\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(g_cells_all, s_cells_all, s = 2, c='orange', alpha=0.1)\n",
    "ax.plot(x_corr, y_corr, 'k')\n",
    "ax.plot([g_NADH[0], g_NADH[1]], [s_NADH[0], s_NADH[1]], 'gray')\n",
    "ax.scatter(g_NADH, s_NADH, s=15, c='r')\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_title('Phasor Plot')\n",
    "# Plot text in x, y coordinates\n",
    "ax.text(g_NADH[0]-0.15, s_NADH[0], '3.4 ns', color='k', fontsize=12)\n",
    "ax.text(g_NADH[1]-0.15, s_NADH[1], '0.4 ns', color='k', fontsize=12)\n",
    "ax.axes.set_aspect('equal') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a 2D histogram of the g and s values\n",
    "hist, xedges, yedges = np.histogram2d(g_cells_all, s_cells_all, bins=50)\n",
    "\n",
    "# Construct arrays for the anchor positions of the bars.\n",
    "xpos, ypos = np.meshgrid(xedges[:-1] + 0.01, yedges[:-1] + 0.01, indexing=\"ij\")\n",
    "xpos = xpos.ravel()\n",
    "ypos = ypos.ravel()\n",
    "zpos = 0\n",
    "\n",
    "# Construct arrays with the dimensions for the bars.\n",
    "dx = dy = 0.02 * np.ones_like(zpos)\n",
    "dz = hist.ravel()\n",
    "\n",
    "# Filter out the zero values\n",
    "nonzero = dz > 0\n",
    "xpos = xpos[nonzero]\n",
    "ypos = ypos[nonzero]\n",
    "dz = dz[nonzero]\n",
    "\n",
    "# Create the 3D plot\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "# Use a colormap to color the bars based on their height\n",
    "colors = plt.cm.viridis(dz / max(dz))\n",
    "\n",
    "ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color=colors, zsort='average')\n",
    "ax.plot(x_corr, y_corr, 'k')\n",
    "ax.scatter(g_NADH, s_NADH, s=15, c='r')\n",
    "ax.set_xlim([0, 1])\n",
    "ax.set_ylim([0, 0.5])\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_zlabel('Frequency')\n",
    "ax.set_title('3D Phasor plot with frequency distribution')\n",
    "ax.set_box_aspect([1, 1, 0.5])  # Aspect ratio is 1:1:0.5\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the center of mass of the phasor plot\n",
    "g_center = np.sum(g_cells_all) / len(g_cells_all)\n",
    "s_center = np.sum(s_cells_all) / len(s_cells_all)\n",
    "\n",
    "# Plot the phasor plot with the NADH points\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(g_cells_all, s_cells_all, s = 2, c='orange', alpha=0.1)\n",
    "ax.scatter(g_center, s_center, s = 35, c='orange', edgecolors='black')\n",
    "ax.plot(x_corr, y_corr, 'k')\n",
    "ax.scatter(g_NADH, s_NADH, s=15, c='r')\n",
    "ax.set_xlabel('G')\n",
    "ax.set_ylabel('S')\n",
    "ax.set_title('Phasor Plot')\n",
    "# Plot text in x, y coordinates\n",
    "ax.text(g_NADH[0]-0.15, s_NADH[0], '3.4 ns', color='k', fontsize=12)\n",
    "ax.text(g_NADH[1]-0.15, s_NADH[1], '0.4 ns', color='k', fontsize=12)\n",
    "ax.axes.set_aspect('equal') \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the ratio of bound NADH\n",
    "# The quantum yield of bound NADH is 6x higher than free NADH\n",
    "qy = 6\n",
    "d_bound = np.sqrt((g_center - g_NADH[0])**2 + (s_center - s_NADH[0])**2)\n",
    "d_free = np.sqrt((g_center - g_NADH[1])**2 + (s_center - s_NADH[1])**2)\n",
    "free_nadh_fraction = (d_bound / qy) / ((d_bound / qy) + d_free)\n",
    "print(f\"The ratio of bound NADH is: {free_nadh_fraction*100:.2f}\", \"%\")   "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}