diff --git a/.gitignore b/.gitignore
index 74ddcb5..fa669d1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -25,8 +25,12 @@ whislerwavelin/
 whislerwaveHLL/
 hallharrisres/
 hallharrisres2/
+hallharrisresCT/
+hallharrisresCT2/
+hallharrisresCT2long/
 subprojects
 
 debug/
 frames/
+framestest/
 pyMHD/pyMHD.cpython*
diff --git a/README.md b/README.md
index e875a8d..9880363 100644
--- a/README.md
+++ b/README.md
@@ -1,12 +1,5 @@
 ## TODO
-
-# By priority order in short term
-- Implement diagnostics and testing (the rough debugging is done).
-- Implement post-processing
-- Make the code more modular
-
-# And then
-- Extend the model (3D, Non-Ideal MHD, more sofisticated Riemann solvers, more sophisticated CT ...).
+- Extend the model (3D, more sofisticated Riemann solvers, more sophisticated CT ...).
 - Clean up the code.
 - Integrate in PHARE.
 - Add coupling (Fluid-Fluid, Fluid-Hybrid)
\ No newline at end of file
diff --git a/diagnostics/Hallharris/colormesh.py b/diagnostics/Hallharris/colormesh.py
index e696b92..597d2ba 100644
--- a/diagnostics/Hallharris/colormesh.py
+++ b/diagnostics/Hallharris/colormesh.py
@@ -12,7 +12,7 @@ def reshape_data(data, ny, nx):
 nx = 250
 ny = 250
 
-file_path = "hallharrisres/PRK2_15_0001697500.txt"
+file_path = "hallharrisresCT/PRK2_0_5002794389.txt"
 
 data = read_data(file_path)
 
@@ -22,8 +22,10 @@ def reshape_data(data, ny, nx):
 dy = 0.1
 
 # Extract Bx and By
+rho = reshaped_data[:, :, 0]
 Bx = reshaped_data[:, :, 4]
 By = reshaped_data[:, :, 5]
+Bz = reshaped_data[:, :, 6]
 
 # Calculate the derivatives
 dBy_dx = np.gradient(By, dx, axis=0)
@@ -32,24 +34,26 @@ def reshape_data(data, ny, nx):
 # Calculate Jz
 Jz = (dBy_dx - dBx_dy)
 
-plt.pcolormesh(Jz.T, cmap='coolwarm')  
-plt.title('Contour Plot of Jz at t=30')
+toPlot = Bz
+
+plt.pcolormesh(toPlot.T, cmap='coolwarm')  
+plt.title('Contour Plot of Bz at t=0.5')
 plt.xlabel('X')
 plt.ylabel('Y')
 plt.show()
 
-x_middle = nx // 2
+# x_middle = nx // 2
 
-Jz_cutx = Jz[x_middle, :]
+# Jz_cutx = Jz[x_middle, :]
 
-y_positions = np.arange(ny)
+# y_positions = np.arange(ny)
 
-# Create the plot
-plt.figure(figsize=(10, 6))
-plt.plot(y_positions, Jz_cutx, marker='o')
-plt.title(f'1D cut of Jz at X = {x_middle/nx}')
-plt.xlabel('y')
-plt.ylabel('Jz')
-plt.grid(True)
+# # Create the plot
+# plt.figure(figsize=(10, 6))
+# plt.plot(y_positions, Jz_cutx, marker='o')
+# plt.title(f'1D cut of Jz at X = {x_middle/nx}')
+# plt.xlabel('y')
+# plt.ylabel('Jz')
+# plt.grid(True)
 
-plt.show()
+# plt.show()
diff --git a/diagnostics/Hallharris/colormesh2.py b/diagnostics/Hallharris/colormesh2.py
new file mode 100644
index 0000000..6725e9e
--- /dev/null
+++ b/diagnostics/Hallharris/colormesh2.py
@@ -0,0 +1,52 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def read_data(file_path):
+    data = np.loadtxt(file_path)
+    return data
+
+def reshape_data(data, ny, nx):
+    reshaped_data = data.reshape((ny, nx, -1), order='F')
+    return reshaped_data
+
+nx = 250
+ny = 250
+
+file_path1 = "hallharrisresCT/PRK2_0_5002794389.txt"
+file_path2 = "hallharrisresCT2/PRK2_0_5000859584.txt"
+# file_path1 = "hallharrisres2/PRK2_2_4174122015.txt"
+# file_path2 = "hallharrisresCT2long/PRK2_2_4039801735.txt"
+
+# Read and process data from the first file
+data1 = read_data(file_path1)
+reshaped_data1 = reshape_data(data1, ny, nx)
+
+# Extract Bz for the first file
+Bz1 = reshaped_data1[:, :, 6]
+
+# Read and process data from the second file
+data2 = read_data(file_path2)
+reshaped_data2 = reshape_data(data2, ny, nx)
+
+# Extract Bz for the second file
+Bz2 = reshaped_data2[:, :, 6]
+
+# Create subplots
+fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+
+# Plot Bz for the first file
+im1 = axes[0].pcolormesh(Bz1.T, cmap='coolwarm')
+axes[0].set_title('Contour Plot of Bz at t=0.5 (NoCT)')
+axes[0].set_xlabel('X')
+axes[0].set_ylabel('Y')
+fig.colorbar(im1, ax=axes[0])
+
+# Plot Bz for the second file
+im2 = axes[1].pcolormesh(Bz2.T, cmap='coolwarm')
+axes[1].set_title('Contour Plot of Bz at t=0.5 (CT)')
+axes[1].set_xlabel('X')
+axes[1].set_ylabel('Y')
+fig.colorbar(im2, ax=axes[1])
+
+plt.tight_layout()
+plt.show()
diff --git a/diagnostics/Hallharris/colormesh_anim.py b/diagnostics/Hallharris/colormesh_anim.py
index 5ef654d..86b30df 100644
--- a/diagnostics/Hallharris/colormesh_anim.py
+++ b/diagnostics/Hallharris/colormesh_anim.py
@@ -24,7 +24,7 @@ def read_times(file_paths):
         times.append(time)
     return times
 
-results_dir = "hallharrisres/"
+results_dir = "hallharrisresCT2long/"
 file_paths = [results_dir + file for file in os.listdir(results_dir) if file.startswith("PRK2_") and file.endswith(".txt")]
 
 nx = 250
@@ -41,8 +41,10 @@ def read_times(file_paths):
 dy = 0.1
 
 # Extract Bx and By
+rho = [reshaped_data[i][:, :, 0] for i in range(len(data))]
 Bx = [reshaped_data[i][:, :, 4] for i in range(len(data))]
 By = [reshaped_data[i][:, :, 5] for i in range(len(data))]
+Bz = [reshaped_data[i][:, :, 6] for i in range(len(data))]
 
 # Calculate the derivatives
 dBy_dx = [np.gradient(By[i], dx, axis=0) for i in range(len(data))]
@@ -51,7 +53,7 @@ def read_times(file_paths):
 # Calculate Jz
 Jz = [(dBy_dx[i] - dBx_dy[i]) for i in range(len(data))]
 
-toPlot = Jz
+toPlot = Bz
 
 data_min = np.min(toPlot[-1])
 data_max = np.max(toPlot[-1])
@@ -70,7 +72,7 @@ def read_times(file_paths):
 
 def update(frame):
     im.set_array(toPlot[frame].T)
-    plt.title(f'Jz at t={times[frame]}')
+    plt.title(f'rho at t={times[frame]}')
     plt.savefig(f'{output_dir}/frame_{frame:04d}.png')
     return im,
 
@@ -83,4 +85,4 @@ def update(frame):
 plt.ylabel('Y')
 plt.show()
 
-#ffmpeg -r 10 -i frames/frame_%04d.png -vcodec mpeg4 -q:v 5 orszagtangP.mp4
\ No newline at end of file
+#ffmpeg -r 10 -i frames/frame_%04d.png -vcodec mpeg4 -q:v 5 hallharris.mp4
\ No newline at end of file
diff --git a/diagnostics/Hallharris/hallharris.py b/diagnostics/Hallharris/hallharris.py
index 4e45adf..d1fa9c9 100644
--- a/diagnostics/Hallharris/hallharris.py
+++ b/diagnostics/Hallharris/hallharris.py
@@ -13,7 +13,7 @@
 Dx = 0.1
 Dy = 0.1
 Dt = 0.0
-FinalTime = 30
+FinalTime = 0.5
 nghost = 2
 
 boundaryconditions = p.BoundaryConditions.Periodic
@@ -104,7 +104,7 @@ def P_(x, y):
 
 #############################################################################################################################################################################
 
-result_dir = 'hallharrisres2/'
+result_dir = 'hallharrisresCT2/'
 if os.path.exists(result_dir):
     shutil.rmtree(result_dir)
 
diff --git a/diagnostics/MHDrotor/MHDrotor.py b/diagnostics/MHDrotor/MHDrotor.py
index 4ae0883..9021ee9 100644
--- a/diagnostics/MHDrotor/MHDrotor.py
+++ b/diagnostics/MHDrotor/MHDrotor.py
@@ -28,12 +28,11 @@
 
 ##############################################################################################################################################################################
 B0 = 5/(np.sqrt(4*np.pi))
-v0 = 1
+v0 = 2
 
 r0 = 0.1
 r1 = 0.115
 
-# Domain [-0.5, 0.5]*[-0.5, 0.5]
 def r(x, y):
     return np.sqrt((x-0.5)**2 + (y-0.5)**2)
 
@@ -51,14 +50,14 @@ def vx_(x, y):
     r_ = r(x, y)
     f_ = f(r_)
     
-    vx_values = np.where(r_ <= r0, -f_ * v0 * (y - 0.5) / r0, np.where(r_ < r1, -f_ * v0 * (y - 0.5) / r_, 0.0))
+    vx_values = np.where(r_ <= r0, - v0 * (y - 0.5) / r0, np.where(r_ < r1, -f_ * v0 * (y - 0.5) / r_, 0.0))
     return vx_values
 
 def vy_(x, y):
     r_ = r(x, y)
     f_ = f(r_)
     
-    vy_values = np.where(r_ <= r0, f_ * v0 * (x - 0.5) / r0, np.where(r_ < r1, f_ * v0 * (x - 0.5) / r_, 0.0))
+    vy_values = np.where(r_ <= r0, v0 * (x - 0.5) / r0, np.where(r_ < r1, f_ * v0 * (x - 0.5) / r_, 0.0))
     return vy_values
 
 def vz_(x, y):
@@ -74,7 +73,7 @@ def Bz_(x, y):
     return 0.0 
 
 def P_(x, y):
-    return 1.0
+    return 0.5
 
 
 x = np.arange(nx) * Dx + 0.5 * Dx
diff --git a/diagnostics/MHDrotor/colormesh.py b/diagnostics/MHDrotor/colormesh.py
new file mode 100644
index 0000000..a610207
--- /dev/null
+++ b/diagnostics/MHDrotor/colormesh.py
@@ -0,0 +1,27 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def read_data(file_path):
+    data = np.loadtxt(file_path)
+    return data
+
+def reshape_data(data, ny, nx):
+    reshaped_data = data.reshape((ny, nx, -1), order='F')
+    return reshaped_data
+
+nx = 100
+ny = 100
+
+file_path = "MHDrotorres/URK2_0_0.txt"
+
+data = read_data(file_path)
+
+reshaped_data = reshape_data(data, nx, ny)
+
+rho = reshaped_data[:, :, 0]
+
+plt.pcolormesh(rho.T, cmap='coolwarm')  
+plt.title('Contour Plot of rho at t=0')
+plt.xlabel('X')
+plt.ylabel('Y')
+plt.show()
\ No newline at end of file
diff --git a/diagnostics/MHDrotor/colormesh_anim.py b/diagnostics/MHDrotor/colormesh_anim.py
index 35c475f..5c995ea 100644
--- a/diagnostics/MHDrotor/colormesh_anim.py
+++ b/diagnostics/MHDrotor/colormesh_anim.py
@@ -25,8 +25,8 @@ def read_times(file_paths):
 results_dir = "MHDrotorres/"
 file_paths = [results_dir + file for file in os.listdir(results_dir) if file.startswith("URK2_") and file.endswith(".txt")]
 
-nx = 200
-ny = 200
+nx = 100
+ny = 100
 
 data = [read_data(file_path) for file_path in file_paths]
 times = read_times(file_paths)
diff --git a/diagnostics/orszag-tang/colormesh_anim.py b/diagnostics/orszag-tang/colormesh_anim.py
index 4cfa429..a707cae 100644
--- a/diagnostics/orszag-tang/colormesh_anim.py
+++ b/diagnostics/orszag-tang/colormesh_anim.py
@@ -2,6 +2,7 @@
 import matplotlib.pyplot as plt
 import os
 from matplotlib.colors import Normalize
+from matplotlib.animation import FuncAnimation
 import shutil
 
 # Function to read data from file
@@ -59,8 +60,7 @@ def update(frame):
     plt.savefig(f'{output_dir}/frame_{frame:04d}.png')
     return im,
 
-for frame in range(len(times)):
-    update(frame)
+ani = FuncAnimation(fig, update, frames=len(times), interval=100)
 
 plt.xlabel('X')
 plt.ylabel('Y')
diff --git a/diagnostics/orszag-tang/orszagtang.py b/diagnostics/orszag-tang/orszagtang.py
index d45a8ef..b68c1e2 100644
--- a/diagnostics/orszag-tang/orszagtang.py
+++ b/diagnostics/orszag-tang/orszagtang.py
@@ -8,13 +8,13 @@
 
 #############################################################################################################################################################################
 
-nx = 128
-ny = 128
+nx = 1024
+ny = 1024
 Dx = 1/nx
 Dy = 1/ny
 Dt = 0.0
-FinalTime = 0.5
-nghost = 1
+FinalTime = 1
+nghost = 2
 
 boundaryconditions = p.BoundaryConditions.Periodic
 
@@ -22,12 +22,12 @@
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
 constainedtransport = p.CTMethod.Average
-timeintegrator = p.Integrator.EulerIntegrator
+timeintegrator = p.Integrator.TVDRK2Integrator
 
 dumpvariables = p.dumpVariables.Primitive
 dumpfrequency = 80
 
-OptionalPhysics = p.OptionalPhysics.HallResHyper
+OptionalPhysics = p.OptionalPhysics.Off
 
 ##############################################################################################################################################################################
 B0 = 1./(np.sqrt(4.*np.pi))
diff --git a/diagnostics/whistlerwave/fft.py b/diagnostics/whistlerwave/fft.py
index 199cdc9..88c1cb6 100644
--- a/diagnostics/whistlerwave/fft.py
+++ b/diagnostics/whistlerwave/fft.py
@@ -8,8 +8,13 @@
 nx = 128
 Dx = 0.1
 
-Dt = 0.0025
-finalTime = 1
+Dt = 0.002
+
+lx=nx*Dx
+m = 4
+kt=2*np.pi/lx * m
+w = (kt**2 /2) *(np.sqrt(1+4/kt**2) + 1)
+finalTime = 2*np.pi / w
 
 column_names = ['rho', 'vx', 'vy', 'vz', 'Bx', 'By', 'Bz', 'P']
 
@@ -54,21 +59,17 @@ def leftAlfvenIonCyclotron(k):
 kplus = k[:halfk]
 wplus = w[:halfw]
 
-"""modes = [4]
-kmodes = 2*np.pi*np.asarray([1/(nx*Dx) * m for m in modes])"""
+modes = [m]
+kmodes = 2*np.pi*np.asarray([1/(nx*Dx) * m for m in modes])
 
-"""
-fig,ax = plt.subplots()
-ax.plot(kplus, whistler(kplus))
-ax.plot(kplus, np.abs(fft(quantities['By'][-1,:] + 1j*fft(quantities['Bz'][-1,:]))[:halfk]))
-for km in kmodes:
-    ax.axvline(km, ls='--', color='k')
-    ax.plot(km, km**2, marker='o')
+B_combined = quantities['By'] + 1j * quantities['Bz']
 
-plt.show()
-"""
+# window_space = np.hanning(nx)
+# window_time = np.hanning(int(finalTime / Dt) + 1)
+# window = np.outer(window_time, window_space)
+
+# B_combined *= window
 
-B_combined = quantities['By'] + 1j * quantities['Bz']
 B_fft = fft2(B_combined)
 B_fft_abs = np.abs(B_fft)[:halfw,:halfk]
 
@@ -80,8 +81,8 @@ def leftAlfvenIonCyclotron(k):
 ax.plot(kplus, whistler(kplus), marker = '+', color='k')
 ax.plot(kplus, leftAlfvenIonCyclotron(kplus), marker = '*', color='g')
 ax.plot(kplus, kplus)
-"""for km in kmodes:
+for km in kmodes:
     ax.axvline(km, ls='--', color='k')
-    ax.plot(km, km**2, marker='o')"""
+    ax.plot(km, whistler(km), marker='o')
 
 plt.show()
\ No newline at end of file
diff --git a/diagnostics/whistlerwave/plot_anim.py b/diagnostics/whistlerwave/plot_anim.py
index fb6adfc..41521d9 100644
--- a/diagnostics/whistlerwave/plot_anim.py
+++ b/diagnostics/whistlerwave/plot_anim.py
@@ -136,3 +136,29 @@ def onclick(event):
 plt.axvline(x[-1]-Dx/2, ls='--', color='k')
 plt.show()
 """
+def calculate_total_energy(quantities, t, gamma):
+    rho = quantities['rho'][t, fixed_index, :]
+    vx = quantities['vx'][t, fixed_index, :]
+    vy = quantities['vy'][t, fixed_index, :]
+    vz = quantities['vz'][t, fixed_index, :]
+    Bx = quantities['Bx'][t, fixed_index, :]
+    By = quantities['By'][t, fixed_index, :]
+    Bz = quantities['Bz'][t, fixed_index, :]
+    P = quantities['P'][t, fixed_index, :]
+
+    # Calculate kinetic energy density
+    kinetic_energy = 0.5 * rho * (vx**2 + vy**2 + vz**2)
+    
+    # Calculate magnetic energy density
+    magnetic_energy = 0.5 * (Bx**2 + By**2 + Bz**2)
+    
+    # Calculate internal energy density using pressure and gamma
+    internal_energy = P / (gamma - 1)
+    
+    # Total energy density
+    total_energy = internal_energy + kinetic_energy + magnetic_energy
+    
+    return np.sum(total_energy)
+
+print(calculate_total_energy(quantities, 0, 5.0/3.0))
+print(calculate_total_energy(quantities, -1, 5.0/3.0))
diff --git a/diagnostics/whistlerwave/whislerwaveX.py b/diagnostics/whistlerwave/whislerwaveX.py
index de8aa36..4937cc6 100644
--- a/diagnostics/whistlerwave/whislerwaveX.py
+++ b/diagnostics/whistlerwave/whislerwaveX.py
@@ -12,8 +12,8 @@
 ny = 1
 Dx = 0.1
 Dy = 1
-Dt = 0.0025
-FinalTime = 1
+Dt = 0.002
+
 nghost = 2
 
 boundaryconditions = p.BoundaryConditions.Periodic
@@ -34,11 +34,16 @@
 ##############################################################################################################################################################################
 lx=nx*Dx
 k=2*np.pi/lx
-cw = (k/2) + (np.sqrt(1+k**2/4))
+
+m = 4
+
+kt=2*np.pi/lx * m
+w = (kt**2 /2) *(np.sqrt(1+4/kt**2) + 1)
+FinalTime = 2*np.pi / w
 
 np.random.seed(0)
 
-modes = [4]#[int(nx/4)]
+modes = [m]#[int(nx/4)]
 phases = np.random.rand(len(modes))
 
 def rho_(x, y):
diff --git a/pyMHD/pyMHD.cpp b/pyMHD/pyMHD.cpp
index 54c9524..478a856 100644
--- a/pyMHD/pyMHD.cpp
+++ b/pyMHD/pyMHD.cpp
@@ -108,7 +108,8 @@ PYBIND11_MODULE(pyMHD, m)
     py::enum_<CTMethod>(m, "CTMethod")
         .value("Average", Average)
         .value("Contact", Contact)
-        .value("UCT_HLL", UCT_HLL);
+        .value("UCT_HLL", UCT_HLL)
+        .value("NoCT", NoCT);
 
     py::enum_<Integrator>(m, "Integrator")
         .value("EulerIntegrator", EulerIntegrator)
diff --git a/src/ComputeJ.hpp b/src/ComputeJ.hpp
index 8cb0791..aaf5079 100644
--- a/src/ComputeJ.hpp
+++ b/src/ComputeJ.hpp
@@ -3,7 +3,7 @@
 
 // J has 1 more ghost cells for laplacian calculations in Flux vector
 template<typename Variables>
-void ComputeJ(Variables& Vn, double Dx, double Dy, int nghost) {
+void ComputeJ(Variables& Vn, double Dx, double Dy, int nghost, CTMethod ct) {
 
     for (int j = nghost; j <= Vn.ny - nghost; ++j)
     {
@@ -26,8 +26,10 @@ void ComputeJ(Variables& Vn, double Dx, double Dy, int nghost) {
         for (int i = nghost; i <= Vn.nx - nghost; ++i)
         {
             // In 1D (no CT), use cell centered Bx and By
-            Vn.Jz[j + 1][i + 1] =  (Vn.By[j][i] - Vn.By[j][i - 1]) / (Dx) - (Vn.Bx[j][i] - Vn.Bx[j - 1][i]) / (Dy);
-            // Vn.Jz[j + 1][i + 1] = (Vn.Byf[j][i] - Vn.Byf[j][i - 1]) / (Dx) - (Vn.Bxf[j][i] - Vn.Bxf[j - 1][i]) / (Dy);
+            if (ct == CTMethod::NoCT)
+                Vn.Jz[j + 1][i + 1] =  (Vn.By[j][i] - Vn.By[j][i - 1]) / (Dx) - (Vn.Bx[j][i] - Vn.Bx[j - 1][i]) / (Dy);
+            else
+                Vn.Jz[j + 1][i + 1] = (Vn.Byf[j][i] - Vn.Byf[j][i - 1]) / (Dx) - (Vn.Bxf[j][i] - Vn.Bxf[j - 1][i]) / (Dy);
         }
     }
 }
diff --git a/src/Enums.hpp b/src/Enums.hpp
index a30f444..f52e083 100644
--- a/src/Enums.hpp
+++ b/src/Enums.hpp
@@ -6,7 +6,7 @@ enum BoundaryConditions { Periodic, ZeroGradient };
 enum Reconstruction { Constant, Linear };
 enum Slope { VanLeer, MinMod };
 enum Riemann { Rusanov, HLL };
-enum CTMethod { Average, Contact, UCT_HLL };
+enum CTMethod { Average, Contact, UCT_HLL, NoCT };
 enum Integrator { EulerIntegrator, TVDRK2Integrator, TVDRK3Integrator };
 
 enum OptionalPhysics { Hall, Res, Hyper, HallRes, HallHyper, ResHyper, HallResHyper, Off };
diff --git a/src/TimeIntegrator.cpp b/src/TimeIntegrator.cpp
index a406e7d..efdda6f 100644
--- a/src/TimeIntegrator.cpp
+++ b/src/TimeIntegrator.cpp
@@ -18,7 +18,7 @@ ConservativeVariables Euler(ConservativeVariables& Un, double Dx, double Dy, dou
     UpdateGhostCells(Un, nghost, bc);
 
     if (OptP == OptionalPhysics::HallResHyper) {
-        ComputeJ(Un, Dx, Dy, nghost);
+        ComputeJ(Un, Dx, Dy, nghost, ct);
         UpdateGhostJ(Un, nghost, bc);
         
         // static int j = 0;
@@ -36,9 +36,12 @@ ConservativeVariables Euler(ConservativeVariables& Un, double Dx, double Dy, dou
 
     Un1 = EulerAdvance(Un, Dx, Dy, Dt, nghost, rec, sl, rs, OptP);
 
-    //ApplyConstrainedTransport(Un1, Un, Dx, Dy, Dt, nghost, rec, sl, rs, ct, OptP); // If Un1 not needed in CT scheme (else update ghost cells before)
+    // Desactivate in 1D
+    if (ct != CTMethod::NoCT) {
+        ApplyConstrainedTransport(Un1, Un, Dx, Dy, Dt, nghost, rec, sl, rs, ct, OptP); // If Un1 not needed in CT scheme (else update ghost cells before)
     
-    //Un1.ReconstructCenteredB(nghost);
+        Un1.ReconstructCenteredB(nghost);
+    }
 
     return Un1;
 }