[Examples] Turbulent cloud collapse

[tdavidcl]

No description provided.

[gemini-code-assist]

Summary of Changes

Hello @tdavidcl, I'm Gemini Code Assist¹! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed!

This pull request significantly enhances the simulation capabilities by introducing a comprehensive example for turbulent cloud collapse. It integrates a new, more advanced Equation of State (Machida 2006) and provides robust Perlin noise generation, which is leveraged to set up realistic turbulent initial velocity fields for the simulation. The changes span core physics definitions, solver configurations, and user-facing Python scripts for setup and analysis.

Highlights

• Turbulent Cloud Collapse Example: A new example script, run_collapse.py, has been added to simulate turbulent cloud collapse, demonstrating the setup, evolution, and analysis of such astrophysical phenomena within the framework.
• Machida 2006 Equation of State (EOS) Implementation: The Machida 2006 EOS has been integrated, allowing for more realistic modeling of gas behavior under varying density conditions, crucial for simulations like star formation.
• Perlin Noise Generation: New C++ and Python functionalities for generating Perlin noise have been introduced, providing a tool for creating natural-looking, pseudo-random initial conditions, such as turbulent velocity fields.

Using Gemini Code Assist

The full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips.

Invoking Gemini

You can request assistance from Gemini at any point by creating a comment using either /gemini <command> or @gemini-code-assist <command>. Below is a summary of the supported commands on the current page.

Feature	Command	Description
Code Review	/gemini review	Performs a code review for the current pull request in its current state.
Pull Request Summary	/gemini summary	Provides a summary of the current pull request in its current state.
Comment	@gemini-code-assist	Responds in comments when explicitly tagged, both in pull request comments and review comments.
Help	/gemini help	Displays a list of available commands.

Customization

To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a .gemini/ folder in the base of the repository. Detailed instructions can be found here.

Limitations & Feedback

Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here.

You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension.

Footnotes

1. Review the Privacy Notices, Generative AI Prohibited Use Policy, Terms of Service, and learn how to configure Gemini Code Assist in GitHub here. Gemini can make mistakes, so double check it and use code with caution. ↩

[github-actions]

Thanks @tdavidcl for opening this PR!

You can do multiple things directly here:
1 - Comment pre-commit.ci run to run pre-commit checks.
2 - Comment pre-commit.ci autofix to apply fixes.
3 - Add label autofix.ci to fix authorship & pre-commit for every commit made.
4 - Add label light-ci to only trigger a reduced & faster version of the CI (need the full one before merge).
5 - Add label trigger-ci to create an empty commit to trigger the CI.

Once the workflow completes a message will appear displaying informations related to the run.

Also the PR gets automatically reviewed by gemini, you can:
1 - Comment /gemini review to trigger a review
2 - Comment /gemini summary for a summary
3 - Tag it using @gemini-code-assist either in the PR or in review comments on files

[gemini-code-assist]

Code Review

This pull request introduces a new example for simulating turbulent cloud collapse, along with the necessary backend support for a new equation of state (Machida06) and Perlin noise generation. The changes are extensive, touching both Python example scripts and the core C++ source code.

My review focuses on improving the Python scripts' efficiency, maintainability, and adherence to best practices by addressing code duplication and inefficient patterns. In the C++ code, I've identified a potentially problematic file that seems to be obsolete and could cause issues, along with a minor cleanup in a lambda capture list. The core C++ changes for the new EOS appear to be well-implemented.

[gemini-code-assist]

This file appears to be either broken or obsolete and should likely be removed. It has several critical issues:

1. Memory Leak: A global int* perm is allocated with new int[512] on line 42 but is never deleted, causing a memory leak.
2. Uninitialized Data: The perm array is populated by an init() function, but it's not clear if or where this function is called. If not called, the noise functions will use uninitialized memory, leading to undefined behavior.
3. Compilation Errors: The functions dot() and fastfloor() are used but not defined within this file or its namespace. They seem to be part of the PerlinNoise class in Perlin.h, but this file uses a different noise namespace and cannot access them directly.
4. Inconsistency with Header: This file's structure (namespace noise with free functions) is inconsistent with Perlin.h (class shammath::PerlinNoise). The header seems to be a modern, self-contained, SYCL-compatible implementation, while this .cpp file is C-style and problematic.

Given that src/shammath/include/shammath/Perlin.h provides a complete inline implementation, this .cpp file seems unnecessary and harmful. Please consider removing it from the project.

[gemini-code-assist]

The use of global variables and the repeated creation of RegularGridInterpolator instances inside vel_field is inefficient and makes the code harder to maintain. Each call to vel_field will re-create three interpolator objects, which can be very slow if this function is called in a loop.

A better approach would be to encapsulate the velocity field data and the interpolators within a class. This avoids globals and ensures the expensive interpolator setup is only done once.

Here's a conceptual example of a refactored approach:

from scipy.interpolate import RegularGridInterpolator

class VelocityFieldInterpolator:
    def __init__(self, vx, vy, vz, domain_size=1.0):
        self.vx = vx
        self.vy = vy
        self.vz = vz
        self.domain_size = domain_size
        self.res = self.vx.shape[0]
        
        coords = numpy.linspace(0, self.domain_size, self.res)
        
        # Create interpolators once during initialization
        self.interp_vx = RegularGridInterpolator((coords, coords, coords), self.vx, 
                                                 bounds_error=False, fill_value=0.0)
        self.interp_vy = RegularGridInterpolator((coords, coords, coords), self.vy, 
                                                 bounds_error=False, fill_value=00)
        self.interp_vz = RegularGridInterpolator((coords, coords, coords), self.vz, 
                                                 bounds_error=False, fill_value=0.0)

    def __call__(self, pos):
        x, y, z = pos
        points = numpy.column_stack([numpy.atleast_1d(x), 
                                      numpy.atleast_1d(y), 
                                      numpy.atleast_1d(z)])
        
        vx_interp = self.interp_vx(points)[0]
        vy_interp = self.interp_vy(points)[0]
        vz_interp = self.interp_vz(points)[0]
        
        return vx_interp, vy_interp, vz_interp

# In the main execution block:
vx, vy, vz = TurbField(128, 2, 64, 0.7, 42)
vel_interpolator = VelocityFieldInterpolator(vx, vy, vz)
vel_x, vel_y, vel_z = vel_interpolator((0.5, 0.5, 0.5))

This change would significantly improve performance and code clarity.

[gemini-code-assist]

The function TurbField is duplicated from exemples/generate_vel_cube.py. This code duplication makes maintenance harder, as changes would need to be applied in two places. Instead of copying the function, you should import it from the generate_vel_cube module.

For example:

from generate_vel_cube import TurbField

# ... later in the code
vx,vy,vz = TurbField(128,2,64,0.7,seed)

This will require ensuring exemples is in the Python path, which is usually the case when running scripts from the root of the repository.

[gemini-code-assist]

The numpy module is imported twice. import numpy as np on line 5 is sufficient. This redundant import can be removed.

[gemini-code-assist]

The classes rho_xy_plot and v_slice_xy_plot share a lot of common functionality, such as initialization of paths, loading analysis data, getting dump IDs, and rendering all plots. This leads to significant code duplication.

Consider creating a base class to encapsulate the common logic, and have rho_xy_plot and v_slice_xy_plot inherit from it. This would improve code reuse and maintainability.

For example, a base class could contain methods like __init__, load_analysis, get_list_dumps_id, and render_all, while the subclasses would implement their specific compute and plot logic.

[gemini-code-assist]

The doxygen comment for the file appears to have a copy-paste error. The filename is Perlin.h, but the comment says @file paving_function.hpp. This should be corrected to reflect the actual filename.

Suggested change

	* @file paving_function.hpp
	* @file Perlin.h

[gemini-code-assist]

The lambda function captures mu, mh, and kb, but they are not used within the lambda for calculating pressure or sound speed with the Machida06 EOS. These unused captures can be removed to improve code clarity. While the compiler will likely optimize them away, it's good practice to only capture what's necessary.

                [rho_c1, rho_c2, rho_c3, cs0 = cs](

[github-actions]

Workflow report

workflow report corresponding to commit b3c4df8
Commiter email is timothee.davidcleris@proton.me

Light CI is enabled. This will only run the basic tests and not the full tests.
Merging a PR require the job "on PR / all" to pass which is disabled in this case.

Pre-commit check report

Some failures were detected in base source checks checks.
Check the On PR / Linting / Base source checks (pull_request) job in the tests for more detailled output

❌ ruff-format

2 files reformatted, 158 files left unchanged

❌ ruff-check

�[1m�[91mNPY002 �[0m�[1mReplace legacy `np.random.seed` call with `np.random.Generator`�[0m
  �[1m�[94m-->�[0m exemples/generate_vel_cube.py:18:9
   �[1m�[94m|�[0m
�[1m�[94m16 |�[0m     # apply ~k^-2 exp(-k^2/kmax^2) filter to white noise to get x, y, and z components of velocity field
�[1m�[94m17 |�[0m     for i in range(3):
�[1m�[94m18 |�[0m         numpy.random.seed(seed + i)
   �[1m�[94m|�[0m         �[1m�[91m^^^^^^^^^^^^^^^^^�[0m
�[1m�[94m19 |�[0m         rand_phase = fftpack.fftn(
�[1m�[94m20 |�[0m             numpy.random.normal(size=kSqr.shape)
   �[1m�[94m|�[0m

�[1m�[91mNPY002 �[0m�[1mReplace legacy `np.random.normal` call with `np.random.Generator`�[0m
  �[1m�[94m-->�[0m exemples/generate_vel_cube.py:20:13
   �[1m�[94m|�[0m
�[1m�[94m18 |�[0m         numpy.random.seed(seed + i)
�[1m�[94m19 |�[0m         rand_phase = fftpack.fftn(
�[1m�[94m20 |�[0m             numpy.random.normal(size=kSqr.shape)
   �[1m�[94m|�[0m             �[1m�[91m^^^^^^^^^^^^^^^^^^^�[0m
�[1m�[94m21 |�[0m         )  # fourier transform of white noise
�[1m�[94m22 |�[0m         vk = rand_phase * (float(minmode) / res) ** 2 / (kSqr + 1e-300)
   �[1m�[94m|�[0m

�[1m�[91mNPY002 �[0m�[1mReplace legacy `np.random.seed` call with `np.random.Generator`�[0m
   �[1m�[94m-->�[0m exemples/run_collapse.py:416:9
    �[1m�[94m|�[0m
�[1m�[94m414 |�[0m     # apply ~k^-2 exp(-k^2/kmax^2) filter to white noise to get x, y, and z components of velocity field
�[1m�[94m415 |�[0m     for i in range(3):
�[1m�[94m416 |�[0m         numpy.random.seed(seed + i)
    �[1m�[94m|�[0m         �[1m�[91m^^^^^^^^^^^^^^^^^�[0m
�[1m�[94m417 |�[0m         rand_phase = fftpack.fftn(
�[1m�[94m418 |�[0m             numpy.random.normal(size=kSqr.shape)
    �[1m�[94m|�[0m

�[1m�[91mNPY002 �[0m�[1mReplace legacy `np.random.normal` call with `np.random.Generator`�[0m
   �[1m�[94m-->�[0m exemples/run_collapse.py:418:13
    �[1m�[94m|�[0m
�[1m�[94m416 |�[0m         numpy.random.seed(seed + i)
�[1m�[94m417 |�[0m         rand_phase = fftpack.fftn(
�[1m�[94m418 |�[0m             numpy.random.normal(size=kSqr.shape)
    �[1m�[94m|�[0m             �[1m�[91m^^^^^^^^^^^^^^^^^^^�[0m
�[1m�[94m419 |�[0m         )  # fourier transform of white noise
�[1m�[94m420 |�[0m         vk = rand_phase * (float(minmode) / res) ** 2 / (kSqr + 1e-300)
    �[1m�[94m|�[0m

Found 9 errors (5 fixed, 4 remaining).

❌ trailing-whitespace

Fixing exemples/run_collapse.py
Fixing exemples/generate_vel_cube.py

❌ end-of-file-fixer

Fixing exemples/run_collapse.py

❌ Check license headers

The pre-commit checks have found some missing or ill formed license header.
All C++ files (headers or sources) should start with :

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2025 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//

Any line break before this header or change in its formatting will trigger the fail of this test
List of files with errors :

• /src/shammath/src/Perlin.cpp

❌ remove-tabs

Substituting tabs in: src/shammath/src/Perlin.cpp by 4 whitespaces

Tabs have been successfully removed. Now aborting the commit.
You can check the changes made. Then simply "git add --update ." and re-commit

❌ check-executables-have-shebangs

exemples/generate_vel_cube.py: marked executable but has no (or invalid) shebang!
  If it isn't supposed to be executable, try: `chmod -x exemples/generate_vel_cube.py`
  If on Windows, you may also need to: `git add --chmod=-x exemples/generate_vel_cube.py`
  If it is supposed to be executable, double-check its shebang.
src/shammath/src/Perlin.cpp: marked executable but has no (or invalid) shebang!
  If it isn't supposed to be executable, try: `chmod -x src/shammath/src/Perlin.cpp`
  If on Windows, you may also need to: `git add --chmod=-x src/shammath/src/Perlin.cpp`
  If it is supposed to be executable, double-check its shebang.

❌ forbid-tabs

Tabs detected in file: src/shammath/src/Perlin.cpp

❌ Check doxygen headers

The pre-commit checks have found some files without or with a wrong doxygen file header

A properly formed doxygen header should be like and placed just below #pragma once in headers or below the license in source files:

/**
 * @file <filename>
 * @author <Author name> (<email>)
 * @brief <Description of the file content>
 *
 */

This test is triggered if the doxygen header is missing or with the wrong file name

At some point we will refer to a guide in the doc about this

List of files with errors :

• /src/shammath/src/Perlin.cpp

Suggested changes

Detailed changes :

diff --git a/exemples/generate_vel_cube.py b/exemples/generate_vel_cube.py
index af0d7c71..2ab5b76e 100755
--- a/exemples/generate_vel_cube.py
+++ b/exemples/generate_vel_cube.py
@@ -1,11 +1,9 @@
-
-import numpy
 import matplotlib.pyplot as plt
-
+import numpy
 
 
 def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
-    from scipy import  fftpack
+    from scipy import fftpack
 
     freqs = fftpack.fftfreq(res)
     freq3d = numpy.array(numpy.meshgrid(freqs, freqs, freqs, indexing="ij"))
@@ -44,18 +42,14 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
         for j in range(3):
             if i == j:
                 vk_new[i] += sol_weight * VK[j]
-            vk_new[i] += (
-                (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j]
-            )
+            vk_new[i] += (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j]
     vk_new[:, kSqr == 0] = 0.0
     VK = vk_new
 
-    vel = numpy.array(
-        [fftpack.ifftn(vk).real for vk in VK]
-    )  # transform back to real space
+    vel = numpy.array([fftpack.ifftn(vk).real for vk in VK])  # transform back to real space
     vel -= numpy.average(vel, axis=(1, 2, 3))[:, numpy.newaxis, numpy.newaxis, numpy.newaxis]
     vel = vel / numpy.sqrt(numpy.sum(vel**2, axis=0).mean())  # normalize so that RMS is 1
-    return numpy.array(vel,dtype='f4')
+    return numpy.array(vel, dtype="f4")
 
 
 # Global variables for velocity field
@@ -64,66 +58,67 @@ vy_global = None
 vz_global = None
 domain_size_global = 1.0
 
+
 def vel_field(pos):
     """
     Interpolate velocity at position (x, y, z) using global velocity fields.
-    
+
     Parameters:
     -----------
     pos : tuple
         (x, y, z) position at which to interpolate the velocity
-    
+
     Returns:
     --------
     tuple
         (vx, vy, vz) velocity components at the given position(s)
     """
     from scipy.interpolate import RegularGridInterpolator
-    
+
     global vx_global, vy_global, vz_global, domain_size_global
-    
+
     x, y, z = pos
     res = vx_global.shape[0]
-    
+
     # Create coordinate arrays for the grid
     coords = numpy.linspace(0, domain_size_global, res)
-    
+
     # Create interpolators for each velocity component
-    interp_vx = RegularGridInterpolator((coords, coords, coords), vx_global, 
-                                         bounds_error=False, fill_value=0.0)
-    interp_vy = RegularGridInterpolator((coords, coords, coords), vy_global, 
-                                         bounds_error=False, fill_value=0.0)
-    interp_vz = RegularGridInterpolator((coords, coords, coords), vz_global, 
-                                         bounds_error=False, fill_value=0.0)
-    
-    points = numpy.column_stack([numpy.atleast_1d(x), 
-                                  numpy.atleast_1d(y), 
-                                  numpy.atleast_1d(z)])
-    
+    interp_vx = RegularGridInterpolator(
+        (coords, coords, coords), vx_global, bounds_error=False, fill_value=0.0
+    )
+    interp_vy = RegularGridInterpolator(
+        (coords, coords, coords), vy_global, bounds_error=False, fill_value=0.0
+    )
+    interp_vz = RegularGridInterpolator(
+        (coords, coords, coords), vz_global, bounds_error=False, fill_value=0.0
+    )
+
+    points = numpy.column_stack([numpy.atleast_1d(x), numpy.atleast_1d(y), numpy.atleast_1d(z)])
+
     vx_interp = interp_vx(points)[0]
     vy_interp = interp_vy(points)[0]
     vz_interp = interp_vz(points)[0]
-    
+
     return vx_interp, vy_interp, vz_interp
 
 
 if __name__ == "__main__":
-
     seed = 42
 
-    #avx,avy,avz = make_turb_field(res,power,seed)
-    vx,vy,vz = TurbField(128,2,64,0.7,seed)
+    # avx,avy,avz = make_turb_field(res,power,seed)
+    vx, vy, vz = TurbField(128, 2, 64, 0.7, seed)
 
     print(vx)
     print(vy)
     print(vz)
-    
+
     # Set global velocity fields
     vx_global = vx
     vy_global = vy
     vz_global = vz
     domain_size_global = 1.0
-    
+
     # Example: Interpolate velocity at a specific position
     test_pos = (0.5, 0.5, 0.5)
     vel_x, vel_y, vel_z = vel_field(test_pos)
diff --git a/exemples/run_collapse.py b/exemples/run_collapse.py
index cff417bf..2a8a1587 100644
--- a/exemples/run_collapse.py
+++ b/exemples/run_collapse.py
@@ -1,9 +1,10 @@
 import glob
 import json
+
 import matplotlib
 import matplotlib.pyplot as plt
-import numpy as np
 import numpy
+import numpy as np
 
 import shamrock
 
@@ -26,7 +27,7 @@ kg_m3_codeu = codeu.get("kg") * codeu.get("m", power=-3)
 print(f"m_s_codeu: {kg_m3_codeu}")
 
 cs = 190.0  # m/s
-rho_g = kg_m3_codeu * 1.e-18
+rho_g = kg_m3_codeu * 1.0e-18
 initial_u = 1.0
 Npart = int(0.1e6)
 ampl = 3
@@ -55,33 +56,25 @@ tsound = sim_radius / cs
 tff = shamrock.phys.free_fall_time(rho_g, codeu)
 
 P_init, _cs_init, T_init = shamrock.phys.eos.eos_Machida06(
-        cs=cs, rho=rho_g, rho_c1=rho_c1, rho_c2=rho_c2, rho_c3=rho_c3, mu=mu, mh=mh, kb=kb
-    )
+    cs=cs, rho=rho_g, rho_c1=rho_c1, rho_c2=rho_c2, rho_c3=rho_c3, mu=mu, mh=mh, kb=kb
+)
 
 print("---------------------------------------------------")
 print(f"Npart                         : {Npart}")
 print(f"sim_radius                    : {sim_radius} [pc]")
-print(f"total mass                    : {rho_g * (2*sim_radius)**3:.3e} [solar mass]")
+print(f"total mass                    : {rho_g * (2 * sim_radius) ** 3:.3e} [solar mass]")
 print(f"rho                           : {rho_g / kg_m3_codeu:.3e} [kg/m^3]")
-print(f"P                             : {P_init / codeu.get("Pa"):.3e} [Pa]")
+print(f"P                             : {P_init / codeu.get('Pa'):.3e} [Pa]")
 print(f"T                             : {T_init:.3e} [K]")
 print(f"cs                            : {cs / m_s_codeu:.3e} [m/s]")
 print(f"tsound (= R/cs)               : {tsound:9.1f} [years]")
 print(f"tff (= sqrt(3*pi/(32*G*rho))) : {tff:9.1f} [years]")
-print(f"tff/tsound                    : {tff/tsound:.4f} (<1 = collapse)")
+print(f"tff/tsound                    : {tff / tsound:.4f} (<1 = collapse)")
 print("---------------------------------------------------")
 
 
-
-
-
-
-
-
-
-
 class rho_xy_plot:
-    def __init__(self, model, center, ext_r, nx,ny, analysis_folder, analysis_prefix):
+    def __init__(self, model, center, ext_r, nx, ny, analysis_folder, analysis_prefix):
         self.model = model
         self.center = center
         self.ext_r = ext_r
@@ -98,7 +91,6 @@ class rho_xy_plot:
         self.img_glob = self.analysis_prefix + "*.png"
 
     def compute_rho_xy(self):
-        
         arr_rho_xy = model.render_cartesian_column_integ(
             "rho",
             "f64",
@@ -108,15 +100,18 @@ class rho_xy_plot:
             nx=self.nx,
             ny=self.ny,
         )
-        
+
         return arr_rho_xy
 
     def analysis_save(self, iplot):
         arr_rho_xy = self.compute_rho_xy()
         arr_rho_xy /= kg_m3_codeu
         if shamrock.sys.world_rank() == 0:
-            cx,cy,cz = self.center
-            metadata = {"extent": [cx-self.ext_r, cx+self.ext_r, cy-self.ext_r, cy+self.ext_r], "time": self.model.get_time()}
+            cx, cy, cz = self.center
+            metadata = {
+                "extent": [cx - self.ext_r, cx + self.ext_r, cy - self.ext_r, cy + self.ext_r],
+                "time": self.model.get_time(),
+            }
             np.save(self.npy_data_filename.format(iplot), arr_rho_xy)
 
             with open(self.json_data_filename.format(iplot), "w") as fp:
@@ -130,7 +125,6 @@ class rho_xy_plot:
     def plot_rho_xy(self, iplot):
         arr_rho_xy, metadata = self.load_analysis(iplot)
         if shamrock.sys.world_rank() == 0:
-
             # Reset the figure using the same memory as the last one
             plt.figure(num=1, clear=True, dpi=200)
 
@@ -150,7 +144,7 @@ class rho_xy_plot:
                 origin="lower",
                 extent=metadata["extent"],
                 norm="log",
-                vmin= min_val
+                vmin=min_val,
             )
 
             plt.xlabel("x [pc]")
@@ -164,7 +158,6 @@ class rho_xy_plot:
             plt.close()
 
     def get_list_dumps_id(self):
-
         list_files = glob.glob(self.json_glob)
         list_files.sort()
         list_json_files = []
@@ -179,9 +172,8 @@ class rho_xy_plot:
                 self.plot_rho_xy(iplot)
 
 
-
 class v_slice_xy_plot:
-    def __init__(self, model, center, ext_r, nx,ny, analysis_folder, analysis_prefix):
+    def __init__(self, model, center, ext_r, nx, ny, analysis_folder, analysis_prefix):
         self.model = model
         self.center = center
         self.ext_r = ext_r
@@ -198,7 +190,6 @@ class v_slice_xy_plot:
         self.img_glob = self.analysis_prefix + "*.png"
 
     def compute_v_slice(self):
-        
         arr_v_slice = model.render_cartesian_slice(
             "vxyz",
             "f64_3",
@@ -208,15 +199,18 @@ class v_slice_xy_plot:
             nx=self.nx,
             ny=self.ny,
         )
-        
+
         return arr_v_slice
 
     def analysis_save(self, iplot):
         arr_v_slice = self.compute_v_slice()
         arr_v_slice /= cs
         if shamrock.sys.world_rank() == 0:
-            cx,cy,cz = self.center
-            metadata = {"extent": [cx-self.ext_r, cx+self.ext_r, cy-self.ext_r, cy+self.ext_r], "time": self.model.get_time()}
+            cx, cy, cz = self.center
+            metadata = {
+                "extent": [cx - self.ext_r, cx + self.ext_r, cy - self.ext_r, cy + self.ext_r],
+                "time": self.model.get_time(),
+            }
             np.save(self.npy_data_filename.format(iplot), arr_v_slice)
 
             with open(self.json_data_filename.format(iplot), "w") as fp:
@@ -230,9 +224,10 @@ class v_slice_xy_plot:
     def plot_v_slice(self, iplot):
         arr_v_slice, metadata = self.load_analysis(iplot)
 
-        v_norm = np.sqrt(arr_v_slice[:, :, 0] ** 2 + arr_v_slice[:, :, 1] ** 2 + arr_v_slice[:, :, 2] ** 2)
+        v_norm = np.sqrt(
+            arr_v_slice[:, :, 0] ** 2 + arr_v_slice[:, :, 1] ** 2 + arr_v_slice[:, :, 2] ** 2
+        )
         if shamrock.sys.world_rank() == 0:
-
             # Reset the figure using the same memory as the last one
             plt.figure(num=1, clear=True, dpi=200)
 
@@ -242,12 +237,7 @@ class v_slice_xy_plot:
             my_cmap.set_bad(color="black")
 
             v_ext = 0.04
-            res = plt.imshow(
-                v_norm,
-                cmap=my_cmap,
-                origin="lower",
-                extent=metadata["extent"]
-            )
+            res = plt.imshow(v_norm, cmap=my_cmap, origin="lower", extent=metadata["extent"])
 
             plt.xlabel("x [pc]")
             plt.ylabel("y [pc]")
@@ -260,7 +250,6 @@ class v_slice_xy_plot:
             plt.close()
 
     def get_list_dumps_id(self):
-
         list_files = glob.glob(self.json_glob)
         list_files.sort()
         list_json_files = []
@@ -275,10 +264,9 @@ class v_slice_xy_plot:
                 self.plot_v_slice(iplot)
 
 
+dump_prefix = sim_directory + "/dump_"
 
 
-
-dump_prefix = sim_directory + "/dump_"
 def get_dump_name(idump):
     return dump_prefix + f"{idump:07}" + ".sham"
 
@@ -314,7 +302,6 @@ def purge_old_dumps():
             os.remove(f)
 
 
-
 bmin = (-sim_radius, -sim_radius, -sim_radius)
 bmax = (sim_radius, sim_radius, sim_radius)
 
@@ -353,7 +340,8 @@ model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4"
 
 def is_in_sphere(pt):
     x, y, z = pt
-    return (x**2 + y**2 + z**2) < sim_radius*sim_radius
+    return (x**2 + y**2 + z**2) < sim_radius * sim_radius
+
 
 def config_model():
     cfg = model.gen_default_config()
@@ -362,21 +350,22 @@ def config_model():
     cfg.set_artif_viscosity_VaryingCD10(
         alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
     )
-    #cfg.set_boundary_periodic()
+    # cfg.set_boundary_periodic()
     cfg.set_eos_machida06(rho_c1=rho_c1, rho_c2=rho_c2, rho_c3=rho_c3, cs=cs, mu=mu, mh=mh, kb=kb)
 
     cfg.set_self_gravity_fmm(order=5, opening_angle=0.5)
 
     cfg.set_units(codeu)
-    #cfg.print_status()
+    # cfg.print_status()
     model.set_solver_config(cfg)
     model.init_scheduler(scheduler_split_val, scheduler_merge_val)
 
-
-    bsim_min = (-sim_radius*2, -sim_radius*2, -sim_radius*2)
-    bsim_max = (sim_radius*2, sim_radius*2, sim_radius*2)
+    bsim_min = (-sim_radius * 2, -sim_radius * 2, -sim_radius * 2)
+    bsim_max = (sim_radius * 2, sim_radius * 2, sim_radius * 2)
     model.resize_simulation_box(bsim_min, bsim_max)
-    cfg.add_kill_sphere(center=(0, 0, 0), radius=sim_radius*2)  # kill particles outside the simulation box
+    cfg.add_kill_sphere(
+        center=(0, 0, 0), radius=sim_radius * 2
+    )  # kill particles outside the simulation box
 
     setup = model.get_setup()
     gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
@@ -395,13 +384,13 @@ def config_model():
 
     model.set_cfl_cour(0.1)
     model.set_cfl_force(0.1)
-    
+
     # if corrector is trigger CFL will become normal again after 1000 steps
-    model.set_cfl_mult_stiffness(5) 
+    model.set_cfl_mult_stiffness(5)
 
 
 def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
-    from scipy import  fftpack
+    from scipy import fftpack
 
     freqs = fftpack.fftfreq(res)
     freq3d = numpy.array(numpy.meshgrid(freqs, freqs, freqs, indexing="ij"))
@@ -440,26 +429,20 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
         for j in range(3):
             if i == j:
                 vk_new[i] += sol_weight * VK[j]
-            vk_new[i] += (
-                (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j]
-            )
+            vk_new[i] += (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j]
     vk_new[:, kSqr == 0] = 0.0
     VK = vk_new
 
-    vel = numpy.array(
-        [fftpack.ifftn(vk).real for vk in VK]
-    )  # transform back to real space
+    vel = numpy.array([fftpack.ifftn(vk).real for vk in VK])  # transform back to real space
     vel -= numpy.average(vel, axis=(1, 2, 3))[:, numpy.newaxis, numpy.newaxis, numpy.newaxis]
     vel = vel / numpy.sqrt(numpy.sum(vel**2, axis=0).mean())  # normalize so that RMS is 1
-    return numpy.array(vel,dtype='f4')
-
-
+    return numpy.array(vel, dtype="f4")
 
 
 seed = 42
 
-#avx,avy,avz = make_turb_field(res,power,seed)
-vx,vy,vz = TurbField(128,2,64,0.7,seed)
+# avx,avy,avz = make_turb_field(res,power,seed)
+vx, vy, vz = TurbField(128, 2, 64, 0.7, seed)
 
 print(f"vx shape: {vx.shape}, dtype: {vx.dtype}")
 print(f"vy shape: {vy.shape}, dtype: {vy.dtype}")
@@ -483,52 +466,54 @@ coords = numpy.linspace(0, domain_size_global, res)
 
 # Create interpolators for each velocity component
 from scipy.interpolate import RegularGridInterpolator
-interp_vx = RegularGridInterpolator((coords, coords, coords), vx_global, 
-                                        bounds_error=False, fill_value=0.0)
-interp_vy = RegularGridInterpolator((coords, coords, coords), vy_global, 
-                                        bounds_error=False, fill_value=0.0)
-interp_vz = RegularGridInterpolator((coords, coords, coords), vz_global, 
-                                        bounds_error=False, fill_value=0.0)
+
+interp_vx = RegularGridInterpolator(
+    (coords, coords, coords), vx_global, bounds_error=False, fill_value=0.0
+)
+interp_vy = RegularGridInterpolator(
+    (coords, coords, coords), vy_global, bounds_error=False, fill_value=0.0
+)
+interp_vz = RegularGridInterpolator(
+    (coords, coords, coords), vz_global, bounds_error=False, fill_value=0.0
+)
 
 
 def vel_field(pos):
     """
     Interpolate velocity at position (x, y, z) using global velocity fields.
-    
+
     Parameters:
     -----------
     pos : tuple
         (x, y, z) position in simulation coordinates
-    
+
     Returns:
     --------
     tuple
         (vx, vy, vz) velocity components at the given position(s)
     """
-    
-    
+
     global interp_vx, interp_vy, interp_vz, domain_size_global, sim_bmin, sim_bmax, ampl, cs
-    
+
     x, y, z = pos
-    
+
     # Transform from simulation coordinates to velocity field coordinates [0, 1]
     x_vf = (x - sim_bmin[0]) / (sim_bmax[0] - sim_bmin[0])
     y_vf = (y - sim_bmin[1]) / (sim_bmax[1] - sim_bmin[1])
     z_vf = (z - sim_bmin[2]) / (sim_bmax[2] - sim_bmin[2])
-    
-    
-    
-    points = numpy.column_stack([numpy.atleast_1d(x_vf), 
-                                  numpy.atleast_1d(y_vf), 
-                                  numpy.atleast_1d(z_vf)])
-    
+
+    points = numpy.column_stack(
+        [numpy.atleast_1d(x_vf), numpy.atleast_1d(y_vf), numpy.atleast_1d(z_vf)]
+    )
+
     vx_interp = interp_vx(points)[0]
     vy_interp = interp_vy(points)[0]
     vz_interp = interp_vz(points)[0]
 
-    #print(f"vx_interp = {vx_interp}, vy_interp = {vy_interp}, vz_interp = {vz_interp}")
-    
-    return vx_interp*ampl *cs, vy_interp*ampl *cs, vz_interp*ampl *cs
+    # print(f"vx_interp = {vx_interp}, vy_interp = {vy_interp}, vz_interp = {vz_interp}")
+
+    return vx_interp * ampl * cs, vy_interp * ampl * cs, vz_interp * ampl * cs
+
 
 # Example: Interpolate velocity at a specific position (in simulation coordinates)
 test_pos = (0.0, 0.0, 0.0)  # Center of simulation domain
@@ -546,28 +531,34 @@ print(f"  vx = {vel_x2:.6f}")
 print(f"  vy = {vel_y2:.6f}")
 print(f"  vz = {vel_z2:.6f}")
 
-def setup_particles():
 
+def setup_particles():
     model.set_field_value_lambda_f64_3("vxyz", vel_field)
 
     model.timestep()
 
 
-rho_plot_large = rho_xy_plot(model, (0.0, 0.0, 0.0), sim_radius *1.2, 2048, 2048, sim_directory, "rho_xy")
-rho_plot_mid = rho_xy_plot(model, (0.0, 0.0, 0.0), sim_radius /1.2, 2048, 2048, sim_directory, "rho_mid")
+rho_plot_large = rho_xy_plot(
+    model, (0.0, 0.0, 0.0), sim_radius * 1.2, 2048, 2048, sim_directory, "rho_xy"
+)
+rho_plot_mid = rho_xy_plot(
+    model, (0.0, 0.0, 0.0), sim_radius / 1.2, 2048, 2048, sim_directory, "rho_mid"
+)
 rho_plot_zoom = rho_xy_plot(model, (0.0, 0.0, 0.0), 0.25, 2048, 2048, sim_directory, "rho_zoom")
 
-v_slice_plot = v_slice_xy_plot(model, (0.0, 0.0, 0.0), sim_radius *1.2, 2048, 2048, sim_directory, "v_slice")
+v_slice_plot = v_slice_xy_plot(
+    model, (0.0, 0.0, 0.0), sim_radius * 1.2, 2048, 2048, sim_directory, "v_slice"
+)
+
 
 def analysis(i):
-    
     rho_plot_large.analysis_save(i)
     rho_plot_mid.analysis_save(i)
     v_slice_plot.analysis_save(i)
 
     dat = ctx.collect_data()
 
-    #get index with max rho
+    # get index with max rho
     max_rho_index = np.argmin(dat["hpart"])
     hpart_min = dat["hpart"][max_rho_index]
     max_rho = model.get_particle_mass() * (model.get_hfact() / hpart_min) ** 3
@@ -575,11 +566,13 @@ def analysis(i):
     max_rho_x = dat["xyz"][max_rho_index, 0]
     max_rho_y = dat["xyz"][max_rho_index, 1]
     max_rho_z = dat["xyz"][max_rho_index, 2]
-    print(f"max rho: {max_rho/kg_m3_codeu:.3e} at ({max_rho_x:.3e}, {max_rho_y:.3e}, {max_rho_z:.3e})")
+    print(
+        f"max rho: {max_rho / kg_m3_codeu:.3e} at ({max_rho_x:.3e}, {max_rho_y:.3e}, {max_rho_z:.3e})"
+    )
 
-    #render around max rho
+    # render around max rho
     rho_plot_zoom.center = (max_rho_x, max_rho_y, max_rho_z)
-    
+
     rho_plot_zoom.analysis_save(i)
 
 
@@ -598,7 +591,7 @@ next_t = 0
 if run_simulation == "1" and idump_last_dump is not None:
     model.load_from_dump(get_dump_name(idump_last_dump))
     i = idump_last_dump + 1
-    next_t =  model.get_time() + tff * 0.1
+    next_t = model.get_time() + tff * 0.1
 elif run_simulation == "1":
     config_model()
     setup_particles()
@@ -610,7 +603,6 @@ rho_plot_zoom.render_all()
 v_slice_plot.render_all()
 
 while i < 10000:
-
     model.evolve_until(next_t, niter_max=250)
 
     next_t = model.get_time() + tff * 1e-2
@@ -630,4 +622,4 @@ while i < 10000:
 rho_plot_large.render_all()
 rho_plot_mid.render_all()
 rho_plot_zoom.render_all()
-v_slice_plot.render_all()
\ No newline at end of file
+v_slice_plot.render_all()
diff --git a/src/shammath/src/Perlin.cpp b/src/shammath/src/Perlin.cpp
old mode 100755
new mode 100644
index bb336dab..00fd5f71
--- a/src/shammath/src/Perlin.cpp
+++ b/src/shammath/src/Perlin.cpp
@@ -5,284 +5,346 @@
  *      Author: tim
  */
 
-
-
-
 #include "Perlin.h"
 
 namespace noise {
 
-	int grad3[][3] = {{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0},
-			{1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},
-			{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}};
-
-	int grad4[][4]= {{0,1,1,1}, {0,1,1,-1}, {0,1,-1,1}, {0,1,-1,-1},
-			{0,-1,1,1}, {0,-1,1,-1}, {0,-1,-1,1}, {0,-1,-1,-1},
-			{1,0,1,1}, {1,0,1,-1}, {1,0,-1,1}, {1,0,-1,-1},
-			{-1,0,1,1}, {-1,0,1,-1}, {-1,0,-1,1}, {-1,0,-1,-1},
-			{1,1,0,1}, {1,1,0,-1}, {1,-1,0,1}, {1,-1,0,-1},
-			{-1,1,0,1}, {-1,1,0,-1}, {-1,-1,0,1}, {-1,-1,0,-1},
-			{1,1,1,0}, {1,1,-1,0}, {1,-1,1,0}, {1,-1,-1,0},
-			{-1,1,1,0}, {-1,1,-1,0}, {-1,-1,1,0}, {-1,-1,-1,0}};
-
-	int p[] = {151,160,137,91,90,15,
-			131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
-			190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
-			88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
-			77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
-			102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
-			135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
-			5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
-			223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
-			129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
-			251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
-			49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
-			138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180};
-
-	int* perm = new int[512];
-
-	void init(){
-		for(int i=0; i<512; i++){
-			perm[i]=p[i & 255];
-		}
-	}
+    int grad3[][3]
+        = {{1, 1, 0},
+           {-1, 1, 0},
+           {1, -1, 0},
+           {-1, -1, 0},
+           {1, 0, 1},
+           {-1, 0, 1},
+           {1, 0, -1},
+           {-1, 0, -1},
+           {0, 1, 1},
+           {0, -1, 1},
+           {0, 1, -1},
+           {0, -1, -1}};
 
-	int simplex[][4] = {
-			{0,1,2,3},{0,1,3,2},{0,0,0,0},{0,2,3,1},{0,0,0,0},{0,0,0,0},{0,0,0,0},{1,2,3,0},
-			{0,2,1,3},{0,0,0,0},{0,3,1,2},{0,3,2,1},{0,0,0,0},{0,0,0,0},{0,0,0,0},{1,3,2,0},
-			{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},
-			{1,2,0,3},{0,0,0,0},{1,3,0,2},{0,0,0,0},{0,0,0,0},{0,0,0,0},{2,3,0,1},{2,3,1,0},
-			{1,0,2,3},{1,0,3,2},{0,0,0,0},{0,0,0,0},{0,0,0,0},{2,0,3,1},{0,0,0,0},{2,1,3,0},
-			{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},
-			{2,0,1,3},{0,0,0,0},{0,0,0,0},{0,0,0,0},{3,0,1,2},{3,0,2,1},{0,0,0,0},{3,1,2,0},
-			{2,1,0,3},{0,0,0,0},{0,0,0,0},{0,0,0,0},{3,1,0,2},{0,0,0,0},{3,2,0,1},{3,2,1,0}};
+    int grad4[][4]
+        = {{0, 1, 1, 1},    {0, 1, 1, -1},  {0, 1, -1, 1},   {0, 1, -1, -1},  {0, -1, 1, 1},
+           {0, -1, 1, -1},  {0, -1, -1, 1}, {0, -1, -1, -1}, {1, 0, 1, 1},    {1, 0, 1, -1},
+           {1, 0, -1, 1},   {1, 0, -1, -1}, {-1, 0, 1, 1},   {-1, 0, 1, -1},  {-1, 0, -1, 1},
+           {-1, 0, -1, -1}, {1, 1, 0, 1},   {1, 1, 0, -1},   {1, -1, 0, 1},   {1, -1, 0, -1},
+           {-1, 1, 0, 1},   {-1, 1, 0, -1}, {-1, -1, 0, 1},  {-1, -1, 0, -1}, {1, 1, 1, 0},
+           {1, 1, -1, 0},   {1, -1, 1, 0},  {1, -1, -1, 0},  {-1, 1, 1, 0},   {-1, 1, -1, 0},
+           {-1, -1, 1, 0},  {-1, -1, -1, 0}};
 
+    int p[]
+        = {151, 160, 137, 91,  90,  15,  131, 13,  201, 95,  96,  53,  194, 233, 7,   225, 140, 36,
+           103, 30,  69,  142, 8,   99,  37,  240, 21,  10,  23,  190, 6,   148, 247, 120, 234, 75,
+           0,   26,  197, 62,  94,  252, 219, 203, 117, 35,  11,  32,  57,  177, 33,  88,  237, 149,
+           56,  87,  174, 20,  125, 136, 171, 168, 68,  175, 74,  165, 71,  134, 139, 48,  27,  166,
+           77,  146, 158, 231, 83,  111, 229, 122, 60,  211, 133, 230, 220, 105, 92,  41,  55,  46,
+           245, 40,  244, 102, 143, 54,  65,  25,  63,  161, 1,   216, 80,  73,  209, 76,  132, 187,
+           208, 89,  18,  169, 200, 196, 135, 130, 116, 188, 159, 86,  164, 100, 109, 198, 173, 186,
+           3,   64,  52,  217, 226, 250, 124, 123, 5,   202, 38,  147, 118, 126, 255, 82,  85,  212,
+           207, 206, 59,  227, 47,  16,  58,  17,  182, 189, 28,  42,  223, 183, 170, 213, 119, 248,
+           152, 2,   44,  154, 163, 70,  221, 153, 101, 155, 167, 43,  172, 9,   129, 22,  39,  253,
+           19,  98,  108, 110, 79,  113, 224, 232, 178, 185, 112, 104, 218, 246, 97,  228, 251, 34,
+           242, 193, 238, 210, 144, 12,  191, 179, 162, 241, 81,  51,  145, 235, 249, 14,  239, 107,
+           49,  192, 214, 31,  181, 199, 106, 157, 184, 84,  204, 176, 115, 121, 50,  45,  127, 4,
+           150, 254, 138, 236, 205, 93,  222, 114, 67,  29,  24,  72,  243, 141, 128, 195, 78,  66,
+           215, 61,  156, 180};
 
+    int *perm = new int[512];
 
+    void init() {
+        for (int i = 0; i < 512; i++) {
+            perm[i] = p[i & 255];
+        }
+    }
 
-	double noise(double xin, double yin){
-		double n0, n1, n2;
+    int simplex[][4]
+        = {{0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1}, {0, 0, 0, 0}, {0, 0, 0, 0},
+           {0, 0, 0, 0}, {1, 2, 3, 0}, {0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1},
+           {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
+           {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
+           {1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
+           {2, 3, 0, 1}, {2, 3, 1, 0}, {1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0},
+           {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
+           {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
+           {2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {3, 0, 1, 2}, {3, 0, 2, 1},
+           {0, 0, 0, 0}, {3, 1, 2, 0}, {2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
+           {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}};
 
-		double F2 = 0.5*(sqrt(3.0)-1.0);
-		double s = (xin+yin)*F2;
+    double noise(double xin, double yin) {
+        double n0, n1, n2;
 
-		int i = fastfloor(xin+s);
-		int j = fastfloor(yin+s);
-		double G2 = (3.0-sqrt(3.0))/6.0;
-		double t = (i+j)*G2;
-		double X0 = i-t;
+        double F2 = 0.5 * (sqrt(3.0) - 1.0);
+        double s  = (xin + yin) * F2;
 
-		double Y0 = j-t;
-		double x0 = xin-X0;
+        int i     = fastfloor(xin + s);
+        int j     = fastfloor(yin + s);
+        double G2 = (3.0 - sqrt(3.0)) / 6.0;
+        double t  = (i + j) * G2;
+        double X0 = i - t;
 
-		double y0 = yin-Y0;
-		int i1, j1;
-		if(x0>y0) {i1=1; j1=0;}
-		else {i1=0; j1=1;}
-		double x1 = x0 - i1 + G2;
+        double Y0 = j - t;
+        double x0 = xin - X0;
 
-		double y1 = y0 - j1 + G2;
-		double x2 = x0 - 1.0 + 2.0 * G2;
-		double y2 = y0 - 1.0 + 2.0 * G2;
-		int ii = i & 255;
-		int jj = j & 255;
-		int gi0 = perm[ii+perm[jj]] % 12;
-		int gi1 = perm[ii+i1+perm[jj+j1]] % 12;
-		int gi2 = perm[ii+1+perm[jj+1]] % 12;
-		double t0 = 0.5 - x0*x0-y0*y0;
-		if(t0<0) n0 = 0.0;
-		else {
-			t0 *= t0;
-			n0 = t0 * t0 * dot(grad3[gi0], x0, y0);
-		}
-		double t1 = 0.5 - x1*x1-y1*y1;
-		if(t1<0) n1 = 0.0;
-		else {
-			t1 *= t1;
-			n1 = t1 * t1 * dot(grad3[gi1], x1, y1);
-		}
+        double y0 = yin - Y0;
+        int i1, j1;
+        if (x0 > y0) {
+            i1 = 1;
+            j1 = 0;
+        } else {
+            i1 = 0;
+            j1 = 1;
+        }
+        double x1 = x0 - i1 + G2;
 
-		double t2 = 0.5 - x2*x2-y2*y2;
-		if(t2<0) n2 = 0.0;
-		else {
-			t2 *= t2;
-			n2 = t2 * t2 * dot(grad3[gi2], x2, y2);
-		}
+        double y1 = y0 - j1 + G2;
+        double x2 = x0 - 1.0 + 2.0 * G2;
+        double y2 = y0 - 1.0 + 2.0 * G2;
+        int ii    = i & 255;
+        int jj    = j & 255;
+        int gi0   = perm[ii + perm[jj]] % 12;
+        int gi1   = perm[ii + i1 + perm[jj + j1]] % 12;
+        int gi2   = perm[ii + 1 + perm[jj + 1]] % 12;
+        double t0 = 0.5 - x0 * x0 - y0 * y0;
+        if (t0 < 0)
+            n0 = 0.0;
+        else {
+            t0 *= t0;
+            n0 = t0 * t0 * dot(grad3[gi0], x0, y0);
+        }
+        double t1 = 0.5 - x1 * x1 - y1 * y1;
+        if (t1 < 0)
+            n1 = 0.0;
+        else {
+            t1 *= t1;
+            n1 = t1 * t1 * dot(grad3[gi1], x1, y1);
+        }
 
-		return 70.0 * (n0 + n1 + n2);
-	}
+        double t2 = 0.5 - x2 * x2 - y2 * y2;
+        if (t2 < 0)
+            n2 = 0.0;
+        else {
+            t2 *= t2;
+            n2 = t2 * t2 * dot(grad3[gi2], x2, y2);
+        }
 
+        return 70.0 * (n0 + n1 + n2);
+    }
 
-	double noise(double xin, double yin, double zin){
-		double n0, n1, n2, n3;
-		double F3 = 1.0/3.0;
-		double s = (xin+yin+zin)*F3;
-		int i = fastfloor(xin+s);
-		int j = fastfloor(yin+s);
-		int k = fastfloor(zin+s);
-		double G3 = 1.0/6.0;
-		double t = (i+j+k)*G3;
-		double X0 = i-t;
-		double Y0 = j-t;
-		double Z0 = k-t;
-		double x0 = xin-X0;
-		double y0 = yin-Y0;
-		double z0 = zin-Z0;
-		int i1, j1, k1;
-		int i2, j2, k2;
-		if(x0>=y0) {
-			if(y0>=z0)
-			{ i1=1; j1=0; k1=0; i2=1; j2=1; k2=0; }
-			else if(x0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=0; k2=1; }
-			else { i1=0; j1=0; k1=1; i2=1; j2=0; k2=1; }
-		}
-		else {
-			if(y0<z0) { i1=0; j1=0; k1=1; i2=0; j2=1; k2=1; }
-			else if(x0<z0) { i1=0; j1=1; k1=0; i2=0; j2=1; k2=1; }
-			else { i1=0; j1=1; k1=0; i2=1; j2=1; k2=0; }
-		}
+    double noise(double xin, double yin, double zin) {
+        double n0, n1, n2, n3;
+        double F3 = 1.0 / 3.0;
+        double s  = (xin + yin + zin) * F3;
+        int i     = fastfloor(xin + s);
+        int j     = fastfloor(yin + s);
+        int k     = fastfloor(zin + s);
+        double G3 = 1.0 / 6.0;
+        double t  = (i + j + k) * G3;
+        double X0 = i - t;
+        double Y0 = j - t;
+        double Z0 = k - t;
+        double x0 = xin - X0;
+        double y0 = yin - Y0;
+        double z0 = zin - Z0;
+        int i1, j1, k1;
+        int i2, j2, k2;
+        if (x0 >= y0) {
+            if (y0 >= z0) {
+                i1 = 1;
+                j1 = 0;
+                k1 = 0;
+                i2 = 1;
+                j2 = 1;
+                k2 = 0;
+            } else if (x0 >= z0) {
+                i1 = 1;
+                j1 = 0;
+                k1 = 0;
+                i2 = 1;
+                j2 = 0;
+                k2 = 1;
+            } else {
+                i1 = 0;
+                j1 = 0;
+                k1 = 1;
+                i2 = 1;
+                j2 = 0;
+                k2 = 1;
+            }
+        } else {
+            if (y0 < z0) {
+                i1 = 0;
+                j1 = 0;
+                k1 = 1;
+                i2 = 0;
+                j2 = 1;
+                k2 = 1;
+            } else if (x0 < z0) {
+                i1 = 0;
+                j1 = 1;
+                k1 = 0;
+                i2 = 0;
+                j2 = 1;
+                k2 = 1;
+            } else {
+                i1 = 0;
+                j1 = 1;
+                k1 = 0;
+                i2 = 1;
+                j2 = 1;
+                k2 = 0;
+            }
+        }
 
-		double x1 = x0 - i1 + G3;
-		double y1 = y0 - j1 + G3;
-		double z1 = z0 - k1 + G3;
-		double x2 = x0 - i2 + 2.0*G3;
-		double y2 = y0 - j2 + 2.0*G3;
-		double z2 = z0 - k2 + 2.0*G3;
-		double x3 = x0 - 1.0 + 3.0*G3;
-		double y3 = y0 - 1.0 + 3.0*G3;
-		double z3 = z0 - 1.0 + 3.0*G3;
-		int ii = i & 255;
-		int jj = j & 255;
-		int kk = k & 255;
-		int gi0 = perm[ii+perm[jj+perm[kk]]] % 12;
-		int gi1 = perm[ii+i1+perm[jj+j1+perm[kk+k1]]] % 12;
-		int gi2 = perm[ii+i2+perm[jj+j2+perm[kk+k2]]] % 12;
-		int gi3 = perm[ii+1+perm[jj+1+perm[kk+1]]] % 12;
-		double t0 = 0.6 - x0*x0 - y0*y0 - z0*z0;
-		if(t0<0) n0 = 0.0;
-		else {
-			t0 *= t0;
-			n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0);
-		}
-		double t1 = 0.6 - x1*x1 - y1*y1 - z1*z1;
-		if(t1<0) n1 = 0.0;
-		else {
-			t1 *= t1;
-			n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1);
-		}
-		double t2 = 0.6 - x2*x2 - y2*y2 - z2*z2;
-		if(t2<0) n2 = 0.0;
-		else {
-			t2 *= t2;
-			n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2);
-		}
-		double t3 = 0.6 - x3*x3 - y3*y3 - z3*z3;
-		if
 ... truncated ...