<jupyter_start><jupyter_text>Table of Contents1. [Density Grid from FFT_DATA](fft)2. [Density Grid from Particle Positions](pmesh)  This notebook shows how to calculate density fields from Indra data on the SciServer.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import indratools as indra
from pmesh.window import CIC
from pmesh.window import Affine
boxlen = 1000.
npart = 1024
runnum = 205 # pick an Indra volume number, between 128 (2_0_0) and 447 (6_7_7) (for now)
# For fun, what is the X,Y,Z of this run?
indra.get_xyz(runnum)<jupyter_output><empty_output><jupyter_text>1. Density Grid from FFT Data A coarse grid of the Fourier modes of the density field are saved as the simulation runs at 505 time steps. These can be inverse Fourier transformed to obtain a 128x128x128 grid of the density field.<jupyter_code>redshift = .5
fft_snaps,fft_redshifts,fft_times = indra.FFT_snapinfo()
tnum = (np.abs(fft_redshifts - redshift)).argmin() # snapshot closest to desired redshift
L = 128 # dimensions of FFT grid
L2 = L//2
ind_kxky = np.concatenate((np.arange(L2,L+1),np.arange(1,L2)))<jupyter_output><empty_output><jupyter_text>First we need to re-index the Indra FFT arrays to correspond to the *k* values expected by `np.fft`.<jupyter_code>kxi,kyi,kzi = indra.getkvals()
kmin = 2.*np.pi/boxlen
kx = kxi[np.ix_(ind_kxky,ind_kxky)]
ky = kyi[np.ix_(ind_kxky,ind_kxky)]
kz = kzi[np.ix_(ind_kxky,ind_kxky)]
k = kmin*np.sqrt(kx**2+ky**2+kz**2)
fft_re,fft_im,a = indra.getfft(runnum,tnum)
fft_re = fft_re[np.ix_(ind_kxky,ind_kxky)]
fft_im = fft_im[np.ix_(ind_kxky,ind_kxky)]<jupyter_output><empty_output><jupyter_text>Now we can inverse Fourier transform to get the density grid.<jupyter_code>dk = np.zeros(fft_re.shape,dtype=complex)
dk.real = fft_re
dk.imag = fft_im
dk *= L**3/npart**3 # normalization
dens = np.fft.ifftn(dk,s=(L,L,L))
delta = np.sqrt(dens.real*dens.real+dens.imag*dens.imag)<jupyter_output><empty_output><jupyter_text>Let's plot a slice!<jupyter_code>plt.figure(figsize=(8,8))
plt.imshow(np.log10(delta[30]+1),origin='lower',extent=(0,L,0,L))
plt.xlim(0,L)
plt.ylim(0,L);<jupyter_output><empty_output><jupyter_text>2. Density Grid from Particle Positions The density field can also be measured directly from the particle positions, at higher resolution than the FFT grid, using the public package [pmesh](https://github.com/rainwoodman/pmesh). In this example we will choose a Cloud-In-Cell interpolation.<jupyter_code>snapnum = 41 # z = 1
ngrid = 256
cell_len = np.float64(boxlen)/np.float64(ngrid)
pos = indra.getpos(runnum,snapnum)
mesh = np.zeros((ngrid,ngrid,ngrid))
posnorm = pos/cell_len # normalize positions
affine = Affine(ndim=3, period=ngrid) # takes care of periodic boundary conditions
CIC.paint(mesh,posnorm,transform=affine)
meanden = np.mean(mesh.flatten())
delta = mesh/meanden-1
plt.figure(figsize=(8,8))
plt.imshow(np.log10(delta[0]+1),origin='lower',extent=(0,ngrid,0,ngrid))
# Note that this puts z on x-axis and y on y-axis
plt.xlim(0,ngrid)
plt.ylim(0,ngrid);<jupyter_output><empty_output><jupyter_text>pmesh has a few other window functions, including Triangular-Shaped-Cloud (see the pmesh [documentation](https://github.com/rainwoodman/pmesh)).<jupyter_code>from pmesh.window import TSC
tmesh = np.zeros((ngrid,ngrid,ngrid))
TSC.paint(tmesh,posnorm,transform=affine)
delta_tsc = tmesh/np.mean(tmesh) - 1
plt.figure(figsize=(8,8))
plt.imshow(np.log10(delta_tsc[0]+1),origin='lower',extent=(0,ngrid,0,ngrid))
plt.xlim(0,ngrid)
plt.ylim(0,ngrid);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1323529411764706, 'fraction_numerical': 0.027997737556561087, 'mean_word_length': 6.80794701986755, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '630081', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1467, 'n_tokens_neox': 1391, 'n_words': 305, 'path': 'density_field_examples.ipynb', 'repo_name': 'bfalck/indra-tools', 'size': '3536'}