diff --git a/Project.toml b/Project.toml
index ec8954c..f2bed14 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,7 +10,6 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 PaddedViews = "5432bcbf-9aad-5242-b902-cca2824c8663"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-TensorCast = "02d47bb6-7ce6-556a-be16-bb1710789e2b"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
@@ -20,7 +19,6 @@ JLD2 = "0.4"
 PaddedViews = "0.5"
 QuadGK = "2"
 SpecialFunctions = "1, 1.3, 1.4, 2"
-TensorCast = "0.4"
 Tullio = "0.2, 0.3"
 UnPack = "1, 1.0"
 julia = "1.8, 1.9, 1.10"
diff --git a/src/QuantumFluidSpectra.jl b/src/QuantumFluidSpectra.jl
index 106b061..72ba49a 100644
--- a/src/QuantumFluidSpectra.jl
+++ b/src/QuantumFluidSpectra.jl
@@ -3,11 +3,10 @@ module QuantumFluidSpectra
 using Tullio
 using Hwloc
 using FFTW
-FFTW.set_num_threads(num_physical_cores())
 using SpecialFunctions
 using PaddedViews
 using UnPack
-using TensorCast
+FFTW.set_num_threads(8)
 
 # fallback since fast_hypot is 2 argument only
 @fastmath hypot(x::Float64, y::Float64, z::Float64)=sqrt(x^2+y^2+z^2)
@@ -31,8 +30,12 @@ export log10range, convolve
 export xk_arrays, fft_differentials
 export gradient, velocity, current
 export energydecomp, helmholtz, kinetic_density
-export incompressible_spectrum, compressible_spectrum, qpressure_spectrum
-export incompressible_density, compressible_density, qpressure_density
+export incompressible_spectrum, 
+compressible_spectrum, 
+qpressure_spectrum
+export incompressible_density, 
+compressible_density, 
+qpressure_density
 export ic_density, iq_density, cq_density
 export density_spectrum, trap_spectrum
 
diff --git a/src/analysis.jl b/src/analysis.jl
index 92b75a3..5978df5 100644
--- a/src/analysis.jl
+++ b/src/analysis.jl
@@ -107,9 +107,9 @@ Returned fields `Wi`, `Wc` are tuples of cartesian components of incompressible
 """
 function helmholtz(wx, wy, kx, ky)
     wxk = fft(wx); wyk = fft(wy)
-    @cast kw[i,j] := (kx[i] * wxk[i,j] + ky[j] * wyk[i,j])/ (kx[i]^2+ky[j]^2)
-    @cast wxkc[i,j] := kw[i,j] * kx[i] 
-    @cast wykc[i,j] := kw[i,j] * ky[j]
+    kw = @. (kx * wxk + ky' * wyk)/ (kx^2+ky'^2)
+    wxkc = @. kw*kx
+    wykc = @. kw*ky'
     wxkc[1] = zero(wxkc[1]); wykc[1] = zero(wykc[1])
     wxki = @. wxk - wxkc
     wyki = @. wyk - wykc
@@ -121,10 +121,11 @@ end
 
 function helmholtz(wx, wy, wz, kx, ky, kz)
     wxk = fft(wx); wyk = fft(wy); wzk = fft(wz)
-    @cast kw[i,j,k] := (kx[i] * wxk[i,j,k] + ky[j] * wyk[i,j,k] + kz[k] * wzk[i,j,k])/ (kx[i]^2 + ky[j]^2 + kz[k]^2)
-    @cast wxkc[i,j,k] := kw[i,j,k] * kx[i]  
-    @cast wykc[i,j,k] := kw[i,j,k] * ky[j] 
-    @cast wzkc[i,j,k] := kw[i,j,k] * kz[k]  
+    kzr = reshape(kz,1,1,length(kz))
+    kw = @. (kx * wxk + ky' * wyk + kzr * wzk)/ (kx^2 + ky'^2 + kzr^2)
+    wxkc = @. kw * kx  
+    wykc = @. kw * ky'
+    wzkc = @. kw * kzr  
     wxkc[1] = zero(wxkc[1]); wykc[1] = zero(wykc[1]); wzkc[1] = zero(wzkc[1])
     wxki = @. wxk - wxkc
     wyki = @. wyk - wykc