In [1]:
import yt
import numpy as np
import matplotlib.pyplot as plt

ppc  = 3.0856776e18

%matplotlib inline
In [2]:
def fft_comp(pf, irho, iu, nindex_rho, Nx, pos, rad, level):

    cube =  pf.covering_grid(level=level, left_edge=np.array([(pos[0]-rad)*ppc, (pos[1]-rad)*ppc, (pos[2]-rad)*ppc]), 
                             dims=np.array([Nx, Nx, Nx]))

    rho = cube[irho].d
    u   = cube[iu].d

    nx, ny, nz = rho.shape
    
    xcenter = (np.min(cube["x"]) + (np.max(cube["x"]) - np.min(cube["x"]))/2.0).in_units("pc").d
    ycenter = (np.min(cube["y"]) + (np.max(cube["y"]) - np.min(cube["y"]))/2.0).in_units("pc").d
    zcenter = (np.min(cube["z"]) + (np.max(cube["z"]) - np.min(cube["z"]))/2.0).in_units("pc").d
    
    print("I'm within the fft calculations. Min dx %2f"%(np.min(cube["dx"].in_units("pc")).d))
    print("The center of the box is at x = %.2f, y=%.2f, z=%.2f"%(xcenter, ycenter, zcenter))

    del cube

    # do the FFTs -- note that since our data is real, there will be
    # too much information here.  fftn puts the positive freq terms in
    # the first half of the axes -- that's what we keep.  Our
    # normalization has an '8' to account for this clipping to one
    # octant.
    ru = np.fft.fftn(rho**nindex_rho * u)[0:nx//2+1,0:ny//2+1,0:nz//2+1]
    ru = 8.0*ru/(nx*ny*nz)
    
    return np.abs(ru)**2
In [3]:
def get_spectrum(jeans_ref=32, Lmax=50, resolution=0.06, snp=9.):
    
    pf   = yt.load("jeans%i/StBx_M3_006pc_jeans%i_hdf5_plt_cnt_00%.2i"%(jeans_ref, jeans_ref, snp))
    t0   = 239.2    
    time = abs(239.2894 - pf.current_time.in_units("Myr").d)
    
    print("Running the calculation of energy spectrum. Box size %i, jeans length resolved with %i zones,"%(Lmax, jeans_ref))
    print("mapping the box in a UG with res %.2f, at time = %.2f"%(resolution, time))
    
    # Normally fixed values but I can modify them if I need
    ppc  = 3.0856776e18
    tnow = snp / 10.
    
    px0, py0, pz0 = 460., -380., 10.
    
    rad = Lmax/2.0   # in pc

    # Mean bulk velocity of the cloud.
    bvx_mean = 0.
    bvy_mean = 1.0
    bvz_mean = -1.

    px = px0+bvx_mean*tnow
    py = py0+bvy_mean*tnow
    pz = pz0+bvz_mean*tnow

    pos = np.array([px, py, pz])
    nx  = int(Lmax / resolution)

    nindex_rho = 1./2.

    Kk = np.zeros( (nx//2+1, nx//2+1, nx//2+1))
    
    dd    = pf.all_data()
    maxdx = np.max(dd["dx"]).d/ppc
    
    level = int(np.log2(maxdx/resolution-2)) + 1

    for vel in [("gas", "velocity_x"), ("gas", "velocity_y"),("gas", "velocity_z")]:
        Kk += 0.5*fft_comp(pf, ("gas", "density"), vel, nindex_rho, nx, pos, rad, level)

    # wavenumbers
    kx = np.fft.rfftfreq(nx)*nx/Lmax
    ky = np.fft.rfftfreq(nx)*nx/Lmax
    kz = np.fft.rfftfreq(nx)*nx/Lmax

    # physical limits to the wavenumbers
    kmin = np.min(1.0/Lmax)
    kmax = np.min(0.5*nx/Lmax)

    kbins = np.arange(kmin, kmax, kmin)
    N = len(kbins)

    # bin the Fourier KE into radial kbins
    kx3d, ky3d, kz3d = np.meshgrid(kx, ky, kz, indexing="ij")
    k = np.sqrt(kx3d**2 + ky3d**2 + kz3d**2)

    whichbin = np.digitize(k.flat, kbins)
    ncount   = np.bincount(whichbin)

    E_spectrum = np.zeros(len(ncount)-1)

    for n in range(1,len(ncount)):
        E_spectrum[n-1] = np.sum(Kk.flat[whichbin==n])

    k = 0.5*(kbins[0:N-1] + kbins[1:N])
    E_spectrum = E_spectrum[1:N]

    index = np.argmax(E_spectrum)
    kmax  = k[index]
    Emax  = E_spectrum[index]

    return k, kmax, E_spectrum, Emax
In [61]:
k, kmax, E_spectrum_J8_0,  Emax = get_spectrum(jeans_ref=8, Lmax=Lmax, resolution=resol, snp=0.)
k, kmax, E_spectrum_J32_0,  Emax = get_spectrum(jeans_ref=32, Lmax=Lmax, resolution=resol, snp=0.)
yt : [INFO     ] 2017-09-12 16:48:51,764 Parameters: current_time              = 7.55060795711e+15
yt : [INFO     ] 2017-09-12 16:48:51,765 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 16:48:51,766 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 16:48:51,767 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 16:48:51,769 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 8 zones,
mapping the box in a UG with res 0.06, at time = 0.03
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
yt : [INFO     ] 2017-09-12 16:50:16,036 Parameters: current_time              = 7.55060795711e+15
yt : [INFO     ] 2017-09-12 16:50:16,112 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 16:50:16,160 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 16:50:16,163 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 16:50:16,167 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 32 zones,
mapping the box in a UG with res 0.06, at time = 0.03
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
In [45]:
Lmax = 25
resol = 0.06

k, kmax, E_spectrum_J4_0,  Emax = get_spectrum(jeans_ref=4, Lmax=Lmax, resolution=resol, snp=0.)
k, kmax, E_spectrum_J4_10, Emax = get_spectrum(jeans_ref=4, Lmax=Lmax, resolution=resol, snp=10.)
k, kmax, E_spectrum_J4_20, Emax = get_spectrum(jeans_ref=4, Lmax=Lmax, resolution=resol, snp=20.)
k, kmax, E_spectrum_J4_30, Emax = get_spectrum(jeans_ref=4, Lmax=Lmax, resolution=resol, snp=30.)

k, kmax, E_spectrum_J8_0,  Emax = get_spectrum(jeans_ref=8, Lmax=Lmax, resolution=resol, snp=0.)
k, kmax, E_spectrum_J8_10, Emax = get_spectrum(jeans_ref=8, Lmax=Lmax, resolution=resol, snp=10.)
k, kmax, E_spectrum_J8_20, Emax = get_spectrum(jeans_ref=8, Lmax=Lmax, resolution=resol, snp=20.)
k, kmax, E_spectrum_J8_30, Emax = get_spectrum(jeans_ref=8, Lmax=Lmax, resolution=resol, snp=30.)

k, kmax, E_spectrum_J32_0,  Emax = get_spectrum(jeans_ref=32, Lmax=Lmax, resolution=resol, snp=0.)
k, kmax, E_spectrum_J32_10, Emax = get_spectrum(jeans_ref=32, Lmax=Lmax, resolution=resol, snp=10.)
yt : [INFO     ] 2017-09-12 15:31:27,232 Parameters: current_time              = 7.5514001772e+15
yt : [INFO     ] 2017-09-12 15:31:27,233 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:31:27,234 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:31:27,235 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:31:27,236 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 4 zones,
mapping the box in a UG with res 0.06, at time = 0.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
yt : [INFO     ] 2017-09-12 15:32:52,232 Parameters: current_time              = 7.58294610627e+15
yt : [INFO     ] 2017-09-12 15:32:52,233 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:32:52,271 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:32:52,346 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:32:52,349 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 4 zones,
mapping the box in a UG with res 0.06, at time = 1.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
yt : [INFO     ] 2017-09-12 15:34:15,370 Parameters: current_time              = 7.61450049514e+15
yt : [INFO     ] 2017-09-12 15:34:15,371 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:34:15,372 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:34:15,373 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:34:15,374 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 4 zones,
mapping the box in a UG with res 0.06, at time = 2.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-378.16, z=7.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-378.16, z=7.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-378.16, z=7.84
yt : [INFO     ] 2017-09-12 15:35:33,351 Parameters: current_time              = 7.64606857208e+15
yt : [INFO     ] 2017-09-12 15:35:33,352 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:35:33,353 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:35:33,354 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:35:33,355 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 4 zones,
mapping the box in a UG with res 0.06, at time = 3.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-377.16, z=6.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-377.16, z=6.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-377.16, z=6.84
yt : [INFO     ] 2017-09-12 15:36:55,157 Parameters: current_time              = 7.55140017822e+15
yt : [INFO     ] 2017-09-12 15:36:55,158 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:36:55,158 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:36:55,159 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:36:55,160 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 8 zones,
mapping the box in a UG with res 0.06, at time = 0.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
yt : [INFO     ] 2017-09-12 15:38:11,864 Parameters: current_time              = 7.58295287081e+15
yt : [INFO     ] 2017-09-12 15:38:11,865 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:38:11,866 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:38:11,867 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:38:11,867 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 8 zones,
mapping the box in a UG with res 0.06, at time = 1.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
yt : [INFO     ] 2017-09-12 15:39:23,611 Parameters: current_time              = 7.6145054263e+15
yt : [INFO     ] 2017-09-12 15:39:23,612 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:39:23,613 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:39:23,614 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:39:23,615 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 8 zones,
mapping the box in a UG with res 0.06, at time = 2.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-378.16, z=7.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-378.16, z=7.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-378.16, z=7.84
yt : [INFO     ] 2017-09-12 15:40:20,450 Parameters: current_time              = 7.64606005105e+15
yt : [INFO     ] 2017-09-12 15:40:20,451 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:40:20,452 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:40:20,453 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:40:20,454 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 8 zones,
mapping the box in a UG with res 0.06, at time = 3.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-377.16, z=6.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-377.16, z=6.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-377.16, z=6.84
yt : [INFO     ] 2017-09-12 15:41:29,750 Parameters: current_time              = 7.55140394303e+15
yt : [INFO     ] 2017-09-12 15:41:29,751 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:41:29,752 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:41:29,753 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:41:29,754 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 32 zones,
mapping the box in a UG with res 0.06, at time = 0.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-380.16, z=9.84
yt : [INFO     ] 2017-09-12 15:42:41,675 Parameters: current_time              = 7.58294922254e+15
yt : [INFO     ] 2017-09-12 15:42:41,676 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 15:42:41,677 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:42:41,678 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 15:42:41,679 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 32 zones,
mapping the box in a UG with res 0.06, at time = 1.00
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.16, z=8.84
In [6]:
deltaE_J4_1  = E_spectrum_J4_0[0] - E_spectrum_J4_10[0]
deltaE_J8_1  = E_spectrum_J8_0[0] - E_spectrum_J8_10[0]
deltaE_J32_1 = E_spectrum_J32_0[0] - E_spectrum_J32_10[0]

deltaE_J4_2  = E_spectrum_J4_0[0] - E_spectrum_J4_20[0]
deltaE_J8_2  = E_spectrum_J8_0[0] - E_spectrum_J8_20[0]

deltaE_J4_3  = E_spectrum_J4_0[0] - E_spectrum_J4_30[0]
deltaE_J8_3  = E_spectrum_J8_0[0] - E_spectrum_J8_30[0]
In [68]:
fig = plt.figure(figsize=(10, 6))

ax = fig.add_axes([0.127, 0.13, 0.85, 0.85])

ax.loglog(k, E_spectrum_J4_0,  "k", linewidth=2,       label="$\lambda_{J} = 4 \Delta$x, t=0 Myr")
ax.loglog(k, E_spectrum_J8_0,  "b", linewidth=2,       label="$\lambda_{J} = 8 \Delta$x, t=0 Myr")
ax.loglog(k, E_spectrum_J32_0, "r", linewidth=2,       label="$\lambda_{J} = 32 \Delta$x, t=0 Myr")


ax.loglog(k, E_spectrum_J4_0[0]*(k/kmax)**(-5./3.), ls="-.", color="0.5")
ax.text(1.7e-1, 1.0e-13, "k$^{-5/3}$", fontsize=15)

ax.loglog([1.0/0.6, 1.0/0.6], [1.0e-16, 1.0e-11], "--k", linewidth=1)
ax.text(1.0/0.6, 1.e-13, "10$\\Delta$x", fontsize=15, )

ax.loglog([1.0/0.30, 1.0/0.30], [1.0e-16, 1.0e-11], ":k", linewidth=1)
ax.text(1.0/0.30, 1.e-13, "5$\\Delta$x", fontsize=15)

ax.set_ylim(1.0e-16, 4e-12)

ax.set_xlabel(r"$k$[pc$^{-1}$]",fontsize=25)
ax.set_ylabel(r"$E(k)dk$",fontsize=25)

ax.xaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.xaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.yaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.yaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.legend(loc=3, fontsize=20)

fig.show()

fig.savefig("/home/jcibanezm/codes/run/AccretionPaper/Figures/TurbulentEnergySpectrum_t0.pdf", format="pdf")
/home/jcibanezm/codes/libs/yt3.4/yt-conda/lib/python3.6/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
In [65]:
k, kmax, E_spectrum_J4_10, Emax = get_spectrum(jeans_ref=4, Lmax=Lmax, resolution=resol, snp=11.)
yt : [INFO     ] 2017-09-12 17:27:16,036 Parameters: current_time              = 7.58611662968e+15
yt : [INFO     ] 2017-09-12 17:27:16,037 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-09-12 17:27:16,038 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-09-12 17:27:16,039 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-09-12 17:27:16,040 Parameters: cosmological_simulation   = 0.0
Running the calculation of energy spectrum. Box size 25, jeans length resolved with 4 zones,
mapping the box in a UG with res 0.06, at time = 1.10
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.06, z=8.74
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.06, z=8.74
I'm within the fft calculations. Min dx 0.059340
The center of the box is at x = 459.84, y=-379.06, z=8.74
In [69]:
fig = plt.figure(figsize=(10, 6))

ax = fig.add_axes([0.127, 0.13, 0.85, 0.85])

ax.loglog(k, E_spectrum_J4_0,  ":k", linewidth=2)
ax.loglog(k, E_spectrum_J8_0,  ":b", linewidth=2)
ax.loglog(k, E_spectrum_J32_0,  ":r", linewidth=2)

ax.loglog(k, E_spectrum_J4_10+deltaE_J4_1*(k/kmax)**(-5./3.), "k", linewidth=2,label="$\lambda_{J} = 4 \Delta$x, t=1 Myr")
ax.loglog(k, E_spectrum_J8_10+deltaE_J8_1*(k/kmax)**(-5./3.), "b", linewidth=2,label="$\lambda_{J} = 8 \Delta$x, t=1 Myr")
ax.loglog(k, E_spectrum_J32_10+deltaE_J32_1*(k/kmax)**(-5./3.), "r", linewidth=2,label="$\lambda_{J} = 32 \Delta$x, t=1 Myr")


ax.loglog(k, E_spectrum_J4_0[0]*(k/kmax)**(-5./3.), ls="-.", color="0.5")
ax.text(1.7e-1, 1.0e-13, "k$^{-5/3}$", fontsize=15)

ax.loglog([1.0/0.6, 1.0/0.6], [1.0e-16, 1.0e-11], "--k", linewidth=1)
ax.text(1.0/0.6, 1.e-13, "10$\\Delta$x", fontsize=15, )

ax.loglog([1.0/0.30, 1.0/0.30], [1.0e-16, 1.0e-11], ":k", linewidth=1)
ax.text(1.0/0.30, 1.e-13, "5$\\Delta$x", fontsize=15)

ax.set_ylim(1.0e-16, 4e-12)

ax.set_xlabel(r"$k$[pc$^{-1}$]",fontsize=25)
ax.set_ylabel(r"$E(k)dk$",fontsize=25)

ax.xaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.xaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.yaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.yaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.legend(loc=3, fontsize=20)

fig.show()

fig.savefig("/home/jcibanezm/codes/run/AccretionPaper/Figures/TurbulentEnergySpectrum_t1.pdf", format="pdf")
/home/jcibanezm/codes/libs/yt3.4/yt-conda/lib/python3.6/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
In [70]:
fig = plt.figure(figsize=(10, 6))

ax = fig.add_axes([0.127, 0.13, 0.85, 0.85])

ax.loglog(k, E_spectrum_J4_0,  ":k", linewidth=2)
ax.loglog(k, E_spectrum_J8_0,  ":b", linewidth=2)

ax.loglog(k, E_spectrum_J4_20+deltaE_J4_2*(k/kmax)**(-5./3.), "k", linewidth=2,label="$\lambda_{J} = 4 \Delta$x, t=2 Myr")
ax.loglog(k, E_spectrum_J8_20+deltaE_J8_2*(k/kmax)**(-5./3.), "b", linewidth=2,label="$\lambda_{J} = 8 \Delta$x, t=2 Myr")

ax.loglog(k, E_spectrum_J4_0[0]*(k/kmax)**(-5./3.), ls="-.", color="0.5")
ax.text(1.7e-1, 1.0e-13, "k$^{-5/3}$", fontsize=15)

ax.loglog([1.0/0.6, 1.0/0.6], [1.0e-16, 1.0e-11], "--k", linewidth=1)
ax.text(1.0/0.6, 1.e-13, "10$\\Delta$x", fontsize=15, )

ax.loglog([1.0/0.30, 1.0/0.30], [1.0e-16, 1.0e-11], ":k", linewidth=1)
ax.text(1.0/0.30, 1.e-13, "5$\\Delta$x", fontsize=15)

ax.set_ylim(1.0e-16, 4e-12)

ax.set_xlabel(r"$k$[pc$^{-1}$]",fontsize=25)
ax.set_ylabel(r"$E(k)dk$",fontsize=25)

ax.xaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.xaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.yaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.yaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.legend(loc=3, fontsize=20)

fig.show()

fig.savefig("/home/jcibanezm/codes/run/AccretionPaper/Figures/TurbulentEnergySpectrum_t2.pdf", format="pdf")
/home/jcibanezm/codes/libs/yt3.4/yt-conda/lib/python3.6/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
In [71]:
fig = plt.figure(figsize=(10, 6))

ax = fig.add_axes([0.127, 0.13, 0.85, 0.85])

ax.loglog(k, E_spectrum_J4_0,  ":k", linewidth=2)
ax.loglog(k, E_spectrum_J8_0,  ":b", linewidth=2)

ax.loglog(k, E_spectrum_J4_30+deltaE_J4_3*(k/kmax)**(-5./3.), "k", linewidth=2,label="$\lambda_{J} = 4 \Delta$x, t=3 Myr")
ax.loglog(k, E_spectrum_J8_30+deltaE_J8_3*(k/kmax)**(-5./3.), "b", linewidth=2,label="$\lambda_{J} = 8 \Delta$x, t=3 Myr")

ax.loglog(k, E_spectrum_J4_0[0]*(k/kmax)**(-5./3.), ls="-.", color="0.5")
ax.text(1.7e-1, 1.0e-13, "k$^{-5/3}$", fontsize=15)

ax.loglog([1.0/0.6, 1.0/0.6], [1.0e-16, 1.0e-11], "--k", linewidth=1)
ax.text(1.0/0.6, 1.e-13, "10$\\Delta$x", fontsize=15, )

ax.loglog([1.0/0.30, 1.0/0.30], [1.0e-16, 1.0e-11], ":k", linewidth=1)
ax.text(1.0/0.30, 1.e-13, "5$\\Delta$x", fontsize=15)

ax.set_ylim(1.0e-16, 4e-12)

ax.set_xlabel(r"$k$[pc$^{-1}$]",fontsize=25)
ax.set_ylabel(r"$E(k)dk$",fontsize=25)

ax.xaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.xaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.yaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.yaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.legend(loc=3, fontsize=20)

fig.show()

fig.savefig("/home/jcibanezm/codes/run/AccretionPaper/Figures/TurbulentEnergySpectrum_t3.pdf", format="pdf")
/home/jcibanezm/codes/libs/yt3.4/yt-conda/lib/python3.6/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
In [21]:
fig = plt.figure(figsize=(12, 10))

ax = fig.add_axes([0.1, 0.1, 0.85, 0.85])

ax.loglog(k, E_spectrum_J4_0,  "k", linewidth=2,       label="$\lambda_{J} = 4 \Delta$x, t=0 Myr")
ax.loglog(k, E_spectrum_J4_10+deltaE_J4_1*(k/kmax)**(-5./3.), "orange", linewidth=2,  label="$\lambda_{J} = 4 \Delta$x, t=1 Myr")
ax.loglog(k, E_spectrum_J4_20+deltaE_J4_2*(k/kmax)**(-5./3.), "b", linewidth=2,       label="$\lambda_{J} = 4 \Delta$x, t=2 Myr")
ax.loglog(k, E_spectrum_J4_30+deltaE_J4_3*(k/kmax)**(-5./3.), "r", linewidth=2,       label="$\lambda_{J} = 4 \Delta$x, t=3 Myr")

ax.loglog(k,  E_spectrum_J4_0[0]*(k/kmax)**(-5./3.), ls=":", color="0.5")

ax.set_xlabel(r"$k$[pc$^{-1}$]",fontsize=25)
ax.set_ylabel(r"$E(k)dk$",fontsize=25)

ax.xaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.xaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.yaxis.set_tick_params(which="minor", length=4, width=2, direction="in")
ax.yaxis.set_tick_params(which="major", length=8, width=2, direction="in", labelsize=20)

ax.legend(loc=0, fontsize=20)

fig.show()

#fig.savefig("/home/jcibanezm/codes/run/AccretionPaper/Figures/TurbulentEnergySpectrum_Jeans4.pdf", format="pdf")
/home/jcibanezm/codes/libs/yt3.4/yt-conda/lib/python3.6/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "