Total_Energy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8638d3b4-533b-4749-a5db-5d2be3127c17",
   "metadata": {},
   "outputs": [],
   "source": [
    "from array import array\n",
    "import numpy as np\n",
    "from array import array\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7fa13e8d-cad0-490c-9888-16e1bf80d4fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "dirin='/home/rdavis39/CLIM470/CLIM-470-Numerical-Weather-Prediction/'\n",
    "energy_file ='totalen.dat'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "51967b9a-491c-4742-8126-54676b8ba566",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Variables\n",
    "ntime = 3    # Number of time steps\n",
    "nx = 48      # Number of grid points in x-direction\n",
    "ny = 16      # Number of grid points in y-direction\n",
    "recl=(nx*ny)*4\n",
    "total_energy=np.zeros((ntime,nx,ny))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "29689b81-058d-4933-bb44-740cda134829",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_total_energy = open(dirin+energy_file,'rb')\n",
    "for t in range(ntime):\n",
    "    tmp = ds_total_energy.read(recl)\n",
    "    tmp1 = array('f',tmp)\n",
    "    #tmp2=tmp1[1:-1] #for fortran sequential\n",
    "    total_energy[t,:,:]=np.reshape(tmp1,(nx,ny))\n",
    "del tmp,tmp1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "105eb157-05e8-4ba1-9cc6-9741920f75f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Total energy at each time step\n",
    "total_energy_all = total_energy.sum(axis=(1, 2))  # Sum over spatial dimensions (ny, nx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cc62ee6e-bac1-4a55-bfba-ad698268a44a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAAEWCAYAAACtyARlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAst0lEQVR4nO3dfZxUdfn/8de1uyx3CyKKZqCChSIgN66YSCYoBCka3iVWJn1L0krNSgszs/RrVuT3q2TxTe2HFoqpiHd5h4KZisqdgBgq5g1gKSiwy82yN9fvj3N2mZ2dnT0Le3YPM+/n4zGPnTm31zVnZq79nPM555i7IyIiIslS0NYBiIiISEMq0CIiIgmkAi0iIpJAKtAiIiIJpAItIiKSQCrQIiIiCaQCLTnBzOab2Td3cd6DzKzczApbOq4I63Yz+3Rrr7c5zOwrZvZEW8fRGDObYWbXtnUcIi1NBVoSw8zeNrNtYbH8d/jDWxLTekbXvnb3d929xN2rY1jXAWZ2m5m9b2ZlZvZPM/u5mXVu6XXFxd1nuvvnmzufmU0Pt2W5me0ws8qU149mma/e9tkdZjbJzKpT1lv7+GRLLF8kTirQkjSnuHsJMAQYCkxp23B2nZl1B14AOgLD3b0LMAboBnyqDUNrFe5+QfiPTwlwHXB37Wt3/0IrhvJCynprH+tacgVmVtSSyxMBFWhJKHf/N/A4QaEGwMyOMbPnzWyjmb1iZiMzzWtmnzKzp81sg5mtN7OZZtYtHPdn4CDgobAldbmZ9Q53NReZ2dlmtjBteZea2YPh8/ZmNtXM3jWz/4StxI6NpPF9oAz4qru/Heb1nrtf4u7LUqYbbWZvhHndbGbWVB7h+LfN7IdmtszMNpnZ3WbWIWX85WHLfZ2ZfTN1d3pz8ghbof9Iee1mdkGmmKMys1PN7NVw/vlmdng4vMH2CYffE+5V2WRmfzezAc1ZX5Y4mnoPx5vZ0jDO581sUNq8PzKzZcCW8PPzNTN7J9xmP63dG2BmnzCzrWa2T8r8R5rZh2bWriVykdyjAi2JZGa9gC8Ab4avewKPANcC3YEfAveZWY9MswO/BD4JHA4cCFwN4O7nAu8SttTd/ddp8z4EHGZmfVOGfRm4M3x+PXAowT8OnwZ6Alc1ksZoYLa71zSR7nhgGDAI+BIwtqk8UnwJGAf0CeefBGBm4wj+QRgdxjkybb7m5NGcmJtkZocCdwHfA3oAfyMoyMVZts+jQF9gP2AxMLMZsTalsfdwKPAn4FvAPsD/AQ+aWfuUec8BTibYK3Io8HvgK8ABwF4E72vtP5zzw3XVOheY5e6VLZiL5BAVaEmaOWZWBrwHfAD8LBz+VeBv7v43d69x9yeBhcBJ6Qtw9zfd/Ul3r3D3D4EbgOOjrNzdtwIPEPzwEhbqfgQ/zAZMBi5194/cvYxg1+3ERha3D/B+hNVe7+4b3f1dYB7hXoOIedzk7uvc/SOCfy6GhMO/BPw/d381zOnq2hl2IY/IMUd0NvBImFslMJXgMMCxjc3g7n9y9zJ3rwhzGWxme0Vc3zFhC7j2sTptfGPv4WTg/9z9RXevdvfbgQrgmLR533P3bcCZwEPu/g9330HwD0/qzQ5uJ/gcY0GHxHOAP0fMQfKQCrQkzYTwWO1IgsK4bzj8YOCs1B9a4LMELZV6zGx/M5tlZmvNbDPwl5TlRHEnYYEmaD3PCYtcD6ATsCglhsfC4ZlsyBRfBv9Oeb4VKGlGHhnnJWh1v5cyLvV5c/OIHHNEnwTeqX0R7mF4j7C1mc7MCs3sejNbHb4Pb4ejom7TBe7eLeWRfvy/sVwOBn6Q9pk7MIy/Vur7Wu89Dz8zG1LGPwD0N7M+BH0RNrn7SxFzkDykAi2J5O7PADMIWlcQ/PD9Oe2HtrO7X59h9usIWi5HuHtXglZL6jHSpm7h9iTQw8yGEBTq2t3b64FtwICUGPYKO0FlMhc4zcx29XvWVB7ZvA/0Snl9YMrz5ubR0tYRFD+grkV/ILA2HJS+fb4MfJFgd/1eQO/aWWONMvjM/XfaZ66Tu9+VMk1qrPXe8/CYft0xZ3ffDvyVYDuei1rP0gQVaEmy/wXGmNlggtbjKWY2NmxRdTCzkeGx6nRdgHJgU3js+rK08f8BDmlspeFu13uA3xAc734yHF4D3AL8j5ntB8GxcTNr7PjrDUBX4HYzOzhl+htSOxtl0VQe2fwV+LqZHW5mnYCfpuTX3Dxa2l+Bk83sxLCD1A8Idh0/H45P3z5dwvEbCFr+17VSnLcAF5jZZyzQ2cxONrMujUx/L8Fn9FgzKybYFZ/+T8QdBMe4T0UFWpqgAi2JFR53vQO4yt3fI2hFXQF8SNC6uYzMn+GfA0cCmwg6ls1OG/9L4Mpwt+UPG1n9nQQttnvcvSpl+I8IOq4tCHe3zgUOayT+jwiOq1YCL4bH1p8K43ozS+pR82iUuz8K3ERwfPhNYEE4qqK5ebQ0d19F0IqcRtCaP4WgU9iOcJL07XMHwS7xtcDKlFyiGm4Nz4MeFiHOhcD5wO+Ajwner0lZpn8VuAiYRdCaLifoR1GRMs1zQA2w2N3fybQckVrm3tTePhHZ01lwGtMKoH3aPxwSEwsusrMR6Ovu/0oZ/jRwp7vf2laxyZ5BLWiRHGVmp1lwvvPewK8IehirOMfIzE4xs04WXCluKrCcnZ3aCFvuRwJ3t02EsieJtUCHJ+kvt+BE/4VNzyEiLehbBLtYVwPVwIVtG05e+CJBJ7h1BOdtT/RwN6WZ3U5wKOF74altIlnFuovbzN4GjnL39bGtREREJAdpF7eIiEgCxd2C/hdB70cnuCLPHzNMM5ngij107Nix9MADD0yfZJfV1NRQUJAb/4PkSi65kgcol6TKlVxyJQ9QLtm8/vrr690980WC3D22B9Az/Lsf8ArwuWzTl5aWekuaN29eiy6vLeVKLrmSh7tySapcySVX8nBXLtkAC72RmhjrvzTuvjb8+wFwP3B0nOsTERHJFbEV6PCqO11qnwOfJzgPU0RERJoQ503G9wfuDy6zSxHBifmPxbg+ERGRnBFbgXb3t4DBcS1fRCQpKisrWbNmDdu3b2+xZe6111689tprLba8tqRcoEOHDvTq1Yt27dpFnifOFrSISF5Ys2YNXbp0oXfv3oR7DXdbWVkZXbo0dl+OPUu+5+LubNiwgTVr1tCnT5/I8+VGv3cRkTa0fft29tlnnxYrzpJbzIx99tmn2XtYVKBFRFqAirNksyufDxVoERGRBFKBFhHJAYWFhQwZMoSBAwdyyimnsHHjxqzTX3311UydOjXrNHPmzGHlypV1r6+66irmzp3bEuEydepU+vXrx5AhQxg2bBh33HEHAL1792b9+mTcvmHJkiV84xvfAOC+++5jwIABjB07lg0bNgCwevVqzj777Lrpd+zYwec+9zmqqlrmpnEq0CIirWzOkrWMuP5p+vz4EUZc/zRzlqzd7WV27NiRpUuXsmLFCrp3787NN9+8+3GmFehf/OIXjB49ereXO336dJ588kleeuklli5dylNPPVV79clEue6667j44osBmDZtGi+//DJf//rXufPOOwG48sorufbaa+umLy4u5sQTT+Tuu1vmbqIq0CIirWjOkrVMmb2ctRu34cDajduYMnt5ixTpWsOHD2ft2mB5q1evZty4cZSWlnLcccfxz3/+s8H0t9xyC8OGDWPw4MGcccYZbN26leeff54HH3yQyy67jCFDhrB69WomTZrEvffey2OPPcZZZ51VN//8+fMZP348AE888QTDhw/nyCOP5KyzzqK8vLzB+q677jr+8Ic/0LVrVwC6du3KeeedVzd+2rRpHHnkkRxxxBF18b700ksMHz6coUOHcuyxx7Jq1SoAZsyYwemnn864cePo27cvl19+ed1ybrvtNg499FCOPvpozj//fL773e8C8OGHH3LGGWcwbNgwhg0bxnPPPdcgxrKyMpYtW8bgwcHZwgUFBVRUVLBt2zbatWvHs88+yyc+8Qn69u1bb74JEyYwc+bMpjZRJDrNSkSkBf38oVdZuW5zo+OXvLuRHdU19YZtq6zm8nuXcddL79YNq66uprCwEID+n+zKz04ZEGn91dXVPPXUU3W7ZidPnsz06dPp27cvL774It/+9rd5+umn681z+umnc/755wNBq/C2227joosu4tRTT2X8+PGceeaZ9aYfPXo0kydPZsuWLXTu3Jm7776biRMnsn79eq699lrmzp1L586d+dWvfsUNN9zApZdeWjfv5s2bKSsr45BDDmk0h3333ZfFixfz+9//nqlTp3LrrbfSr18/nn32WYqKipg7dy5XXHEF9913HwBLly5lyZIltG/fnsMOO4yLLrqIwsJCrrnmGhYvXkyXLl044YQT6ortJZdcwqWXXspnP/tZ3n33XcaOHdvg3OaFCxcycODAutdTpkxh9OjR7LfffsyaNYuzzjqLWbNmNYh94MCBvPzyy01upyhUoEVEWlF6cW5qeFTbtm1jyJAhrF27lsMPP5wxY8ZQXl7O888/X6+1W1FR0WDeFStWcOWVV7Jx40bKy8sZO3Zs1nUVFRUxbtw4HnroIc4880weeeQRfv3rX/PMM8+wcuVKRowYEeS0YwfDhw9vdi6nn346AKWlpcyePRuATZs2cd555/HGG29gZlRWVtZNf+KJJ7LXXnsB0L9/f9555x3Wr1/P8ccfT/fu3QE466yzeP311wGYO3duvV33mzdvpry8nJKSkrph77//Pj167LzJ1JgxYxgzZgxlZWXcf//9nHTSSbz++utMnTqVvffemxtvvJFOnTpRWFhIcXFxi5z7rQItItKCmmrpjrj+adZu3NZgeM9uHbn7WzuLWXN/4GuPQW/dupWxY8dy8803M2nSJLp168bSpUuzzjtp0iTmzJnD4MGDmTFjBvPnz29yfRMnTuR3v/sd3bt356ijjqJLly64O2PGjOGuu+6qN21ZWVnd865du1JSUsJbb73VaCu6ffv2QNDxrbbD1U9/+lNGjRrF/fffz9tvv83IkSMbTJ8+T2NqampYsGABHTp0aHSajh07ZjxveevWrcyYMYPHH3+c8ePHM3v2bO69915mzpxZtxeioqIi67Kj0jFoEZFWdNnYw+jYrrDesI7tCrls7GEtsvxOnTpx00038dvf/pZOnTrRp08f7rnnHiC4otUrr7zSYJ6ysjIOOOAAKisr6x0/7dKlS73imur4449n8eLF3HLLLUycOBGAY445hueee44333wTgC1bttS1WlNNmTKF73znO2zeHBwKKC8vr+vF3ZhNmzbRs2dPIDju3JRhw4bxzDPP8PHHH1NVVVW3Oxzg85//PNOmTat7nekfmMMPP7wuj1Q33ngjF198Me3atWPbtm2YGQUFBWzduhWADRs2sO+++zbrkp6NUYEWEWlFE4b25JenH0HPbh0xgpbzL08/gglDe7bYOoYOHcqgQYO46667mDlzJrfddhuDBw9mwIABPPDAAw2mv+aaa/jMZz7DiBEj6NevX93wiRMn8pvf/IahQ4eyevXqevMUFhYyfvx4Hn300boOYj169GDGjBmcc845DBo0iOHDh2fslHbhhRcyatQohg0bxsCBAznuuOMoKMheji6//HKmTJnC0KFDI53G1LNnT6644gqOPvpoRowYQe/evet2g990000sXLiQQYMG0b9/f6ZPn95g/n79+rFp06Z6/6CsW7eORYsWMWHCBAAuuugihg0bxvTp0/nyl78MwLx58zj55JObjC+Sxm4U3RaP0tLSFr0Rtm4Snjy5koe7ckmqtshl5cqVLb7MzZs3t/gy20pb5VJWVubu7pWVlT5+/HifPXt2s+a/4YYb/JZbbqk3rKlcTjvtNF+1alXGcZk+J8BCb6QmqgUtIiI56eqrr667eEufPn3qWr5RXXjhhfWObzdlx44dTJgwgUMPPbSZkWamTmIiIpKTmrpSWlM6dOjAueeeG3n64uJivva1r+3WOlOpBS0i0gI8gVfCkuTYlc+HCrSIyG7q0KEDGzZsUJGWjDy8H3RzT73SLm4Rkd3Uq1cv1qxZw4cffthiy9y+fXuLnEubBMol+CeuV69ezZpHBVpEZDe1a9eOPn36tOgy58+fz9ChQ1t0mW1Fuewa7eIWERFJIBVoERGRBFKBFhERSSAVaBERkQRSgRYREUkgFWgREZEEUoEWERFJIBVoERGRBFKBFhERSSAVaBERkQRSgRYREUkgFWgREZEEUoEWERFJIBVoERGRBFKBFhERSSAVaBERkQSKvUCbWaGZLTGzh+Nel4iISK5ojRb0JcBrrbAeERGRnBFrgTazXsDJwK1xrkdERCTXmLvHt3Cze4FfAl2AH7r7+AzTTAYmA+y///6ls2bNarH1l5eXU1JS0mLLa0u5kkuu5AHKJalyJZdcyQOUSzajRo1a5O5HZRzp7rE8gPHA78PnI4GHm5qntLTUW9K8efNadHltKVdyyZU83JVLUuVKLrmSh7tyyQZY6I3UxDh3cY8ATjWzt4FZwAlm9pcY1yciIpIzYivQ7j7F3Xu5e29gIvC0u381rvWJiIjkEp0HLSIikkBFrbESd58PzG+NdYmIiOQCtaBFREQSSAVaREQkgVSgRUREEijrMWgzGw58FTgOOADYBqwAHgH+4u6bYo9QREQkDzXagjazR4FvAo8D4wgKdH/gSqAD8ICZndoaQYqIiOSbbC3oc919fdqwcmBx+Pitme0bW2QiIiJ5rNEWdIbijJmdaGanmFm7xqYRERGR3Re5k5iZ/Zbg8p2DgQdii0hEREQa38UdFuRr3H1jOOgg4Evh8+UxxyUiIpLXsrWgZwOzzOxiMysE7gDmAS8At7RGcCIiIvkq2zHo59x9HPARQU9uc/eR7n6Mu9/YahGKiIjkoWynWRWZ2cnAB8AEYLCZPWhmg1srOBERkXyV7TSrOQS7szsBX3H388zsk8AvzMzd/fzWCFBERCQfZSvQB7v7eDMrBhYAuPs64JtmNqQ1ghMREclX2Qr0H83shfD5Dakj3H1pbBGJiIhI4wXa3acB01oxFhEREQll6yR2pZntnWX8CWY2Pp6wRERE8lu2XdzLgYfNbDvBtbc/JLhJRl9gCDAXuC7uAEVERPJRtl3cDxDcsaovwSU+DwA2A38BJrv7ttYJUUREJP9kvR80gLu/AbzRCrGIiIhIKPLNMkRERKT1qECLiIgkkAq0iIhIAjVZoM3sUDN7ysxWhK8HmdmV8YcmIiKSv6K0oG8BpgCVAO6+DJgYZ1AiIiL5LkqB7uTuL6UNq4ojGBEREQlEKdDrzexTgAOY2ZnA+7FGJSIikueaPA8a+A7wR6Cfma0F/gV8NdaoRERE8lyUC5W8BYw2s85AgbuXxR+WiIhIfmuyQJvZ99NeA2wCFum2kyIiIvGIcgz6KOACoGf4+BYwDrjFzC6PMTYREZG8FeUYdC/gSHcvBzCznwGPAJ8DFgG/ji88ERGR/BSlBb0fUJHyuhLYP7ybVUXmWURERGR3RGlBzwReNLMHwtenAHeGncZWxhaZiIhIHovSi/saM3sMODYcdIG7Lwyff6Wx+cysA/B3oH24nnvd/We7Ga+IiEheiNKCBlgMrK2d3swOcvd3m5inAjjB3cvNrB3wDzN71N0X7Hq4IiIi+SHKaVYXAT8D/gNUA0ZwVbFB2eZzdwfKw5ftwofvTrAiIiL5woI6mmUCszeBz7j7hmYv3KyQoKf3p4Gb3f1HGaaZDEwG2H///UtnzZrV3NU0qry8nJKSkhZbXlvKlVxyJQ9QLkmVK7nkSh6gXLIZNWrUInc/KuNId8/6AOYBRU1N18QyuoXLGZhtutLSUm9J8+bNa9HltaVcySVX8nBXLkmVK7nkSh7uyiUbYKE3UhOjHIN+C5hvZo+QclqVu98Q9T8Ed99oZvMILnCyIup8IiIi+SrKedDvAk8CxUCXlEdWZtbDzLqFzzsCY4B/7nKkIiIieSTKaVY/38VlHwDcHh6HLgD+6u4P7+KyRERE8kqUXtw9gMuBAUCH2uHufkK2+dx9GTB0dwMUERHJR1F2cc8k2DXdB/g58DbwcowxiYiI5L0oBXofd78NqHT3Z9z9v4CsrWcRERHZPVF6cVeGf983s5OBdUD3+EISERGRKAX6WjPbC/gBMA3oClwaa1QiIiJ5Lkov7tqe15uAUfGGIyIiIhC9F/f5QO/U6cNj0SIiIhKDKLu4HwCeBeYS3CxDREREYhalQHfyDDe5EBERkfhEOc3qYTM7KfZIREREpE6jLWgzKyO4f7MBV5hZBcEpV0Zwu+eurROiiIhI/mm0QLt7kzfEEBERkXg0uovbzMaa2ZkZhp9hZmPiDUtERCS/ZTsGfRXwTIbhzwC/iCccERERgewFur27f5g+0N3XA53jC0lERESyFeiuZtbgGLWZtQM6xheSiIiIZCvQs4FbzKyutWxmJcD0cJyIiIjEJFuBvhL4D/COmS0ys0XAv4APw3EiIiISk2ynWVUBPzaznwOfDge/6e7bWiUyERGRPBblblbbgOWtEIuIiIiEolzqU0RERFqZCrSIiEgCNVmgLfBVM7sqfH2QmR0df2giIiL5K0oL+vfAcOCc8HUZcHNsEYmIiEik+0F/xt2PNLMlAO7+sZkVxxyXiIhIXovSgq40s0KCW09iZj2AmlijEhERyXNRCvRNwP3Afmb238A/gOtijUpERCTPRTkPemZ4FbETAQMmuPtrsUcmIiKSx5os0GbWHfgAuCtlWDt3r4wzMBERkXwWZRf3YoLrb78OvBE+f9vMFptZaZzBiYiI5KsoBfpJ4CR339fd9wG+ADwMfJvgFCwRERFpYVEK9DHu/njtC3d/Ahju7guA9rFFJiIikseinAf9vpn9CJgVvj4b+E946pVOtxIREYlBlBb0l4FewJzwcVA4rBD4UlyBiYiI5LMop1mtBy5qZPSbLRuOiIiIQLTTrHoAlwMDgA61w939hBjjEhERyWtRdnHPBP4J9AF+DrwNvNzUTGZ2oJnNM7OVZvaqmV2yW5GKiIjkkSgFeh93vw2odPdn3P2/gCit5yrgB+7eHzgG+I6Z9d+NWEVERPJGlF7ctVcMe9/MTgbWAd2bmsnd3wfeD5+XmdlrQE9g5S7GKiIikjfM3bNPYDYeeBY4EJgGdAWudveHIq/ErDfwd2Cgu29OGzcZmAyw//77l86aNavhAnZReXk5JSUlLba8tpQrueRKHqBckipXcsmVPEC5ZDNq1KhF7n5UxpHunvUBjIgyLMv8JcAi4PSmpi0tLfWWNG/evBZdXlvKlVxyJQ935ZJUuZJLruThrlyyARZ6IzUxyjHoaRGHNWBm7YD7gJnuPjvKPCIiIpLlGLSZDQeOBXqY2fdTRnUluEhJVmZmwG3Aa+5+w+4GKiIikk+ytaCLCXZPFwFdUh6bgTMjLHsEcC5wgpktDR8n7Wa8IiIieaHRFrS7PwM8Y2Yz3P2d5i7Y3f8B2O4EJyIikq+inGbV3sz+CPROnd51JTEREZHYRCnQ9wDTgVuB6njDEREREYhWoKvc/Q+xRyIiIiJ1opxm9ZCZfdvMDjCz7rWP2CMTERHJY1Fa0OeFfy9LGebAIS0fjoiIiEC0+0H3aY1AREREZKcmd3GbWSczuzLsyY2Z9Q2vzy0iIiIxiXIM+v8BOwiuKgawFrg2tohEREQkUoH+lLv/mvC2k+6+FV2AREREJFZRCvQOM+tI0DEMM/sUUBFrVCIiInkuSi/unwGPAQea2UyCa2xPijMoERGRfBelF/eTZrYYOIZg1/Yl7r4+9shERETyWJRe3KcRXE3sEXd/GKgyswmxRyYiIpLHohyD/pm7b6p94e4bCXZ7i4iISEyiFOhM00Q5di0iIiK7KEqBXmhmN5jZp8LHDcCiuAMTERHJZ1EK9EUEFyq5G5gFbAe+E2dQIiIi+S7rrmozKwQedvdRrRSPiIiI0EQL2t2rgRoz26uV4hERERGidfYqB5ab2ZPAltqB7n5xbFGJiIjkuSgFenb4EBERkVYS5Upit4fX4j7I3Ve1QkwiIiJ5L8qVxE4BlhJcjxszG2JmD8Ycl4iISF6LcprV1cDRwEYAd18KHBJbRCIiIhKpQFemXuozVBNHMCIiIhKI0knsVTP7MlBoZn2Bi4Hn4w1LREQkv0W9ktgAoAK4E9gEfC/GmERERPJeoy1oM+sAXAB8GlgODHf3qtYKTEREJJ9la0HfDhxFUJy/AExtlYhEREQk6zHo/u5+BICZ3Qa81DohiYiISLYWdGXtE+3aFhERaV3ZWtCDzWxz+NyAjuFrA9zdu8YenYiISJ5qtEC7e2FrBiIiIiI7RTnNSkRERFpZbAXazP5kZh+Y2Yq41iEiIpKr4mxBzwDGxbh8ERGRnBVbgXb3vwMfxbV8ERGRXGbuHt/CzXoDD7v7wCzTTAYmA+y///6ls2bNarH1l5eXU1JS0mLLa0u5kkuu5AHKJalyJZdcyQOUSzajRo1a5O5HZRzp7rE9gN7AiqjTl5aWekuaN29eiy6vLeVKLrmSh7tySapcySVX8nBXLtkAC72Rmqhe3CIiIgmkAi0iIpJAcZ5mdRfwAnCYma0xs2/EtS4REZFck+1Sn7vF3c+Ja9kiIiK5Tru4RUREEkgFWkREJIFUoEVERBJIBVpERCSBVKBFREQSSAVaREQkgVSgRUREEkgFWkREJIFUoEVERBJIBVpERCSBVKBFREQSSAVaREQkgVSgRUREEkgFWkREJIFUoEVERBJIBVpERCSBVKBFREQSSAVaREQkgVSgRUREEkgFWkREJIFUoEVERBJIBVpERCSBVKBFREQSqKitAxAREUm6OUvW8pvHV7F24zZ6Lniay8YexoShPWNdpwq0iIhIFnOWrGXK7OVsq6wGYO3GbUyZvRwg1iKtAi0iIjnJ3dleWUN5RRVbKqooDx9b6v5WpzyvYsuOKsorqinfXsmWiupg+I4q3vtoKzVef9nbKqv5zeOrVKBFRCQ/VFXXBMVxR1W94lm+PbWQVqcV2qDYltU9D4Zv3VFNdXplbUSn4kI6ty+iS/siOrcvonP7Qj7ZrQOd2xfxzoatGedZt3FbS6begAq0iIjsMndnW2V1xhZpbYt12b8qWVz5er3iWVtU01u1FVU1kdbbrtCCQlpcRElYUPfq2I6e3TrQuTgosiXta/8WhkU3tQDvnK9zcREFBdbouha+/TFrMxTjT3bruMvvWxQq0CIieaayuqbebt7ytJbnzuJaXX/4jrRh24NhkRqpq96gc9hKLUlrpZa0L0obXr+olqQ8audrX1QY+/tU67Kxh9U7Bg3QsV0hl409LNb1qkCLiCRcXSt1e4aiuiO1qAbTbKmoqttFvCVDoW1OKzW9eKa3Urt0SCuqxfULbZcORSx+6QXGnjAyays1yWqPM9f14u7WUb24RUT2VKmt1PK0Qlm/xVrdYLotFVV88PFW/Pmn6opwxEOpdC4upKRDSlEtLqJnt+K6FmlJWqu0tkXasBXbcq3UjkW2xxbnWhOG9mTC0J7Mnz+fkSNHtso6VaBFRAhaqVt3VNfb9VtWUdmwp29qUd1RtbPFGhbS2tbtjoit1OLCguA4aErR7NapmOKqAvr02rfBrt8uWYpqp3aFe3whlJ1UoEVkj7WjqiatOFax/MMqtix7v/7x1LrdvdWkd2CqHb5lRxUeoZVqRrh7t7BeK/XA7p12djpqX0RJccMWad3u4JTdwMVFmS/oGLTUBrfwOyZ7EhVokTzW2ldHqqlxtlY27Olb20pNP02mqaK6o7qRVuqixfVeFhcV1OuxW9K+iL07FXPg3p3Clmi7Bp2SUlupqbuDO6qVKq1EBVokT0W9OlJFVXXG02e2ZDh2mn6stV4Hpu1VbK2sjtxKrW2Bpu7KPbBz/VZqw1Nminj91WUcN3xYpFaqSJLFWqDNbBxwI1AI3Oru18e5vlptcc1USQ53xx1q3KkJ/+58HQ6r2fncU6ZLnba6Jvv4tzZW0+WdjxsdXzesJjWO1Glrl+1U12QZX1N/+fVySlt2dRPjaxyqw+cPLF1X77QRCK6O9MN7XmHqE6uabqWmSW2l1rZIu3cuDnb91u3uTTkftUNR2vmqOwtxp+JCzHatlerrCun3ia67NK9IksRWoM2sELgZGAOsAV42swfdfWVc64SmWwWZfrzr/ejV1B9W70evpuEPfd2PXk36j3OGH/ea7D/knvJDnT7+tTWV/Ofld+vW1dgPfeZ1NRzvTcUcxlHdSJw1NZnfw7r11FBv2bXv4eZN25i6/NkG4xssqya9MKZsi5qmC3CrWfB8K64sOzMoMKOg7u/O52ZQULBzmJmxdUd1xuVU1ThH9+6e9XzUTEW1XaFaqSItKc4W9NHAm+7+FoCZzQK+CMRaoH/z+KqMrYLv3b2US/+6NNLutcRasXyXZzWDwvBHO/2HvP6Pd+oPfPBDXlCw8wc/43IKapcTzJM6vrDAaFdgdeN9m7Fvlw7111MQrictptrl7BxfP44CMwoLMuWz83lT4wsyFK4G7024nPTxK1YsZ/CgQTvft4L6hTFTzIUN1tWM8eHyCxvZhs1tcY64/umMV0fq2a0jN5w9ZJc/ayLSMuIs0D2B91JerwE+kz6RmU0GJocvy81s1e6stPgTny5tbNyOf7+5aHeW3cb2Bda3dRAtIFfygD08l4KOXbsXde1xMGY7m77uNe9t/vAdm7L5ozYMbXft0dslRa7kAcolm4MbG9HmncTc/Y/AH+NYtpktdPej4lh2a8uVXHIlD1AuSZUrueRKHqBcdlWcB43WAgemvO4VDhMREZEmxFmgXwb6mlkfMysGJgIPxrg+ERGRnBHbLm53rzKz7wKPE5xm9Sd3fzWu9TUill3nbSRXcsmVPEC5JFWu5JIreYBy2SXme3S3ZhERkdykExdFREQSSAVaREQkgfbIAm1m48xslZm9aWY/zjC+vZndHY5/0cx6p4ybEg5fZWZjWzXwDCLk8n0zW2lmy8zsKTM7OGVctZktDR9t3gEvQi6TzOzDlJi/mTLuPDN7I3yc17qRNxQhl/9JyeN1M9uYMi4x28XM/mRmH5jZikbGm5ndFOa5zMyOTBmXmG0SIY+vhPEvN7PnzWxwyri3w+FLzWxh60WdWYRcRprZppTP0FUp47J+LltbhFwuS8ljRfjd6B6OS8x2MbMDzWxe+Fv7qpldkmGa1v+ueHjpxj3lQdDhbDVwCFAMvAL0T5vm28D08PlE4O7wef9w+vZAn3A5hQnPZRTQKXx+YW0u4evytt4ezcxlEvC7DPN2B94K/+4dPt87ybmkTX8RQSfIJG6XzwFHAisaGX8S8ChgwDHAiwndJk3lcWxtfMAXavMIX78N7NvW26IZuYwEHs4wvFmfyyTkkjbtKcDTSdwuwAHAkeHzLsDrGX6/Wv27sie2oOsuIeruO4DaS4im+iJwe/j8XuBEM7Nw+Cx3r3D3fwFvhstrK03m4u7z3H1r+HIBwfnkSRRluzRmLPCku3/k7h8DTwLjYooziubmcg5wV6tE1kzu/ncg21XBvgjc4YEFQDczO4CEbZOm8nD358M4IdnfkyjbpDG78x2LRTNzSfL35H13Xxw+LwNeI7gaZqpW/67siQU60yVE09/IumncvQrYBOwTcd7W1Nx4vkHwH1ytDma20MwWmNmEGOJrjqi5nBHuHrrXzGovZLPHbpfwkEMf4OmUwUnaLk1pLNekbZPmSP+eOPCEmS2y4NLCe4LhZvaKmT1qZgPCYXvsNjGzTgRF676UwYncLhYcEh0KvJg2qtW/K21+qU+Jxsy+ChwFHJ8y+GB3X2tmhwBPm9lyd1/dNhFG8hBwl7tXmNm3CPZynNDGMe2uicC97p56h5Y9bbvkDDMbRVCgP5sy+LPh9tgPeNLM/hm2/JJqMcFnqNzMTgLmAH3bNqTddgrwnLuntrYTt13MrITgn4jvufvmtowF9swWdJRLiNZNY2ZFwF7AhojztqZI8ZjZaOAnwKnuXlE73N3Xhn/fAuYT/NfXVprMxd03pMR/K1Aadd5W1px4JpK22y5h26UpjeWatG3SJDMbRPC5+qK7b6gdnrI9PgDup20PazXJ3Te7e3n4/G9AOzPblz1wm6TI9j1JxHYxs3YExXmmu8/OMEnrf1fa+uB8cx8Erf63CHYr1naUGJA2zXeo30nsr+HzAdTvJPYWbdtJLEouQwk6hvRNG7430D58vi/wBm3YYSRiLgekPD8NWBA+7w78K8xp7/B59yTnEk7Xj6CjiyV1u4Rx9KbxDkknU7/jy0tJ3CYR8jiIoE/JsWnDOwNdUp4/D4xryzwi5PKJ2s8UQdF6N9w+kT6XScolHL8XwXHqzkndLuH7ewfwv1mmafXvSptu2N14M08i6GW3GvhJOOwXBC1MgA7APeEX9iXgkJR5fxLOtwr4wh6Qy1zgP8DS8PFgOPxYYHn4JV0OfGMPyOWXwKthzPOAfinz/le4vd4Evp70XMLXVwPXp82XqO1C0Gp5H6gkODb2DeAC4IJwvAE3h3kuB45K4jaJkMetwMcp35OF4fBDwm3xSvjZ+0kCPltN5fLdlO/JAlL+6cj0uUxyLuE0kwg656bOl6jtQnBIxIFlKZ+hk9r6u6JLfYqIiCTQnngMWkREJOepQIuIiCSQCrSIiEgCqUCLiIgkkAq0iIhIAqlAi+whzGyflDsD/dvM1obPy83s9zGud6SZHRvX8kUkM13qU2QP4cHVsYYAmNnVBHfNmtoKqx4JlBNcTEJEWola0CJ7uLCF+3D4/Gozu93MnjWzd8zsdDP7dXjf3cfCyxliZqVm9kx4o4LHw7vyYGYX2877j88KbxxwAXBp2Fo/zsx6mNl9ZvZy+BiRsu4/m9kL4X1xz2+jt0QkJ6gFLZJ7PkVwH/H+wAvAGe5+uZndD5xsZo8A0wiuWf2hmZ0N/DfB1ZB+DPTx4IYm3dx9o5lNJ6W1bmZ3Av/j7v8ws4OAx4HDw3UPIrgMYmdgiZk94u7rWi1zkRyiAi2Sex5190ozWw4UAo+Fw5cTXDf5MGAgwR2ECKd5P5xmGTDTzOYQ3EUpk9FA/3BegK7hXYAAHnD3bcA2M5tHcC3pxpYjIlmoQIvkngoAd68xs0rfeT3fGoLvvAGvuvvwDPOeDHyO4PaAPzGzIzJMUwAc4+7bUweGBTv92sG6lrDILtIxaJH8swroYWbDIbjNnpkNMLMC4EB3nwf8iOAuRCVAGdAlZf4ngItqX5jZkJRxXzSzDma2D0HnspfjTEQkl6lAi+QZd98BnAn8ysxeIbhzz7EEu7r/Eu4aXwLc5O4bgYeA02o7iQEXA0eFHclWEnQiq7WM4E5lC4BrdPxZZNfpblYi0iJa+dQvkZynFrSIiEgCqQUtIiKSQGpBi4iIJJAKtIiISAKpQIuIiCSQCrSIiEgCqUCLiIgk0P8HJTwv26Lr3DcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Step 3: Compute relative change in total energy\n",
    "energy_relative_change = (total_energy_all - total_energy_all[0]) / total_energy_all[0] * 100\n",
    "\n",
    "# Step 4: Plot relative change in total energy over time\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.set_ylim(0, 5)\n",
    "ax.plot(range(ntime), energy_relative_change, marker='o', label='Relative Change (%)')\n",
    "ax.set_title('Relative Change in Total Energy')\n",
    "ax.set_xlabel('Timestep')\n",
    "ax.set_ylabel('Percentage Change (%)')\n",
    "ax.axhline(0, color='k', linestyle='--', linewidth=1)  # Reference line\n",
    "ax.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "704e513c-6112-4f8a-b54b-6caa4043934d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ORC)",
   "language": "python",
   "name": "sys_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.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}