Synthetic_field.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ad8a883e",
   "metadata": {},
   "source": [
    "# Synthetic turbulent magnetic field with CRPropa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dac55f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from crpropa import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "267eeb04",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pylab import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21c59b7b",
   "metadata": {},
   "source": [
    "## 1. Creation of field"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f8697ba",
   "metadata": {},
   "source": [
    "### Properties of grid and field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "248df946",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grid\n",
    "spacing = 10 * pc\n",
    "grid_number = 100\n",
    "\n",
    "# Field\n",
    "Brms = 5 * nG\n",
    "lmin = 50 * pc \n",
    "lmax = 1000 * pc \n",
    "index = 5./3. # Kolmogorov index = 5/3, Kraichnan index = 3/2\n",
    "random_seed = 42\n",
    "\n",
    "# Creation of turbulent field\n",
    "spectrum = SimpleTurbulenceSpectrum(Brms, lmin, lmax, index)\n",
    "grid = GridProperties(Vector3d(0), grid_number, spacing)\n",
    "BField = SimpleGridTurbulence(spectrum, grid, random_seed)\n",
    "\n",
    "# Save field\n",
    "dumpGridToTxt(BField.getGrid(), 'kolmogorov_field.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3043ed73",
   "metadata": {},
   "source": [
    "### Check of field properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c0b75ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Lc = {:.1f} kpc'.format(BField.getCorrelationLength() / pc))  \n",
    "print('sqrt(<B^2>) = {:.1f} nG'.format(BField.getBrms() / nG))   \n",
    "print('<|B|> = {:.1f} nG'.format(BField.getMeanFieldStrength() / nG))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27169e7d",
   "metadata": {},
   "source": [
    "## 2. Diagnostic of field"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a4ea417",
   "metadata": {},
   "source": [
    "### Extraction of field components of grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "951275a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "field = np.genfromtxt('kolmogorov_field.txt')\n",
    "\n",
    "# Field components\n",
    "Bx = field[:,0] / nG\n",
    "By = field[:,1] / nG\n",
    "Bz = field[:,2] / nG\n",
    "\n",
    "# Field magntiudes\n",
    "B = np.sqrt(Bx**2 + By**2 + Bz**2)\n",
    "\n",
    "# Squared values of field\n",
    "Bx_sqr = Bx**2\n",
    "By_sqr = By**2\n",
    "Bz_sqr = Bz**2\n",
    "B_sqr = B**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "361f402d",
   "metadata": {},
   "source": [
    "### Check of mean values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85638580",
   "metadata": {},
   "outputs": [],
   "source": [
    "Bx_mean = np.mean(Bx)\n",
    "By_mean = np.mean(By)\n",
    "Bz_mean = np.mean(Bz)\n",
    "\n",
    "print('Mean Bx', Bx_mean)\n",
    "print('Mean By', By_mean)\n",
    "print('Mean Bz', Bz_mean)\n",
    "print('Mean B', B_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6480500",
   "metadata": {},
   "source": [
    "### Plotting distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "861964c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot of field components\n",
    "plt.figure(figsize = (15, 10))\n",
    "#\n",
    "plt.subplot(2, 2, 1)\n",
    "counts, bins = np.histogram(Bx, bins=50)\n",
    "plt.hist(bins[:-1], bins, weights=counts)\n",
    "plt.xlabel(r'$\\delta B_x$ [nG]', fontsize=17)\n",
    "plt.grid()\n",
    "#\n",
    "plt.subplot(2, 2, 2)\n",
    "counts, bins = np.histogram(By, bins=50)\n",
    "plt.hist(bins[:-1], bins, weights=counts)\n",
    "plt.xlabel(r'$\\delta B_y$ [nG]', fontsize=17)\n",
    "plt.grid()\n",
    "#\n",
    "plt.subplot(2, 2, 3)\n",
    "counts, bins = np.histogram(Bz, bins=50)\n",
    "plt.hist(bins[:-1], bins, weights=counts)\n",
    "plt.xlabel(r'$\\delta B_z$ [nG]', fontsize=17)\n",
    "plt.grid()\n",
    "#\n",
    "plt.subplot(2, 2, 4)\n",
    "counts, bins = np.histogram( B, bins=50)\n",
    "plt.hist(bins[:-1], bins, weights=counts)\n",
    "plt.xlabel(r'$\\delta B [nG]$', fontsize=17)\n",
    "plt.grid()\n",
    "#\n",
    "plt.savefig('mean_field_components_kolmogorov.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa24031b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot of squared field components\n",
    "plt.figure(figsize = (15, 10))\n",
    "#\n",
    "plt.subplot(2, 2, 1)\n",
    "counts, bins = np.histogram(Bx_sqr, bins=50)\n",
    "plt.hist(bins[:-1], bins, weights=counts)\n",
    "plt.xlabel(r'$(\\delta B_x)^2$', fontsize=17)\n",
    "plt.yscale('log')\n",
    "plt.grid()\n",
    "#\n",
    "plt.subplot(2, 2, 2)\n",
    "counts, bins = np.histogram(By_sqr, bins=50)\n",
    "plt.hist(bins[:-1], bins, weights=counts)\n",
    "plt.xlabel(r'$(\\delta B_y)^2$', fontsize=17)\n",
    "plt.yscale('log')\n",
    "plt.grid()\n",
    "#\n",
    "plt.subplot(2, 2, 3)\n",
    "counts, bins = np.histogram(Bz_sqr, bins=50)\n",
    "plt.hist(bins[:-1], bins, weights=counts)\n",
    "plt.xlabel(r'$(\\delta B_z)^2$', fontsize=17)\n",
    "plt.yscale('log')\n",
    "plt.grid()\n",
    "#\n",
    "plt.subplot(2, 2, 4)\n",
    "counts, bins = np.histogram( B**2, bins=50)\n",
    "plt.hist(bins[:-1], bins, weights=counts)\n",
    "plt.xlabel(r'$(\\delta B)^2$', fontsize=17)\n",
    "plt.yscale('log')\n",
    "plt.grid()\n",
    "#\n",
    "plt.savefig('squared_field_components_kolmogorov.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba0e0c53",
   "metadata": {},
   "source": [
    "## 3. Power spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f772e403",
   "metadata": {},
   "source": [
    "### FFT of field components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4db2f6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3D FFT to go from x-space to k-space\n",
    "B_ifft = np.array((Bx, By, Bz))\n",
    "Bk = np.fft.fftn(B_ifft)\n",
    "Bkx = Bk[0]\n",
    "Bky = Bk[1]\n",
    "Bkz = Bk[2]\n",
    "Bk_total = np.sqrt( (Bkx.real**2 + Bkx.imag**2) + (Bky.real**2 + Bky.imag**2) + (Bkz.real**2 + Bkz.imag**2) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe6d5c73",
   "metadata": {},
   "source": [
    "### Computation of k-values of $\\delta B^2 (k)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e66e9bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generating N_grids-possible discrete wave numbers\n",
    "K = np.zeros(grid_number)\n",
    "for i in range(0, grid_number):\n",
    "    K[i] = (i / grid_number) - (i / (grid_number / 2.) )\n",
    "\n",
    "# Computing vector with corresponding k values to each grid point\n",
    "grid_number_z = int(np.floor(grid_number/2) + 1)\n",
    "k = np.zeros(grid_number * grid_number * grid_number_z)\n",
    "k_max = (spacing / pc) / (lmin / pc)\n",
    "k_min = (spacing / pc) / (lmax / pc)\n",
    "for ix in range(0, grid_number):\n",
    "    for iy in range(0, grid_number):\n",
    "        for iz in range(0, grid_number_z):\n",
    "            i = ix * grid_number * grid_number_z + iy * grid_number_z + iz\n",
    "            k[i] = np.sqrt(K[ix]**2 + K[iy]**2 + K[iz]**2)\n",
    "\n",
    "range_k = len(k)\n",
    "Bk_total = Bk_total[0:range_k]\n",
    "\n",
    "# Binning of k-values and respective B^2\n",
    "counts, bins = np.histogram(k, bins=50)\n",
    "k_values = np.zeros(len(bins) - 1)\n",
    "Bk_sqrd_values = np.zeros(len(bins) - 1)\n",
    "for j in range(0, len(bins) - 1):\n",
    "    k_values[j] = (bins[j] + bins[j+1]) / 2\n",
    "    index_where = np.where(np.logical_and(k >= bins[j], k <= bins[j+1]))\n",
    "    Bk_sqrd_values[j] = np.mean(Bk_total[index_where]**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d766cb71",
   "metadata": {},
   "source": [
    "### Plotting spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "767b7fad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting power spectrum\n",
    "plt.figure(figsize = (7, 5))\n",
    "plt.plot(k_values, Bk_sqrd_values, color='blue', label=r'Test')\n",
    "plt.loglog()\n",
    "plt.ylim(1e9, 1e13)\n",
    "plt.xlim(k_min, k_max)\n",
    "plt.xlabel(r'k', fontsize=17)\n",
    "plt.ylabel(r'$\\delta B^2(k)$', fontsize=17)\n",
    "plt.grid()\n",
    "#plt.legend(fontsize=17)\n",
    "plt.savefig('Plots/power_spectrum_kolmogorov.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a95dc21",
   "metadata": {},
   "source": [
    "### Corroborating spectrum slope"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9901d6eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check slope\n",
    "index_slope = np.where(np.logical_and(k_values >= k_min, k_values <= k_max))\n",
    "linear = linregress(np.log10(k_values[index_slope]), np.log10(Bk_sqrd_values[index_slope]))\n",
    "print('Numerical slope of turbulence spectrum:', linear[0])\n",
    "print('Theoretical slope of turbulence spectrum:', -index-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "138c3453",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}