{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Correlation functions\n\nGenerate the velocity for a Ornstein-Uhlenbeck process and compute its\nautocorrelation function.\n\nWe also display the force velocity correlation as an example of using the\nroutine `correlation`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport tidynamics\nimport matplotlib.pyplot as plt\n\nplt.rcParams['figure.figsize'] = 5.12, 3.84\nplt.rcParams['figure.subplot.bottom'] = 0.18\nplt.rcParams['figure.subplot.left'] = 0.16\n\n# Generate data for a Ornstein-Uhlenbeck process\n\ngamma = 2.7\nT = 0.1\ndt = 0.02\nv_factor = np.sqrt(2*T*gamma*dt)\n\nN = 32768\nv = 0\nfor i in range(100):\n    noise_force = v_factor*np.random.normal()\n    v = v - gamma*v*dt + noise_force\n\nv_data = []\nnoise_data = []\n\nfor i in range(N):\n    noise_force = v_factor*np.random.normal()\n    v = v - gamma*v*dt + noise_force\n    v_data.append(v)\n    noise_data.append(noise_force)\nv_data = np.array(v_data)\nnoise_data = np.array(noise_data)/np.sqrt(dt)\n\n# Compute the autocorrelation function and plot the result\n\nacf = tidynamics.acf(v_data)[:N//64]\n\ntime = np.arange(N//64)*dt\n\nplt.plot(time, acf, label='VACF (num.)')\n\nplt.plot(time, T*np.exp(-gamma*time), label='VACF (theo.)')\n\nplt.legend()\nplt.title('Velocity autocorrelation')\nplt.xlabel(r'$\\tau$')\nplt.ylabel(r'$\\langle v(t) v(t+\\tau) \\rangle$')\n\n# Compute the force velocity correlation and plot the result\n\nplt.figure()\n\ntime = np.arange(N)*dt\ntwotimes = np.concatenate((-time[1:][::-1], time))\n\nplt.plot(twotimes, tidynamics.correlation(noise_data, v_data),\n         label='num.')\nplt.plot(time, 2*T/gamma*np.exp(-gamma*time), label='theo.')\n\nplt.xlim(-5/gamma, 5/gamma)\n\nplt.ylim(-2*T/gamma/10, 1.1*2*T/gamma)\n\nplt.legend()\nplt.title('Force-velocity correlation')\nplt.xlabel(r'$\\tau$')\nplt.ylabel(r'$\\langle F(t) v(t+\\tau) \\rangle$')\n\nplt.show()"
      ]
    }
  ],
  "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.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}