diff --git a/GRAMSES_tutorial.ipynb b/GRAMSES_tutorial.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..868de4362fcf7af3090553d6afc5b93d47f5a59d
--- /dev/null
+++ b/GRAMSES_tutorial.ipynb
@@ -0,0 +1,724 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "------------------------------------------\n",
+ "# GRAMSES Tutorial\n",
+ "------------------------------------------\n",
+ "\n",
+ "## Introduction\n",
+ "\n",
+ "This Jupyter Notebook will (try to) guide you step-by-step through some GRAMSES examples.\n",
+ "\n",
+ "\n",
+ "## Definitions and conventions\n",
+ "\n",
+ "### Metric\n",
+ "\n",
+ "In GRAMSES the metric is given in the 3+1 form:\n",
+ "\\begin{equation}\n",
+ "{\\rm d}s^2 = g_{\\mu\\nu}{\\rm d}x^{\\mu}{\\rm d}x^{\\nu} = -\\alpha^2 {\\rm d}t^2+\\gamma_{ij}\\left(\\beta^i{\\rm d}t+{\\rm d}x^i\\right)\\left(\\beta^j{\\rm d}t+{\\rm d}x^j\\right)\n",
+ "\\end{equation}\n",
+ "\n",
+ "where $\\alpha$ is the lapse function, $\\beta^i$ the shift vector and $\\gamma_{ij}$ the induced metric on the spatial hypersurfaces, which in the constrained formulation adopted by GRAMSES is approximated by a conformally-flat metric, i.e.\n",
+ "\n",
+ "\\begin{equation}\n",
+ "\\gamma_{ij}=\\psi^4\\delta_{ij}\n",
+ "\\end{equation}\n",
+ "with $\\psi$ being the conformal factor and $\\delta_{ij}$ the Kronecker delta. \n",
+ "\n",
+ "### Matter sources\n",
+ "\n",
+ "The usual matter source terms defined in the 3+1 formalism are given by the following projections of the energy-momentum tensor:\n",
+ "\n",
+ "\\begin{align}\n",
+ "\\rho &\\equiv n_\\mu n_\\nu T^{\\mu\\nu}\\,,\\label{rho-def}\\\\\n",
+ "S_i&\\equiv-\\gamma_{i\\mu}n_\\nu T^{\\mu\\nu}\\,,\\label{S_i-def}\\\\\n",
+ "S_{ij}&\\equiv\\gamma_{i\\mu}\\gamma_{j\\nu}T^{\\mu\\nu}\\,, \\qquad\\qquad S=\\gamma^{ij}S_{ij}\\,.\n",
+ "\\end{align}\n",
+ "\n",
+ "\n",
+ "GRAMSES uses and outputs the following conformal matter source terms:\n",
+ "\n",
+ "\\begin{align}\n",
+ "s_0({\\bf x})&\\equiv\\sqrt{\\gamma}\\rho&&\\propto{m\\alpha{u}^0}\\,,\\label{eq:s0-cic}\\\\\n",
+ "s_i({\\bf x})&\\equiv\\sqrt{\\gamma}S_i&&\\propto m{u}_i\\,,\\\\\n",
+ "s_{ij}({\\bf x})&\\equiv\\sqrt{\\gamma}S_{ij}&&\\propto m\\frac{{u}_i{u}_j}{\\alpha{u}^0}\\,.\\label{eq:s-cic}\n",
+ "\\end{align}\n",
+ "In these, ${\\bf x}$ is a (discrete) position vector on the cartesian simulation grid and the proportionality symbol in each equation stands for the standard cloud-in-cell (CIC) weights used for the particle-mesh projection. From these we have the following useful relations:\n",
+ "\\begin{align}\n",
+ " s_0&=\\rho\\Gamma\\,,\\\\\n",
+ " s_i&=\\frac{\\rho}{\\Gamma}u_i\\,,\\label{eq:s_i-rel}\\\\\n",
+ " s_{ij}&=\\frac{\\rho}{\\Gamma^2}u_iu_j\\,, \\qquad\\implies s=\\rho(1-\\Gamma^{-2})\\,,\\\\\n",
+ " u_i&=\\Gamma^2\\frac{s_i}{s_0}\\,,\n",
+ "\\end{align}\n",
+ "where $\\Gamma\\equiv\\alpha{u}^0=\\sqrt{1+\\gamma^{ij}u_iu_j}$ is the Lorentz factor."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 0. Converting Fortran data to numpy format\n",
+ "\n",
+ "The first step is to run the readgrav script provided."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# You can run it from here where the arguments are: \n",
+ "#\n",
+ "# output_000xx/ ncpus ilevelmin\n",
+ "#\n",
+ "# where 2^ilevelmin is the cells number along 1D.\n",
+ "\n",
+ "%run -i 'readgrav_gr.py' output_00004/ 4 7"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Import python libraries needed\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import math\n",
+ "from matplotlib import gridspec\n",
+ "from matplotlib.colors import LogNorm\n",
+ "import matplotlib.ticker as ticker\n",
+ "\n",
+ "from nbodykit.lab import *\n",
+ "\n",
+ "# Plots style\n",
+ "font = {'size':14}\n",
+ "plt.rc('font', **font)\n",
+ "plt.rc('font', family='serif')\n",
+ "plt.rc('text', usetex=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Read snapshot in numpy format"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Simulations & snapshot info\n",
+ "a_out = 1.0\n",
+ "box_size = 256. # Mpc/h\n",
+ "grid_size= 128\n",
+ "\n",
+ "# load snapshot file in npy format\n",
+ "#grav_raw = np.load('output_00004_B256_PM128/grav_00004_new.out.npy')\n",
+ "grav_raw = np.load('output_00004/grav_00004.out.npy')\n",
+ "\n",
+ "grad_b = [grav_raw[:,1],grav_raw[:,2],grav_raw[:,3]]\n",
+ "gr_pot = [grav_raw[:,4],grav_raw[:,5],grav_raw[:,6],grav_raw[:,7],grav_raw[:,8],grav_raw[:,9]]\n",
+ "gr_mat = [grav_raw[:,10],grav_raw[:,11],grav_raw[:,12],grav_raw[:,13],grav_raw[:,14]]\n",
+ "x = grav_raw[:,15]\n",
+ "y = grav_raw[:,16]\n",
+ "z = grav_raw[:,17]\n",
+ "\n",
+ "# Conversion from internal code units to physical units -- do not modify!\n",
+ "ctilde = 3*10**3/box_size\n",
+ "norm_B = 1/a_out**2/ctilde\n",
+ "norm_si= 2/(ctilde*box_size)**2/a_out**4\n",
+ "\n",
+ "B_x = gr_pot[2] * norm_B\n",
+ "B_y = gr_pot[3] * norm_B\n",
+ "B_z = gr_pot[4] * norm_B\n",
+ "\n",
+ "dx_b=grad_b[0] * norm_B\n",
+ "dy_b=grad_b[1] * norm_B\n",
+ "dz_b=grad_b[2] * norm_B\n",
+ "\n",
+ "s_x = gr_mat[1] * norm_si\n",
+ "s_y = gr_mat[2] * norm_si\n",
+ "s_z = gr_mat[3] * norm_si"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Construct 3D box of fields\n",
+ "xx_box= grid_size*x\n",
+ "yy_box= grid_size*y\n",
+ "zz_box= grid_size*z\n",
+ "x_box = xx_box.astype(int)\n",
+ "y_box = yy_box.astype(int)\n",
+ "z_box = zz_box.astype(int)\n",
+ "\n",
+ "# Shift vector components\n",
+ "B_x_box = np.zeros((grid_size,grid_size,grid_size))\n",
+ "B_y_box = np.zeros((grid_size,grid_size,grid_size))\n",
+ "B_z_box = np.zeros((grid_size,grid_size,grid_size))\n",
+ "B_x_box[(x_box,y_box,z_box)] = B_x\n",
+ "B_y_box[(x_box,y_box,z_box)] = B_y\n",
+ "B_z_box[(x_box,y_box,z_box)] = B_z\n",
+ "\n",
+ "dx_b_box = np.zeros((grid_size,grid_size,grid_size))\n",
+ "dy_b_box = np.zeros((grid_size,grid_size,grid_size))\n",
+ "dz_b_box = np.zeros((grid_size,grid_size,grid_size))\n",
+ "dx_b_box[(x_box,y_box,z_box)] = dx_b\n",
+ "dy_b_box[(x_box,y_box,z_box)] = dy_b\n",
+ "dz_b_box[(x_box,y_box,z_box)] = dz_b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Matter source terms\n",
+ "s0_box = np.zeros((grid_size,grid_size,grid_size))\n",
+ "s_x_box= np.zeros((grid_size,grid_size,grid_size))\n",
+ "s_y_box= np.zeros((grid_size,grid_size,grid_size))\n",
+ "s_z_box= np.zeros((grid_size,grid_size,grid_size))\n",
+ "\n",
+ "s0_box [(x_box,y_box,z_box)]= gr_mat[0]\n",
+ "s_x_box[(x_box,y_box,z_box)]= s_x\n",
+ "s_y_box[(x_box,y_box,z_box)]= s_y\n",
+ "s_z_box[(x_box,y_box,z_box)]= s_z\n",
+ "\n",
+ "# Magnitude of vector field\n",
+ "s_i_box = np.sqrt(s_x_box**2+s_y_box**2+s_z_box**2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 1. Plot density and momentum maps"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Select a slice of the data\n",
+ "s0_slice = s0_box[:,:,50]\n",
+ "s0_slice+= 1E-10 # Add a small number to visualise in log scale.\n",
+ "\n",
+ "s_i_slice = s_i_box[:,:,50]\n",
+ "s_i_slice+= 1E-10 # Add a small number to visualise in log scale.\n",
+ "\n",
+ "#Figs\n",
+ "fig= plt.figure(figsize=(10, 6)) \n",
+ "gs = gridspec.GridSpec(1,2, wspace=0.5)\n",
+ "\n",
+ "ax0 = plt.subplot(gs[0,0])\n",
+ "ax1 = plt.subplot(gs[0,1])\n",
+ "\n",
+ "fig_s0 = ax0.imshow(s0_slice , origin='lower', cmap='inferno', aspect='equal',\n",
+ " norm=LogNorm(vmin=2E-1,vmax=1*s0_slice.max()),interpolation='None')\n",
+ "\n",
+ "fig_si = ax1.imshow(s_i_slice , origin='lower', cmap='jet', aspect='equal',\n",
+ " norm=LogNorm(vmin=1E-10,vmax=1*s_i_slice.max()),interpolation='None')\n",
+ "\n",
+ "#Colorbars\n",
+ "cx1 = fig.add_axes([0.45,0.25,0.02,.5])\n",
+ "cx2 = fig.add_axes([0.92,0.25,0.02,.5])\n",
+ "cb_s0 = plt.colorbar(fig_s0, cax = cx1, orientation='vertical')\n",
+ "cb_si = plt.colorbar(fig_si, cax = cx2, orientation='vertical')\n",
+ "cb_s0.ax.tick_params(labelsize=14, color='k')\n",
+ "cb_si.ax.tick_params(labelsize=14, color='k')\n",
+ "\n",
+ "plt.show()\n",
+ "#fig.savefig(\"2D_map_s0_z-1_dpi_200.pdf\",bbox_inches='tight',dpi=200)\n",
+ "plt.close()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 1.1 Measure $P(k)$ using N-body kit"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Select field to measure\n",
+ "\n",
+ "k_fun = 2*np.pi/box_size # fundamental mode of the box\n",
+ "k_Nyq = np.pi*grid_size/box_size # Nyquist mode\n",
+ "\n",
+ "sim_field = s0_box\n",
+ "\n",
+ "mesh_var = ArrayMesh( sim_field, Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_var = FFTPower( mesh_var, mode='1d', kmin= k_fun )\n",
+ "k_bins= np.array( r_var.power['k'] )\n",
+ "Pk_var= np.array( r_var.power['power'].real )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ " # CAMB z=0 matter P(k)\n",
+ "Pk_lin_z_0 = np.loadtxt('cmc_gauge_matterpower_z_0.dat',usecols=[1],skiprows=0,unpack=True,dtype=np.float32)[0:]\n",
+ "k_lin = np.loadtxt('cmc_gauge_matterpower_z_0.dat',usecols=[0],skiprows=0,unpack=True,dtype=np.float32)[0:]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plots\n",
+ "font = {'size':14}\n",
+ "plt.rc('font', **font)\n",
+ "plt.rc('font', family='serif')\n",
+ "plt.rc('text', usetex=True)\n",
+ "fig = plt.figure(figsize=(4, 5)) \n",
+ "gs = gridspec.GridSpec(2,1, height_ratios=[4,1]) \n",
+ "\n",
+ "ax0 = plt.subplot(gs[0])\n",
+ "fig_lin= ax0.plot(k_lin,Pk_lin_z_0,'gray' , linewidth=1.0,zorder=-10, label='linear')\n",
+ "fig_Pk = ax0.scatter( k_bins, Pk_var, edgecolors='b', s=20,marker='.',facecolors='b',zorder=-2 )\n",
+ "plt.tick_params(axis='y', direction='in', which='both', labelleft='on', labelright='off', right='on')\n",
+ "plt.tick_params(axis='x', direction='in', which='both', labeltop='off', labelbottom='on', top='on')\n",
+ "ax0.set_ylabel(r'$P(k)$')\n",
+ "ax0.set_xlabel(r'$k$ [$h$ Mpc$^{-1}$]')\n",
+ "ax0.set_xscale('log')\n",
+ "ax0.set_yscale('log')\n",
+ "ax0.set_xlim( [ k_fun, 1.5*k_Nyq ] )\n",
+ "plt.fill_betweenx(np.array([1E-20,1E20]),k_Nyq,20 , alpha=0.2,color='red',interpolate=True)\n",
+ "ax0.set_ylim([1E1, 4E4])\n",
+ "\n",
+ "plt.show()\n",
+ "#fig.savefig(\"Pk_gramses.pdf\", bbox_inches='tight')\n",
+ "plt.close()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 2. Vector potential\n",
+ "\n",
+ "The vector mode of the shift vector is given by the linear combination below:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Bv_x = B_x_box -4.*dx_b_box\n",
+ "Bv_y = B_y_box -4.*dy_b_box\n",
+ "Bv_z = B_z_box -4.*dz_b_box\n",
+ "\n",
+ "# Magnitude of vector field\n",
+ "Bv_box = np.sqrt(Bv_x**2+Bv_y**2+Bv_z**2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Select a slice of the data\n",
+ "s0_slice = s0_box[:,:,50]\n",
+ "s0_slice+= 1E-10 # Add a small number to visualise in log scale.\n",
+ "\n",
+ "Bv_slice = Bv_box [:,:,50]\n",
+ "\n",
+ "#Figs\n",
+ "fig= plt.figure(figsize=(10, 6)) \n",
+ "gs = gridspec.GridSpec(1,2, wspace=0.5) \n",
+ "\n",
+ "ax0 = plt.subplot(gs[0,0])\n",
+ "ax1 = plt.subplot(gs[0,1])\n",
+ "\n",
+ "fig_s0 = ax0.imshow(s0_slice , origin='lower', cmap='inferno', aspect='equal',\n",
+ " norm=LogNorm(vmin=2E-1,vmax=1*s0_slice.max()),interpolation='None')\n",
+ "\n",
+ "fig_Bv = ax1.imshow(Bv_slice , origin='lower', cmap='jet', aspect='equal',\n",
+ " norm=LogNorm(vmin=1E-10,vmax=1*Bv_slice.max()),interpolation='bicubic')\n",
+ "\n",
+ "#Colorbars\n",
+ "cx1 = fig.add_axes([0.45,0.25,0.02,.5])\n",
+ "cx2 = fig.add_axes([0.92,0.25,0.02,.5])\n",
+ "cb_s0 = plt.colorbar(fig_s0, cax = cx1, orientation='vertical')\n",
+ "cb_Bv = plt.colorbar(fig_Bv, cax = cx2, orientation='vertical')\n",
+ "cb_s0.ax.tick_params(labelsize=14, color='k')\n",
+ "cb_Bv.ax.tick_params(labelsize=14, color='k')\n",
+ "\n",
+ "plt.show()\n",
+ "#fig.savefig(\"2D_map_s0_z-1_dpi_200.pdf\",bbox_inches='tight',dpi=200)\n",
+ "plt.close()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2.1 Power spectrum"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Measure P(k) from Bv\n",
+ "\n",
+ "mesh_B_x = ArrayMesh( Bv_x, Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_B_x_Pk = FFTPower( mesh_B_x, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_B_x = np.array( r_B_x_Pk.power['k'] )\n",
+ "Pk_B_x= np.array( r_B_x_Pk.power['power'].real )\n",
+ "\n",
+ "mesh_B_y = ArrayMesh( Bv_y, Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_B_y_Pk = FFTPower( mesh_B_y, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_B_y = np.array( r_B_y_Pk.power['k'] )\n",
+ "Pk_B_y= np.array( r_B_y_Pk.power['power'].real )\n",
+ "\n",
+ "mesh_B_z = ArrayMesh( Bv_z, Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_B_z_Pk = FFTPower( mesh_B_z, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_B_z = np.array( r_B_z_Pk.power['k'] )\n",
+ "Pk_B_z= np.array( r_B_z_Pk.power['power'].real )\n",
+ "\n",
+ "Pk_Bv = Pk_B_x+Pk_B_y+Pk_B_z"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2.2 Calculate Curl: $[\\nabla\\times V]_i=\\epsilon_{ijk}\\partial_jV_k$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def curl_array(array_x, array_y, array_z, grid_len, box_len):\n",
+ " Vx1_z = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vx1_y = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vy1_z = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vy1_x = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vz1_y = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vz1_x = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vx2_z = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vx2_y = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vy2_z = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vy2_x = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vz2_y = np.zeros((grid_size,grid_size,grid_size))\n",
+ " Vz2_x = np.zeros((grid_size,grid_size,grid_size))\n",
+ " \n",
+ " # Shifted positions & PBC\n",
+ " xbox1 = x_box + 1\n",
+ " xbox2 = x_box - 1\n",
+ " ybox1 = y_box + 1\n",
+ " ybox2 = y_box - 1\n",
+ " zbox1 = z_box + 1\n",
+ " zbox2 = z_box - 1\n",
+ " for i in (range(len(xbox1))):\n",
+ " if(xbox1[i]==grid_len): xbox1[i]=0\n",
+ " if(xbox2[i]==-1): xbox2[i]=grid_len-1\n",
+ " if(ybox1[i]==grid_len): ybox1[i]=0\n",
+ " if(ybox2[i]==-1): ybox2[i]=grid_len-1\n",
+ " if(zbox1[i]==grid_len): zbox1[i]=0\n",
+ " if(zbox2[i]==-1): zbox2[i]=grid_len-1\n",
+ " \n",
+ " # Define nodes for FD:\n",
+ " Vx1_z[(x_box,y_box,zbox1)] = array_x\n",
+ " Vx1_y[(x_box,ybox1,z_box)] = array_x\n",
+ " Vy1_z[(x_box,y_box,zbox1)] = array_y\n",
+ " Vy1_x[(xbox1,y_box,z_box)] = array_y\n",
+ " Vz1_y[(x_box,ybox1,z_box)] = array_z\n",
+ " Vz1_x[(xbox1,y_box,z_box)] = array_z\n",
+ " \n",
+ " Vx2_z[(x_box,y_box,zbox2)] = array_x\n",
+ " Vx2_y[(x_box,ybox2,z_box)] = array_x\n",
+ " Vy2_z[(x_box,y_box,zbox2)] = array_y\n",
+ " Vy2_x[(xbox2,y_box,z_box)] = array_y\n",
+ " Vz2_y[(x_box,ybox2,z_box)] = array_z\n",
+ " Vz2_x[(xbox2,y_box,z_box)] = array_z\n",
+ " \n",
+ " # Calculate FD using shifted arrays\n",
+ " dx = box_len/float(grid_size)\n",
+ " curl_V_x= ((Vz1_y-Vz2_y)-(Vy1_z-Vy2_z))/(2*dx)\n",
+ " curl_V_y= ((Vx1_z-Vx2_z)-(Vz1_x-Vz2_x))/(2*dx)\n",
+ " curl_V_z= ((Vy1_x-Vy2_x)-(Vx1_y-Vx2_y))/(2*dx)\n",
+ " return [curl_V_x,curl_V_y,curl_V_z]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "curl_B = curl_array(B_x, B_y, B_z, grid_size, box_size )\n",
+ "\n",
+ "# curl of momentum field s_i\n",
+ "curl_s = curl_array(s_x, s_y, s_z, grid_size, box_size )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Measure Power spectra of curls\n",
+ "\n",
+ "mesh_cB_x = ArrayMesh( curl_B[0], Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_cB_x_Pk = FFTPower( mesh_cB_x, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_cB_x = np.array( r_cB_x_Pk.power['k'] )\n",
+ "Pk_cB_x= np.array( r_cB_x_Pk.power['power'].real )\n",
+ "\n",
+ "mesh_cB_y = ArrayMesh( curl_B[1], Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_cB_y_Pk = FFTPower( mesh_cB_y, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_cB_y = np.array( r_cB_y_Pk.power['k'] )\n",
+ "Pk_cB_y= np.array( r_cB_y_Pk.power['power'].real )\n",
+ "\n",
+ "mesh_cB_z = ArrayMesh( curl_B[2], Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_cB_z_Pk = FFTPower( mesh_cB_z, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_cB_z = np.array( r_cB_z_Pk.power['k'] )\n",
+ "Pk_cB_z= np.array( r_cB_z_Pk.power['power'].real )\n",
+ "\n",
+ "# Spectrum from curl\n",
+ "Pk_cBB = Pk_cB_x+Pk_cB_y+Pk_cB_z\n",
+ "Pk_BB = Pk_cBB/k_cB_x**2\n",
+ "k_cB = k_cB_x\n",
+ "\n",
+ "# Curl of density field\n",
+ "mesh_cs_x = ArrayMesh( curl_s[0], Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_cs_x_Pk = FFTPower( mesh_cs_x, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_cs_x = np.array( r_cs_x_Pk.power['k'] )\n",
+ "Pk_cs_x= np.array( r_cs_x_Pk.power['power'].real )\n",
+ "\n",
+ "mesh_cs_y = ArrayMesh( curl_s[1], Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_cs_y_Pk = FFTPower( mesh_cs_y, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_cs_y = np.array( r_cs_y_Pk.power['k'] )\n",
+ "Pk_cs_y= np.array( r_cs_y_Pk.power['power'].real )\n",
+ "\n",
+ "mesh_cs_z = ArrayMesh( curl_s[2], Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_cs_z_Pk = FFTPower( mesh_cs_z, mode='1d', kmin=2*np.pi/box_size )\n",
+ "k_cs_z = np.array( r_cs_z_Pk.power['k'] )\n",
+ "Pk_cs_z= np.array( r_cs_z_Pk.power['power'].real )\n",
+ "\n",
+ "# Spectrum from curl\n",
+ "Pk_csi= Pk_cs_x+Pk_cs_y+Pk_cs_z\n",
+ "Pk_si = Pk_csi/k_cs_x**2\n",
+ "k_cs = k_cs_x"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Plot all the vector modes $P(k)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plots\n",
+ "font = {'size':14}\n",
+ "plt.rc('font', **font)\n",
+ "plt.rc('font', family='serif')\n",
+ "plt.rc('text', usetex=True)\n",
+ "fig = plt.figure(figsize=(4, 5)) \n",
+ "gs = gridspec.GridSpec(2,1, height_ratios=[4,1]) \n",
+ "\n",
+ "ax0 = plt.subplot(gs[0])\n",
+ "fig_Pk = ax0.scatter( k_B_x, k_B_x**3*Pk_Bv/(2*np.pi**2) \n",
+ " , edgecolors='b', s=20,marker='.',facecolors='b',zorder=-2 )\n",
+ "fig_Pk = ax0.scatter( k_cB , k_cB**3*Pk_BB*1/(2*np.pi**2)\n",
+ " , edgecolors='r', s=20,marker='.',facecolors='r',zorder=-2 )\n",
+ "fig_Pk = ax0.scatter( k_cs , k_cs**3*Pk_si/k_cs**4/(2*np.pi**2) \n",
+ " , edgecolors='gray', s=20,marker='.',facecolors='gray',zorder=-2 )\n",
+ "\n",
+ "plt.tick_params(axis='y', direction='in', which='both', labelleft='on', labelright='off', right='on')\n",
+ "plt.tick_params(axis='x', direction='in', which='both', labeltop='off', labelbottom='on', top='on')\n",
+ "ax0.set_ylabel(r'$k^3P(k)/(2\\pi^2)$')\n",
+ "ax0.set_xlabel(r'$k$ [$h$ Mpc$^{-1}$]')\n",
+ "ax0.set_xscale('log')\n",
+ "ax0.set_yscale('log')\n",
+ "ax0.set_xlim( [ k_fun, 1.5*k_Nyq ] )\n",
+ "ax0.set_ylim([1E-20, 1E-15])\n",
+ "plt.fill_betweenx(np.array([1E-20,1E-10]),k_Nyq, 20, alpha=0.2,color='red',interpolate=True)\n",
+ "\n",
+ "plt.show()\n",
+ "#fig.savefig(\"Pk_gramses.pdf\", bbox_inches='tight')\n",
+ "plt.close()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 3. Gravitational slip: $\\chi=\\vert \\Phi-\\Psi \\vert$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define nonlinear scalar perturbations\n",
+ "\n",
+ "# Conformal factor perturbation\n",
+ "Psi = gr_pot[0]/(a_out*ctilde)**2/2\n",
+ "\n",
+ "# Lapse function perturbation\n",
+ "lapse = 1. + gr_pot[1]/(a_out*ctilde)**2/(1. - Psi)\n",
+ "Phi = 0.5*(lapse**2-1.)\n",
+ "\n",
+ "# Define gravitational slip\n",
+ "chi = np.abs(Phi - Psi)\n",
+ "\n",
+ "# Map to box\n",
+ "chi_box = np.zeros((grid_size,grid_size,grid_size))\n",
+ "chi_box[(x_box,y_box,z_box)] = chi\n",
+ "\n",
+ "print(Phi.max(),Phi.min())\n",
+ "print(Psi.max(),Psi.min())\n",
+ "print(chi.max(),chi.min(),np.mean((chi)) )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Select a slice of the data\n",
+ "s0_slice = s0_box[:,:,50]\n",
+ "s0_slice+= 1E-10 # Add a small number to visualise in log scale.\n",
+ "\n",
+ "chi_slice= chi_box[:,:,50]\n",
+ "\n",
+ "#Figs\n",
+ "fig= plt.figure(figsize=(10, 6)) \n",
+ "gs = gridspec.GridSpec(1,2, wspace=0.5) \n",
+ "\n",
+ "ax0 = plt.subplot(gs[0,0])\n",
+ "ax1 = plt.subplot(gs[0,1])\n",
+ "\n",
+ "fig_s0 = ax0.imshow(s0_slice , origin='lower', cmap='inferno', aspect='equal',\n",
+ " norm=LogNorm(vmin=2E-1,vmax=1*s0_slice.max()),interpolation='None')\n",
+ "\n",
+ "\n",
+ "fig_chi= ax1.imshow(chi_slice , origin='lower', cmap='jet', aspect='equal',\n",
+ " norm=LogNorm(vmin=1E-7,vmax=1*chi_slice.max()),interpolation='bicubic')\n",
+ "\n",
+ "#Colorbars\n",
+ "cx1 = fig.add_axes([0.45,0.25,0.02,.5])\n",
+ "cx2 = fig.add_axes([0.92,0.25,0.02,.5])\n",
+ "cb_s0= plt.colorbar(fig_s0, cax = cx1, orientation='vertical')\n",
+ "cb_Bv = plt.colorbar(fig_chi, cax = cx2, orientation='vertical')\n",
+ "cb_s0.ax.tick_params(labelsize=14, color='k')\n",
+ "cb_Bv.ax.tick_params(labelsize=14, color='k')\n",
+ "\n",
+ "plt.show()\n",
+ "#fig.savefig(\"2D_map_s0_z-1_dpi_200.pdf\",bbox_inches='tight',dpi=200)\n",
+ "plt.close()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Measure spectrum\n",
+ "mesh_chi = ArrayMesh( chi_box, Nmesh=grid_size, compensated=False, BoxSize=box_size )\n",
+ "r_chi = FFTPower( mesh_chi, mode='1d', kmin= k_fun )\n",
+ "k_bins= np.array( r_chi.power['k'] )\n",
+ "Pk_chi= np.array( r_chi.power['power'].real )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plots\n",
+ "font = {'size':14}\n",
+ "plt.rc('font', **font)\n",
+ "plt.rc('font', family='serif')\n",
+ "plt.rc('text', usetex=True)\n",
+ "fig = plt.figure(figsize=(4, 5)) \n",
+ "gs = gridspec.GridSpec(2,1, height_ratios=[4,1]) \n",
+ "\n",
+ "ax0 = plt.subplot(gs[0])\n",
+ "fig_Pk = ax0.scatter( k_bins, Pk_chi, edgecolors='b', s=20,marker='.',facecolors='b',zorder=-2 )\n",
+ "plt.tick_params(axis='y', direction='in', which='both', labelleft='on', labelright='off', right='on')\n",
+ "plt.tick_params(axis='x', direction='in', which='both', labeltop='off', labelbottom='on', top='on')\n",
+ "ax0.set_ylabel(r'$P(k)$')\n",
+ "ax0.set_xlabel(r'$k$ [$h$ Mpc$^{-1}$]')\n",
+ "ax0.set_xscale('log')\n",
+ "ax0.set_yscale('log')\n",
+ "ax0.set_xlim( [ k_fun, 1.5*k_Nyq ] )\n",
+ "plt.fill_betweenx(np.array([1E-15,1E1]),k_Nyq,20 , alpha=0.2,color='red',interpolate=True)\n",
+ "ax0.set_ylim([1E-15, 1E-4])\n",
+ "\n",
+ "plt.show()\n",
+ "#fig.savefig(\"Pk_gramses.pdf\", bbox_inches='tight')\n",
+ "plt.close()"
+ ]
+ }
+ ],
+ "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.8.5"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/cmc_gauge_matterpower_z_0.dat b/cmc_gauge_matterpower_z_0.dat
new file mode 100644
index 0000000000000000000000000000000000000000..b36349d8cd38d4c28b05a28f5789853d58022ec5
--- /dev/null
+++ b/cmc_gauge_matterpower_z_0.dat
@@ -0,0 +1,695 @@
+# k/h P
+ 0.100000E-03 0.495744E+06
+ 0.102020E-03 0.467647E+06
+ 0.104081E-03 0.441188E+06
+ 0.106184E-03 0.416272E+06
+ 0.108329E-03 0.392807E+06
+ 0.110517E-03 0.370708E+06
+ 0.112750E-03 0.349895E+06
+ 0.115027E-03 0.330290E+06
+ 0.117351E-03 0.311824E+06
+ 0.119722E-03 0.294430E+06
+ 0.122140E-03 0.278043E+06
+ 0.124608E-03 0.262604E+06
+ 0.127125E-03 0.248059E+06
+ 0.129693E-03 0.234353E+06
+ 0.132313E-03 0.221438E+06
+ 0.134986E-03 0.209266E+06
+ 0.137713E-03 0.197795E+06
+ 0.140495E-03 0.186983E+06
+ 0.143333E-03 0.176791E+06
+ 0.146228E-03 0.167183E+06
+ 0.149182E-03 0.158125E+06
+ 0.152196E-03 0.149586E+06
+ 0.155271E-03 0.141535E+06
+ 0.158407E-03 0.133944E+06
+ 0.161607E-03 0.126787E+06
+ 0.164872E-03 0.120039E+06
+ 0.168203E-03 0.113675E+06
+ 0.171601E-03 0.107674E+06
+ 0.175067E-03 0.102015E+06
+ 0.178604E-03 0.966773E+05
+ 0.182212E-03 0.916425E+05
+ 0.185893E-03 0.868932E+05
+ 0.189648E-03 0.824126E+05
+ 0.193479E-03 0.781851E+05
+ 0.197388E-03 0.741961E+05
+ 0.201375E-03 0.704318E+05
+ 0.205443E-03 0.668790E+05
+ 0.209594E-03 0.635258E+05
+ 0.213828E-03 0.603607E+05
+ 0.218147E-03 0.573729E+05
+ 0.222554E-03 0.545523E+05
+ 0.227050E-03 0.518893E+05
+ 0.231637E-03 0.493749E+05
+ 0.236316E-03 0.470005E+05
+ 0.241090E-03 0.447582E+05
+ 0.245960E-03 0.426402E+05
+ 0.250929E-03 0.406397E+05
+ 0.255998E-03 0.387497E+05
+ 0.261170E-03 0.369640E+05
+ 0.266446E-03 0.352766E+05
+ 0.271828E-03 0.336820E+05
+ 0.277319E-03 0.321748E+05
+ 0.282922E-03 0.307500E+05
+ 0.288637E-03 0.294031E+05
+ 0.294468E-03 0.281296E+05
+ 0.300417E-03 0.269253E+05
+ 0.306485E-03 0.257865E+05
+ 0.312677E-03 0.247094E+05
+ 0.318993E-03 0.236906E+05
+ 0.325437E-03 0.227268E+05
+ 0.332012E-03 0.218150E+05
+ 0.338719E-03 0.209523E+05
+ 0.345561E-03 0.201360E+05
+ 0.352542E-03 0.193636E+05
+ 0.359664E-03 0.186326E+05
+ 0.366930E-03 0.179408E+05
+ 0.374342E-03 0.172861E+05
+ 0.381904E-03 0.166664E+05
+ 0.389619E-03 0.160800E+05
+ 0.397490E-03 0.155249E+05
+ 0.405520E-03 0.149996E+05
+ 0.413712E-03 0.145022E+05
+ 0.422070E-03 0.140314E+05
+ 0.430596E-03 0.135857E+05
+ 0.439295E-03 0.131637E+05
+ 0.448169E-03 0.127642E+05
+ 0.457222E-03 0.123861E+05
+ 0.466459E-03 0.120281E+05
+ 0.475882E-03 0.116892E+05
+ 0.485496E-03 0.113686E+05
+ 0.495303E-03 0.110652E+05
+ 0.505309E-03 0.107783E+05
+ 0.515517E-03 0.105070E+05
+ 0.525931E-03 0.102506E+05
+ 0.536556E-03 0.100084E+05
+ 0.547395E-03 0.977977E+04
+ 0.558453E-03 0.956402E+04
+ 0.569734E-03 0.936050E+04
+ 0.581244E-03 0.916860E+04
+ 0.592986E-03 0.898774E+04
+ 0.604965E-03 0.881738E+04
+ 0.617186E-03 0.865702E+04
+ 0.629654E-03 0.850620E+04
+ 0.642374E-03 0.836449E+04
+ 0.655350E-03 0.823147E+04
+ 0.668589E-03 0.810679E+04
+ 0.682096E-03 0.799008E+04
+ 0.695875E-03 0.788102E+04
+ 0.709933E-03 0.777932E+04
+ 0.724274E-03 0.768470E+04
+ 0.738906E-03 0.759689E+04
+ 0.753832E-03 0.751566E+04
+ 0.769061E-03 0.744079E+04
+ 0.784597E-03 0.737204E+04
+ 0.800447E-03 0.730915E+04
+ 0.816617E-03 0.725191E+04
+ 0.833114E-03 0.720011E+04
+ 0.849944E-03 0.715355E+04
+ 0.867114E-03 0.711205E+04
+ 0.884631E-03 0.707544E+04
+ 0.902501E-03 0.704356E+04
+ 0.920733E-03 0.701628E+04
+ 0.939333E-03 0.699345E+04
+ 0.958309E-03 0.697494E+04
+ 0.977668E-03 0.696065E+04
+ 0.997418E-03 0.695046E+04
+ 0.101757E-02 0.694427E+04
+ 0.103812E-02 0.694200E+04
+ 0.105910E-02 0.694355E+04
+ 0.108049E-02 0.694884E+04
+ 0.110232E-02 0.695780E+04
+ 0.112459E-02 0.697035E+04
+ 0.114730E-02 0.698641E+04
+ 0.117048E-02 0.700593E+04
+ 0.119413E-02 0.702884E+04
+ 0.121825E-02 0.705508E+04
+ 0.124286E-02 0.708460E+04
+ 0.126797E-02 0.711735E+04
+ 0.129358E-02 0.715328E+04
+ 0.131971E-02 0.719235E+04
+ 0.134637E-02 0.723451E+04
+ 0.137357E-02 0.727972E+04
+ 0.140132E-02 0.732793E+04
+ 0.142963E-02 0.737912E+04
+ 0.145851E-02 0.743325E+04
+ 0.148797E-02 0.749027E+04
+ 0.151803E-02 0.755016E+04
+ 0.154870E-02 0.761291E+04
+ 0.157998E-02 0.767849E+04
+ 0.161190E-02 0.774691E+04
+ 0.164446E-02 0.781813E+04
+ 0.167768E-02 0.789216E+04
+ 0.171158E-02 0.796898E+04
+ 0.174615E-02 0.804858E+04
+ 0.178143E-02 0.813094E+04
+ 0.181741E-02 0.821604E+04
+ 0.185413E-02 0.830389E+04
+ 0.189158E-02 0.839445E+04
+ 0.192980E-02 0.848771E+04
+ 0.196878E-02 0.858366E+04
+ 0.200855E-02 0.868228E+04
+ 0.204913E-02 0.878354E+04
+ 0.209052E-02 0.888742E+04
+ 0.213276E-02 0.899391E+04
+ 0.217584E-02 0.910303E+04
+ 0.221979E-02 0.921476E+04
+ 0.226464E-02 0.932911E+04
+ 0.231039E-02 0.944608E+04
+ 0.235706E-02 0.956566E+04
+ 0.240467E-02 0.968785E+04
+ 0.245325E-02 0.981263E+04
+ 0.250281E-02 0.994001E+04
+ 0.255337E-02 0.100700E+05
+ 0.260495E-02 0.102025E+05
+ 0.265758E-02 0.103376E+05
+ 0.271126E-02 0.104752E+05
+ 0.276603E-02 0.106154E+05
+ 0.282191E-02 0.107580E+05
+ 0.287892E-02 0.109031E+05
+ 0.293708E-02 0.110507E+05
+ 0.299641E-02 0.112007E+05
+ 0.305694E-02 0.113531E+05
+ 0.311870E-02 0.115079E+05
+ 0.318170E-02 0.116651E+05
+ 0.324597E-02 0.118247E+05
+ 0.331155E-02 0.119867E+05
+ 0.337844E-02 0.121510E+05
+ 0.344669E-02 0.123176E+05
+ 0.351632E-02 0.124866E+05
+ 0.358735E-02 0.126577E+05
+ 0.365982E-02 0.128311E+05
+ 0.373376E-02 0.130067E+05
+ 0.380918E-02 0.131845E+05
+ 0.388613E-02 0.133643E+05
+ 0.396464E-02 0.135462E+05
+ 0.404473E-02 0.137302E+05
+ 0.412644E-02 0.139161E+05
+ 0.420980E-02 0.141039E+05
+ 0.429484E-02 0.142936E+05
+ 0.438160E-02 0.144852E+05
+ 0.447012E-02 0.146785E+05
+ 0.456042E-02 0.148736E+05
+ 0.465255E-02 0.150704E+05
+ 0.474653E-02 0.152687E+05
+ 0.484242E-02 0.154686E+05
+ 0.494024E-02 0.156699E+05
+ 0.504004E-02 0.158727E+05
+ 0.514186E-02 0.160767E+05
+ 0.524573E-02 0.162820E+05
+ 0.535170E-02 0.164885E+05
+ 0.545981E-02 0.166960E+05
+ 0.557011E-02 0.169045E+05
+ 0.568263E-02 0.171139E+05
+ 0.579743E-02 0.173241E+05
+ 0.591455E-02 0.175350E+05
+ 0.603403E-02 0.177465E+05
+ 0.615592E-02 0.179585E+05
+ 0.628028E-02 0.181708E+05
+ 0.640715E-02 0.183833E+05
+ 0.653658E-02 0.185960E+05
+ 0.666863E-02 0.188087E+05
+ 0.680335E-02 0.190213E+05
+ 0.694078E-02 0.192336E+05
+ 0.708100E-02 0.194455E+05
+ 0.722404E-02 0.196568E+05
+ 0.736998E-02 0.198674E+05
+ 0.751886E-02 0.200772E+05
+ 0.767075E-02 0.202860E+05
+ 0.782571E-02 0.204936E+05
+ 0.798380E-02 0.206998E+05
+ 0.814509E-02 0.209045E+05
+ 0.830963E-02 0.211075E+05
+ 0.847750E-02 0.213086E+05
+ 0.864875E-02 0.215077E+05
+ 0.882347E-02 0.217044E+05
+ 0.900171E-02 0.218987E+05
+ 0.918356E-02 0.220904E+05
+ 0.936908E-02 0.222791E+05
+ 0.955835E-02 0.224647E+05
+ 0.975144E-02 0.226470E+05
+ 0.994843E-02 0.228257E+05
+ 0.101494E-01 0.230007E+05
+ 0.103544E-01 0.231716E+05
+ 0.105636E-01 0.233383E+05
+ 0.107770E-01 0.235005E+05
+ 0.109947E-01 0.236579E+05
+ 0.112168E-01 0.238103E+05
+ 0.114434E-01 0.239574E+05
+ 0.116746E-01 0.240990E+05
+ 0.119104E-01 0.242348E+05
+ 0.121510E-01 0.243645E+05
+ 0.123965E-01 0.244879E+05
+ 0.126469E-01 0.246046E+05
+ 0.129024E-01 0.247145E+05
+ 0.131631E-01 0.248171E+05
+ 0.134290E-01 0.249124E+05
+ 0.137003E-01 0.249999E+05
+ 0.139770E-01 0.250793E+05
+ 0.142594E-01 0.251505E+05
+ 0.145474E-01 0.252131E+05
+ 0.148413E-01 0.252668E+05
+ 0.151411E-01 0.253115E+05
+ 0.154470E-01 0.253467E+05
+ 0.157590E-01 0.253723E+05
+ 0.160774E-01 0.253880E+05
+ 0.164022E-01 0.253935E+05
+ 0.167335E-01 0.253886E+05
+ 0.170716E-01 0.253731E+05
+ 0.174164E-01 0.253467E+05
+ 0.177683E-01 0.253093E+05
+ 0.181272E-01 0.252606E+05
+ 0.184934E-01 0.252005E+05
+ 0.188670E-01 0.251289E+05
+ 0.192481E-01 0.250455E+05
+ 0.196370E-01 0.249504E+05
+ 0.200337E-01 0.248433E+05
+ 0.204384E-01 0.247242E+05
+ 0.208513E-01 0.245931E+05
+ 0.212725E-01 0.244499E+05
+ 0.217022E-01 0.242948E+05
+ 0.221406E-01 0.241276E+05
+ 0.225879E-01 0.239486E+05
+ 0.230442E-01 0.237579E+05
+ 0.235097E-01 0.235555E+05
+ 0.239847E-01 0.233418E+05
+ 0.244692E-01 0.231169E+05
+ 0.249635E-01 0.228812E+05
+ 0.254678E-01 0.226349E+05
+ 0.259823E-01 0.223784E+05
+ 0.265072E-01 0.221122E+05
+ 0.270426E-01 0.218367E+05
+ 0.275889E-01 0.215525E+05
+ 0.281463E-01 0.212601E+05
+ 0.287149E-01 0.209601E+05
+ 0.292949E-01 0.206531E+05
+ 0.298867E-01 0.203399E+05
+ 0.304905E-01 0.200212E+05
+ 0.311064E-01 0.196978E+05
+ 0.317348E-01 0.193705E+05
+ 0.323759E-01 0.190401E+05
+ 0.330299E-01 0.187075E+05
+ 0.336972E-01 0.183736E+05
+ 0.343779E-01 0.180393E+05
+ 0.350724E-01 0.177057E+05
+ 0.357809E-01 0.173737E+05
+ 0.365038E-01 0.170431E+05
+ 0.372412E-01 0.167137E+05
+ 0.379935E-01 0.163869E+05
+ 0.387610E-01 0.160645E+05
+ 0.395440E-01 0.157482E+05
+ 0.403429E-01 0.154389E+05
+ 0.411579E-01 0.151375E+05
+ 0.419893E-01 0.148447E+05
+ 0.428375E-01 0.145613E+05
+ 0.437029E-01 0.142878E+05
+ 0.445858E-01 0.140245E+05
+ 0.454865E-01 0.137718E+05
+ 0.464053E-01 0.135295E+05
+ 0.473428E-01 0.132977E+05
+ 0.482992E-01 0.130761E+05
+ 0.492749E-01 0.128641E+05
+ 0.502703E-01 0.126615E+05
+ 0.512858E-01 0.124674E+05
+ 0.523219E-01 0.122812E+05
+ 0.533789E-01 0.121023E+05
+ 0.544572E-01 0.119299E+05
+ 0.555573E-01 0.117633E+05
+ 0.566796E-01 0.116017E+05
+ 0.578246E-01 0.114479E+05
+ 0.589927E-01 0.113040E+05
+ 0.601845E-01 0.111630E+05
+ 0.614003E-01 0.110170E+05
+ 0.626407E-01 0.108651E+05
+ 0.639061E-01 0.107097E+05
+ 0.651971E-01 0.105506E+05
+ 0.665142E-01 0.103850E+05
+ 0.678578E-01 0.102112E+05
+ 0.692287E-01 0.100281E+05
+ 0.706272E-01 0.983478E+04
+ 0.720539E-01 0.963055E+04
+ 0.735095E-01 0.941505E+04
+ 0.749945E-01 0.918842E+04
+ 0.765095E-01 0.895125E+04
+ 0.780551E-01 0.870447E+04
+ 0.796319E-01 0.844946E+04
+ 0.812406E-01 0.818789E+04
+ 0.828817E-01 0.792172E+04
+ 0.845561E-01 0.765309E+04
+ 0.862642E-01 0.738430E+04
+ 0.880068E-01 0.711775E+04
+ 0.897847E-01 0.685590E+04
+ 0.915985E-01 0.660126E+04
+ 0.934489E-01 0.635629E+04
+ 0.953367E-01 0.612332E+04
+ 0.972626E-01 0.590445E+04
+ 0.992274E-01 0.570141E+04
+ 0.101232E+00 0.551539E+04
+ 0.103277E+00 0.534698E+04
+ 0.105363E+00 0.519605E+04
+ 0.107492E+00 0.506182E+04
+ 0.109663E+00 0.494290E+04
+ 0.111879E+00 0.483737E+04
+ 0.114139E+00 0.474291E+04
+ 0.116445E+00 0.465683E+04
+ 0.118797E+00 0.457608E+04
+ 0.121197E+00 0.449733E+04
+ 0.123645E+00 0.441699E+04
+ 0.126143E+00 0.433146E+04
+ 0.128691E+00 0.423748E+04
+ 0.131291E+00 0.413263E+04
+ 0.133943E+00 0.401572E+04
+ 0.136649E+00 0.388710E+04
+ 0.139409E+00 0.374848E+04
+ 0.142226E+00 0.360265E+04
+ 0.145099E+00 0.345290E+04
+ 0.148030E+00 0.330267E+04
+ 0.151020E+00 0.315528E+04
+ 0.154071E+00 0.301400E+04
+ 0.157184E+00 0.288201E+04
+ 0.160359E+00 0.276214E+04
+ 0.163598E+00 0.265646E+04
+ 0.166903E+00 0.256570E+04
+ 0.170275E+00 0.248900E+04
+ 0.173715E+00 0.242401E+04
+ 0.177224E+00 0.236742E+04
+ 0.180804E+00 0.231551E+04
+ 0.184457E+00 0.226468E+04
+ 0.188183E+00 0.221159E+04
+ 0.191985E+00 0.215342E+04
+ 0.195863E+00 0.208825E+04
+ 0.199820E+00 0.201547E+04
+ 0.203856E+00 0.193610E+04
+ 0.207974E+00 0.185255E+04
+ 0.212176E+00 0.176810E+04
+ 0.216462E+00 0.168620E+04
+ 0.220835E+00 0.160999E+04
+ 0.225296E+00 0.154184E+04
+ 0.229847E+00 0.148293E+04
+ 0.234490E+00 0.143276E+04
+ 0.239227E+00 0.138982E+04
+ 0.244060E+00 0.135114E+04
+ 0.248990E+00 0.131382E+04
+ 0.254020E+00 0.127514E+04
+ 0.259152E+00 0.123331E+04
+ 0.264387E+00 0.118772E+04
+ 0.269728E+00 0.113899E+04
+ 0.275177E+00 0.108902E+04
+ 0.280736E+00 0.104004E+04
+ 0.286407E+00 0.994506E+03
+ 0.292193E+00 0.953689E+03
+ 0.298096E+00 0.917754E+03
+ 0.304118E+00 0.885937E+03
+ 0.310261E+00 0.856706E+03
+ 0.316529E+00 0.828334E+03
+ 0.322923E+00 0.799056E+03
+ 0.329447E+00 0.768318E+03
+ 0.336102E+00 0.736522E+03
+ 0.342892E+00 0.704753E+03
+ 0.349819E+00 0.674412E+03
+ 0.356885E+00 0.646350E+03
+ 0.364095E+00 0.620972E+03
+ 0.371450E+00 0.597761E+03
+ 0.378954E+00 0.575753E+03
+ 0.386609E+00 0.554105E+03
+ 0.394419E+00 0.532177E+03
+ 0.402387E+00 0.510140E+03
+ 0.410516E+00 0.488510E+03
+ 0.418809E+00 0.467881E+03
+ 0.427269E+00 0.448724E+03
+ 0.435901E+00 0.430857E+03
+ 0.444706E+00 0.413884E+03
+ 0.453690E+00 0.397390E+03
+ 0.462855E+00 0.381072E+03
+ 0.472206E+00 0.365086E+03
+ 0.481745E+00 0.349682E+03
+ 0.491477E+00 0.335042E+03
+ 0.501405E+00 0.321244E+03
+ 0.511534E+00 0.308076E+03
+ 0.521868E+00 0.295329E+03
+ 0.532411E+00 0.282952E+03
+ 0.543166E+00 0.270959E+03
+ 0.554139E+00 0.259463E+03
+ 0.565333E+00 0.248529E+03
+ 0.576753E+00 0.238073E+03
+ 0.588404E+00 0.228013E+03
+ 0.600291E+00 0.218314E+03
+ 0.612418E+00 0.208969E+03
+ 0.624789E+00 0.200013E+03
+ 0.637411E+00 0.191456E+03
+ 0.650288E+00 0.183257E+03
+ 0.663424E+00 0.175386E+03
+ 0.676826E+00 0.167829E+03
+ 0.690499E+00 0.160577E+03
+ 0.704448E+00 0.153618E+03
+ 0.718679E+00 0.146943E+03
+ 0.733197E+00 0.140541E+03
+ 0.748009E+00 0.134403E+03
+ 0.763119E+00 0.128520E+03
+ 0.778536E+00 0.122882E+03
+ 0.794263E+00 0.117480E+03
+ 0.810308E+00 0.112306E+03
+ 0.826677E+00 0.107350E+03
+ 0.843378E+00 0.102604E+03
+ 0.860415E+00 0.980600E+02
+ 0.877797E+00 0.937089E+02
+ 0.895529E+00 0.895431E+02
+ 0.913620E+00 0.855551E+02
+ 0.932076E+00 0.817376E+02
+ 0.950906E+00 0.780836E+02
+ 0.970115E+00 0.745865E+02
+ 0.989713E+00 0.712396E+02
+ 0.100971E+01 0.680370E+02
+ 0.103010E+01 0.649726E+02
+ 0.105091E+01 0.620408E+02
+ 0.107214E+01 0.592362E+02
+ 0.109380E+01 0.565535E+02
+ 0.111590E+01 0.539878E+02
+ 0.113844E+01 0.515341E+02
+ 0.116144E+01 0.491880E+02
+ 0.118490E+01 0.469448E+02
+ 0.120884E+01 0.448003E+02
+ 0.123326E+01 0.427504E+02
+ 0.125817E+01 0.407911E+02
+ 0.128359E+01 0.389186E+02
+ 0.130952E+01 0.371293E+02
+ 0.133597E+01 0.354195E+02
+ 0.136296E+01 0.337859E+02
+ 0.139049E+01 0.322252E+02
+ 0.141858E+01 0.307344E+02
+ 0.144724E+01 0.293103E+02
+ 0.147648E+01 0.279503E+02
+ 0.150630E+01 0.266514E+02
+ 0.153673E+01 0.254111E+02
+ 0.156778E+01 0.242268E+02
+ 0.159945E+01 0.230960E+02
+ 0.163176E+01 0.220165E+02
+ 0.166472E+01 0.209861E+02
+ 0.169835E+01 0.200025E+02
+ 0.173266E+01 0.190637E+02
+ 0.176766E+01 0.181677E+02
+ 0.180337E+01 0.173127E+02
+ 0.183980E+01 0.164968E+02
+ 0.187697E+01 0.157184E+02
+ 0.191489E+01 0.149757E+02
+ 0.195357E+01 0.142672E+02
+ 0.199304E+01 0.135914E+02
+ 0.203330E+01 0.129467E+02
+ 0.207437E+01 0.123319E+02
+ 0.211628E+01 0.117455E+02
+ 0.215903E+01 0.111863E+02
+ 0.220265E+01 0.106531E+02
+ 0.224714E+01 0.101447E+02
+ 0.229254E+01 0.965994E+01
+ 0.233885E+01 0.919783E+01
+ 0.238610E+01 0.875730E+01
+ 0.243430E+01 0.833739E+01
+ 0.248348E+01 0.793715E+01
+ 0.253365E+01 0.755570E+01
+ 0.258483E+01 0.719217E+01
+ 0.263705E+01 0.684574E+01
+ 0.269032E+01 0.651563E+01
+ 0.274466E+01 0.620110E+01
+ 0.280011E+01 0.590143E+01
+ 0.285668E+01 0.561593E+01
+ 0.291439E+01 0.534396E+01
+ 0.297326E+01 0.508488E+01
+ 0.303332E+01 0.483811E+01
+ 0.309460E+01 0.460307E+01
+ 0.315712E+01 0.437922E+01
+ 0.322090E+01 0.416604E+01
+ 0.328596E+01 0.396303E+01
+ 0.335234E+01 0.376973E+01
+ 0.342006E+01 0.358567E+01
+ 0.348916E+01 0.341043E+01
+ 0.355964E+01 0.324359E+01
+ 0.363155E+01 0.308476E+01
+ 0.370491E+01 0.293356E+01
+ 0.377976E+01 0.278964E+01
+ 0.385611E+01 0.265265E+01
+ 0.393401E+01 0.252226E+01
+ 0.401348E+01 0.239818E+01
+ 0.409456E+01 0.228008E+01
+ 0.417727E+01 0.216771E+01
+ 0.426166E+01 0.206077E+01
+ 0.434775E+01 0.195902E+01
+ 0.443558E+01 0.186221E+01
+ 0.452519E+01 0.177011E+01
+ 0.461660E+01 0.168248E+01
+ 0.470987E+01 0.159912E+01
+ 0.480501E+01 0.151983E+01
+ 0.490208E+01 0.144440E+01
+ 0.500111E+01 0.137266E+01
+ 0.510214E+01 0.130442E+01
+ 0.520521E+01 0.123952E+01
+ 0.531036E+01 0.117780E+01
+ 0.541763E+01 0.111911E+01
+ 0.552708E+01 0.106330E+01
+ 0.563873E+01 0.101023E+01
+ 0.575264E+01 0.959765E+00
+ 0.586885E+01 0.911787E+00
+ 0.598741E+01 0.866172E+00
+ 0.610836E+01 0.822806E+00
+ 0.623176E+01 0.781580E+00
+ 0.635765E+01 0.742390E+00
+ 0.648609E+01 0.705137E+00
+ 0.661711E+01 0.669728E+00
+ 0.675079E+01 0.636072E+00
+ 0.688716E+01 0.604083E+00
+ 0.702629E+01 0.573682E+00
+ 0.716823E+01 0.544790E+00
+ 0.731304E+01 0.517334E+00
+ 0.746077E+01 0.491243E+00
+ 0.761149E+01 0.466451E+00
+ 0.776525E+01 0.442893E+00
+ 0.792212E+01 0.420510E+00
+ 0.808217E+01 0.399244E+00
+ 0.824543E+01 0.379039E+00
+ 0.841200E+01 0.359844E+00
+ 0.858193E+01 0.341609E+00
+ 0.875530E+01 0.324287E+00
+ 0.893217E+01 0.307832E+00
+ 0.911261E+01 0.292202E+00
+ 0.929670E+01 0.277356E+00
+ 0.948451E+01 0.263255E+00
+ 0.967610E+01 0.249863E+00
+ 0.987158E+01 0.237144E+00
+ 0.100710E+02 0.225064E+00
+ 0.102744E+02 0.213593E+00
+ 0.104820E+02 0.202700E+00
+ 0.106938E+02 0.192356E+00
+ 0.109098E+02 0.182534E+00
+ 0.111302E+02 0.173208E+00
+ 0.113550E+02 0.164353E+00
+ 0.115844E+02 0.155946E+00
+ 0.118184E+02 0.147964E+00
+ 0.120572E+02 0.140386E+00
+ 0.123007E+02 0.133192E+00
+ 0.125492E+02 0.126363E+00
+ 0.128027E+02 0.119880E+00
+ 0.130614E+02 0.113727E+00
+ 0.133252E+02 0.107886E+00
+ 0.135944E+02 0.102341E+00
+ 0.138691E+02 0.970791E-01
+ 0.141492E+02 0.920846E-01
+ 0.144351E+02 0.873445E-01
+ 0.147267E+02 0.828460E-01
+ 0.150242E+02 0.785768E-01
+ 0.153277E+02 0.745254E-01
+ 0.156373E+02 0.706809E-01
+ 0.159532E+02 0.670327E-01
+ 0.162755E+02 0.635710E-01
+ 0.166043E+02 0.602863E-01
+ 0.169397E+02 0.571698E-01
+ 0.172819E+02 0.542128E-01
+ 0.176310E+02 0.514073E-01
+ 0.179872E+02 0.487456E-01
+ 0.183505E+02 0.462204E-01
+ 0.187212E+02 0.438249E-01
+ 0.190994E+02 0.415523E-01
+ 0.194853E+02 0.393966E-01
+ 0.198789E+02 0.373516E-01
+ 0.202805E+02 0.354118E-01
+ 0.206902E+02 0.335719E-01
+ 0.211081E+02 0.318267E-01
+ 0.215346E+02 0.301714E-01
+ 0.219696E+02 0.286015E-01
+ 0.224134E+02 0.271125E-01
+ 0.228662E+02 0.257004E-01
+ 0.233281E+02 0.243611E-01
+ 0.237994E+02 0.230911E-01
+ 0.242802E+02 0.218867E-01
+ 0.247707E+02 0.207445E-01
+ 0.252710E+02 0.196615E-01
+ 0.257816E+02 0.186345E-01
+ 0.263024E+02 0.176608E-01
+ 0.268337E+02 0.167374E-01
+ 0.273758E+02 0.158620E-01
+ 0.279288E+02 0.150319E-01
+ 0.284930E+02 0.142450E-01
+ 0.290686E+02 0.134989E-01
+ 0.296558E+02 0.127915E-01
+ 0.302549E+02 0.121209E-01
+ 0.308661E+02 0.114852E-01
+ 0.314897E+02 0.108825E-01
+ 0.321258E+02 0.103112E-01
+ 0.327748E+02 0.976971E-02
+ 0.334369E+02 0.925637E-02
+ 0.341123E+02 0.876980E-02
+ 0.348014E+02 0.830859E-02
+ 0.355045E+02 0.787145E-02
+ 0.362217E+02 0.745712E-02
+ 0.369535E+02 0.706442E-02
+ 0.377000E+02 0.669224E-02
+ 0.384616E+02 0.633951E-02
+ 0.392386E+02 0.600522E-02
+ 0.400312E+02 0.568842E-02
+ 0.408399E+02 0.538820E-02
+ 0.416649E+02 0.510370E-02
+ 0.425066E+02 0.483410E-02
+ 0.433653E+02 0.457863E-02
+ 0.442413E+02 0.433655E-02
+ 0.451350E+02 0.410717E-02
+ 0.460469E+02 0.388983E-02
+ 0.469770E+02 0.368390E-02
+ 0.479261E+02 0.348878E-02
+ 0.488942E+02 0.330390E-02
+ 0.498820E+02 0.312875E-02
+ 0.508896E+02 0.296280E-02
+ 0.519177E+02 0.280559E-02
+ 0.529665E+02 0.265664E-02
+ 0.540365E+02 0.251554E-02
+ 0.551281E+02 0.238187E-02
+ 0.562417E+02 0.225524E-02
+ 0.573779E+02 0.213528E-02
+ 0.585370E+02 0.202166E-02
+ 0.597195E+02 0.191403E-02
+ 0.609260E+02 0.181208E-02
+ 0.621568E+02 0.171552E-02
+ 0.634124E+02 0.162406E-02
+ 0.646934E+02 0.153744E-02
+ 0.660003E+02 0.145540E-02
+ 0.673336E+02 0.137770E-02
+ 0.686938E+02 0.130412E-02
+ 0.700816E+02 0.123443E-02
+ 0.714973E+02 0.116843E-02
+ 0.729416E+02 0.110593E-02
+ 0.744151E+02 0.104674E-02
+ 0.759184E+02 0.990688E-03
+ 0.774521E+02 0.937601E-03
+ 0.790167E+02 0.887326E-03
+ 0.806129E+02 0.839714E-03
+ 0.822415E+02 0.794626E-03
+ 0.839028E+02 0.751927E-03
+ 0.855978E+02 0.711493E-03
+ 0.873269E+02 0.673205E-03
+ 0.890911E+02 0.636955E-03
+ 0.908908E+02 0.602638E-03
+ 0.927270E+02 0.570153E-03
+ 0.946001E+02 0.539406E-03
+ 0.965112E+02 0.510306E-03
+ 0.984608E+02 0.482769E-03
+ 0.100450E+03 0.456712E-03
+ 0.102479E+03 0.432057E-03
+ 0.104549E+03 0.408731E-03