diff --git a/.github/workflows/continuousintegration.yaml b/.github/workflows/continuousintegration.yaml
index 8b183cb..aff85c5 100644
--- a/.github/workflows/continuousintegration.yaml
+++ b/.github/workflows/continuousintegration.yaml
@@ -22,6 +22,7 @@ jobs:
         run: |  
           python -m pip install --upgrade pip  
           pip install -r test_requirements.txt
+          python tests/importtest.py
 #      - name: Test with pytest  
 #        run: |  
 #         coverage run -m pytest  -v -s  
diff --git a/README.md b/README.md
index 13352a7..e81b626 100644
--- a/README.md
+++ b/README.md
@@ -1,13 +1,29 @@
 # exptool: galactic dynamics in python using basis function expansions
 ## msp, 2014-present
 
-### Version 0.3
-exptool is primarily collection of analysis methods used to analyze _n_-body simulations, primarily those running with _EXP_, though any _n_-body
-code could use the methodology (Generation 1).
+### Version 0.4
+exptool is primarily collection of analysis methods used to analyze _n_-body simulations, primarily those running with `EXP`, though any _n_-body
+code could use the methodology (Generation 0.1).
 
-As of 2018, exptool has expanded to be an all-purpose dynamical framework for studying disc galaxies (Generation 2).
+As of 2018, exptool has expanded to be an all-purpose dynamical framework for studying disc galaxies (Generation 0.2).
+
+As of 2020, exptool includes support for analysing stellar halos (Generation 0.3).
+
+As of 2024, `exptool` supports analysis of cosmological simulation outputs for bar measurements and a range of basic dynamical models (Generation 0.4).
+
+### Installation
+
+`exptool` is not a registered package, so the best installation is
+
+```
+git clone git@github.com:michael-petersen/exptool.git
+cd exptool
+pip install . --user
+python tests/importtest.py
+```
+
+and then you will be ready to use it! For any updates, you'll need to `git pull` and rebuild using `pip install . --user`.
 
-As of 2020, exptool includes support for analysing stellar halos (Generation 3).
 
 The tests/ directory contains simulation outputs to play with, as well as some usage examples for methods.
 
diff --git a/exptool/__init__.py b/exptool/__init__.py
index b10897f..8ecaedb 100644
--- a/exptool/__init__.py
+++ b/exptool/__init__.py
@@ -1,7 +1,2 @@
 
-
-# this is EXPtool
-
-#http://python-guide-pt-br.readthedocs.io/en/latest/writing/structure/
-
-
+__version__ = "0.4.0"
\ No newline at end of file
diff --git a/exptool/analysis/__init__.py b/exptool/analysis/__init__.py
index 18802fb..09346af 100644
--- a/exptool/analysis/__init__.py
+++ b/exptool/analysis/__init__.py
@@ -1,7 +1,14 @@
 
 
-# this is EXPtool
 
-__all__ = ['commensurability','directdetection','ellipsetools','forcetrace','globalquantities','pattern','trapping']
+__all__ = ['centering',
+           'commensurability',
+           'directdetection',
+           'ellipsetools',
+           'globalquantities',
+           'pattern',
+           'resonancemodel',
+           'rotationcurves',
+           'trapping']
 
 
diff --git a/exptool/analysis/centering.py b/exptool/analysis/centering.py
index 7e44d42..035c671 100644
--- a/exptool/analysis/centering.py
+++ b/exptool/analysis/centering.py
@@ -13,6 +13,7 @@
     
 09 Jul 2021: First introduction
 18 Sep 2024: Class structure overhaul
+23 Oct 2024: Added radius masking to alignment
 """
 import numpy as np
 
@@ -85,7 +86,7 @@ def __init__(self, x, y, z, vx, vy, vz, mass):
         self.vz = vz
         self.mass = mass
 
-    def recentre_pos_and_vel(self, r_max):
+    def recentre_pos_and_vel(self, r_max=np.inf):
         """
         Recenter the position and velocity of the particles based on the center of mass 
         and mass-weighted velocity within r_max.
@@ -95,78 +96,28 @@ def recentre_pos_and_vel(self, r_max):
             Maximum radius to consider for recentering.
 
         Returns:
-        Recentered positions and velocities.
+            None. The operation is performed in place.
         """
         mask = (self.x**2 + self.y**2 + self.z**2) < r_max**2
         mass_tot = np.sum(self.mass[mask])
+
         x_cm = np.sum(self.x[mask]*self.mass[mask])/mass_tot
         y_cm = np.sum(self.y[mask]*self.mass[mask])/mass_tot
         z_cm = np.sum(self.z[mask]*self.mass[mask])/mass_tot
+
         vx_cm = np.sum(self.vx[mask]*self.mass[mask])/mass_tot
         vy_cm = np.sum(self.vy[mask]*self.mass[mask])/mass_tot
         vz_cm = np.sum(self.vz[mask]*self.mass[mask])/mass_tot
-        return self.x - x_cm, self.y - y_cm, self.z - z_cm, self.vx - vx_cm, self.vy - vy_cm, self.vz - vz_cm
-
-    def rotate(self, axis, angle):
-        """
-        Rotate the particles around a given axis by a specified angle.
-
-        Parameters:
-        axis : array
-            Axis of rotation.
-        angle : float
-            Angle of rotation in radians.
-
-        Returns:
-        Rotated positions.
-        """
-        axisx, axisy, axisz = axis
-        dot = self.x * axisx + self.y * axisy + self.z * axisz
-        crossx = axisy * self.z - axisz * self.y
-        crossy = axisz * self.x - axisx * self.z
-        crossz = axisx * self.y - axisy * self.x
-
-        cosa = np.cos(angle)
-        sina = np.sin(angle)
-
-        x_rot = self.x * cosa + crossx * sina + axisx * dot * (1 - cosa)
-        y_rot = self.y * cosa + crossy * sina + axisy * dot * (1 - cosa)
-        z_rot = self.z * cosa + crossz * sina + axisz * dot * (1 - cosa)
-
-        return x_rot, y_rot, z_rot
-
-    def compute_rotation_to_vec(self, vec):
-        """
-        Compute the rotation axis and angle to align the angular momentum of the system 
-        with a given vector.
-
-        Parameters:
-        vec : array
-            Target vector to align with.
-
-        Returns:
-        axis : array
-            Rotation axis.
-        angle : float
-            Angle to rotate.
-        """
-        Lxtot = np.sum(self.y * self.vz * self.mass - self.z * self.vy * self.mass)
-        Lytot = np.sum(self.z * self.vx * self.mass - self.x * self.vz * self.mass)
-        Lztot = np.sum(self.x * self.vy * self.mass - self.y * self.vx * self.mass)
-        
-        L = np.array([Lxtot, Lytot, Lztot])
-        L /= np.linalg.norm(L)
-        vec /= np.linalg.norm(vec)
-
-        axis = np.cross(L, vec)
-        axis /= np.linalg.norm(axis)
-
-        c = np.dot(L, vec)
-        angle = np.arccos(np.clip(c, -1, 1))
 
-        return axis, angle
+        # apply changes in place
+        self.x -= x_cm
+        self.y -= y_cm
+        self.z -= z_cm
+        self.vx -= vx_cm
+        self.vy -= vy_cm
+        self.vz -= vz_cm
 
-    def shrinking_sphere(self, w, rmin=1., stepsize=0.5, tol=0.001, verbose=0):
+    def shrinking_sphere(self, w=None, rmin=1., stepsize=0.5, tol=0.001, verbose=0):
         """
         Apply the shrinking sphere method to recentre the particle system.
 
@@ -183,8 +134,12 @@ def shrinking_sphere(self, w, rmin=1., stepsize=0.5, tol=0.001, verbose=0):
             Level of verbosity for output (default is 0).
 
         Returns:
-        Recentered positions and velocities.
+            None. The operation is performed in place.
         """
+        # If w is None, default to self.mass
+        if w is None:
+            w = self.mass
+
         tshiftx = np.nanmedian(np.nansum(w * self.x) / np.nansum(w))
         tshifty = np.nanmedian(np.nansum(w * self.y) / np.nansum(w))
         tshiftz = np.nanmedian(np.nansum(w * self.z) / np.nansum(w))
@@ -193,14 +148,16 @@ def shrinking_sphere(self, w, rmin=1., stepsize=0.5, tol=0.001, verbose=0):
         self.y -= tshifty
         self.z -= tshiftz
 
-        rval = np.sqrt(self.x**2 + self.y**2 + self.z**2)
-        rmax = np.nanmax(rval)
+        rval  = np.sqrt(self.x**2 + self.y**2 + self.z**2)
+        rmax  = np.nanmax(rval)
         rmax0 = np.nanmax(rval)
 
         if verbose:
             print(f'Initial guess: {tshiftx:.0f}, {tshifty:.0f}, {tshiftz:.0f}')
 
+        nstep = 0
         while rmax > rmin:
+            nstep += 1
             u = np.where(rval < stepsize * rmax)[0]
             if float(u.size) / float(self.x.size) < tol:
                 print(f'Too few particles to continue at radius ratio {stepsize * rmax / rmax0}')
@@ -221,6 +178,9 @@ def shrinking_sphere(self, w, rmin=1., stepsize=0.5, tol=0.001, verbose=0):
             rval = np.sqrt(self.x**2 + self.y**2 + self.z**2)
             rmax *= stepsize
 
+        if nstep == 0:
+            print('No shrinking done. Check that rmin is larger than the minimum particle radius.')
+
         comvx = np.nanmedian(np.nansum(w[u] * self.vx[u]) / np.nansum(w[u]))
         comvy = np.nanmedian(np.nansum(w[u] * self.vy[u]) / np.nansum(w[u]))
         comvz = np.nanmedian(np.nansum(w[u] * self.vz[u]) / np.nansum(w[u]))
@@ -233,9 +193,94 @@ def shrinking_sphere(self, w, rmin=1., stepsize=0.5, tol=0.001, verbose=0):
         self.vy -= comvy
         self.vz -= comvz
 
-        return self.x, self.y, self.z, self.vx, self.vy, self.vz
 
-    def compute_density_profile(self, R, W, rbins=10.**np.linspace(-3.7, 0.3, 100)):
+    def compute_rotation_to_vec(self, vec, r_max=np.inf):
+        """
+        Compute the rotation axis and angle to align the angular momentum of the system 
+        with a given vector.
+
+        Parameters:
+        vec : array
+            Target vector to align with.
+        r_max : float, optional
+            Maximum radius to consider for the alignment (default is infinity).
+
+        Returns:
+        axis : array
+            Rotation axis.
+        angle : float
+            Angle to rotate.
+        """
+        mask = (self.x**2 + self.y**2 + self.z**2) < r_max**2
+
+        Lxtot = np.sum((self.mass * (self.y * self.vz - self.z * self.vy))[mask])
+        Lytot = np.sum((self.mass * (self.z * self.vx - self.x * self.vz))[mask])
+        Lztot = np.sum((self.mass * (self.x * self.vy - self.y * self.vx))[mask])
+        
+        L = np.array([Lxtot, Lytot, Lztot])
+        L /= np.linalg.norm(L)
+
+        vec /= np.linalg.norm(vec)
+
+        axis = np.cross(L, vec)
+        axis /= np.linalg.norm(axis)
+
+        c = np.dot(L, vec)
+        angle = np.arccos(np.clip(c, -1, 1))
+
+        return axis, angle
+    
+    def rotate(self, vec=[0.,0.,1.], r_max=np.inf):
+        """
+        Rotate the particles around a given axis by a specified angle.
+
+        Parameters:
+        vec : array, optional
+            Target vector to align with (default is [0, 0, 1]; i.e. angular momentum aligned with z-axis).
+        r_max : float, optional
+            Maximum radius to consider for the alignment (default is infinity).
+
+        Returns:
+            None. The operation is performed in place.
+        """
+        axis, angle = self.compute_rotation_to_vec(vec,r_max=r_max)
+        
+        axisx, axisy, axisz = axis
+
+        # Rotation for position vector (x, y, z)
+        dot_pos = self.x * axisx + self.y * axisy + self.z * axisz  # Dot product for position
+        crossx_pos = axisy * self.z - axisz * self.y  # Cross product for position
+        crossy_pos = axisz * self.x - axisx * self.z
+        crossz_pos = axisx * self.y - axisy * self.x
+
+        cosa = np.cos(angle)
+        sina = np.sin(angle)
+
+        x_rot = self.x * cosa + crossx_pos * sina + axisx * dot_pos * (1 - cosa)
+        y_rot = self.y * cosa + crossy_pos * sina + axisy * dot_pos * (1 - cosa)
+        z_rot = self.z * cosa + crossz_pos * sina + axisz * dot_pos * (1 - cosa)
+
+        # Rotation for velocity vector (vx, vy, vz)
+        dot_vel = self.vx * axisx + self.vy * axisy + self.vz * axisz  # Dot product for velocity
+        crossx_vel = axisy * self.vz - axisz * self.vy  # Cross product for velocity
+        crossy_vel = axisz * self.vx - axisx * self.vz
+        crossz_vel = axisx * self.vy - axisy * self.vx
+
+        vx_rot = self.vx * cosa + crossx_vel * sina + axisx * dot_vel * (1 - cosa)
+        vy_rot = self.vy * cosa + crossy_vel * sina + axisy * dot_vel * (1 - cosa)
+        vz_rot = self.vz * cosa + crossz_vel * sina + axisz * dot_vel * (1 - cosa)
+
+        # perform operation in place
+        self.x = x_rot
+        self.y = y_rot
+        self.z = z_rot
+        self.vx = vx_rot
+        self.vy = vy_rot
+        self.vz = vz_rot
+
+        #return x_rot, y_rot, z_rot, vx_rot, vy_rot, vz_rot
+
+    def compute_density_profile(self, R, W, rbins=10.**np.linspace(-3.7, 0.3, 100),astronomicalG = 0.0000043009125):
         """
         Compute the density profile, enclosed mass, and potential for the particle system.
 
@@ -246,6 +291,8 @@ def compute_density_profile(self, R, W, rbins=10.**np.linspace(-3.7, 0.3, 100)):
             Weights of the particles.
         rbins : array, optional
             Radial bins to compute the profile (default is logarithmic bins from 10^-3.7 to 10^0.3).
+        astronomicalG : float, optional
+            Gravitational constant (default is 0.0000043009125 km/s/kpc/Msun).
 
         Returns:
         dens : array
@@ -259,7 +306,6 @@ def compute_density_profile(self, R, W, rbins=10.**np.linspace(-3.7, 0.3, 100)):
         menc = np.zeros(rbins.size)
         potp = np.zeros(rbins.size)
 
-        astronomicalG = 0.0000043009125
         rbinstmp = np.concatenate([rbins, [2. * rbins[-1] - rbins[-2]]])
 
         for indx, val in enumerate(rbinstmp[:-1]):
diff --git a/exptool/analysis/resonancemodel.py b/exptool/analysis/resonancemodel.py
index a66b94e..5456589 100644
--- a/exptool/analysis/resonancemodel.py
+++ b/exptool/analysis/resonancemodel.py
@@ -59,15 +59,15 @@ def denom(r,E,J,model):
 
 def make_orbit(orbit,E,K,model):
     """make an orbit"""
-      orbit.ee = E
-      #
-      # this should work, the boundaries are in radius...
-      orbit.r_circ = brentq(Ecirc,np.min(model.rcurve),np.max(model.rcurve),args=(orbit.ee,model))
-      orbit.kappa = K
-      orbit.jj = find_j(orbit.r_circ,orbit.kappa,model)
-      orbit.r_apo = brentq(denom,orbit.r_circ,np.max(model.rcurve),args=(orbit.ee,orbit.jj,model))
-      orbit.r_peri = brentq(denom,np.min(model.rcurve),orbit.r_circ,args=(orbit.ee,orbit.jj,model))
-      return orbit
+    orbit.ee = E
+    #
+    # this should work, the boundaries are in radius...
+    orbit.r_circ = brentq(Ecirc,np.min(model.rcurve),np.max(model.rcurve),args=(orbit.ee,model))
+    orbit.kappa = K
+    orbit.jj = find_j(orbit.r_circ,orbit.kappa,model)
+    orbit.r_apo = brentq(denom,orbit.r_circ,np.max(model.rcurve),args=(orbit.ee,orbit.jj,model))
+    orbit.r_peri = brentq(denom,np.min(model.rcurve),orbit.r_circ,args=(orbit.ee,orbit.jj,model))
+    return orbit
 
 
 
diff --git a/exptool/analysis/rotationcurves.py b/exptool/analysis/rotationcurves.py
index 3b6a1f1..9767d83 100644
--- a/exptool/analysis/rotationcurves.py
+++ b/exptool/analysis/rotationcurves.py
@@ -2,6 +2,9 @@
 # tools to fit rotation curves. NEEDS MUCH CONSTRUCTION
 # from NewHorizon bars project
 
+import numpy as np
+astronomicalG = 0.0000043009125 # gravitational constant, (km/s)^2 * kpc / Msun
+
 
 def kuzmin_rotation(R,c,M,G=astronomicalG):
     """see BT08
diff --git a/exptool/basis/eof.py b/exptool/basis/eof.py
index 8b3e2ab..b686b48 100644
--- a/exptool/basis/eof.py
+++ b/exptool/basis/eof.py
@@ -682,12 +682,12 @@ def map_basis(eof_file):
     rmin,rmax,numx,numy,mmax,norder,ascale,hscale,cmap,dens = eof_params(eof_file)
 
     if (dens):
-        mC = np.memmap(eof_file, dtype=np.float64, offset=76, shape=(mmax+1,norder,4,numx+1,numy+1))
-        mS = np.memmap(eof_file, dtype=np.float64, offset=76+(8*4*(mmax+1)*norder*(numx+1)*(numy+1)), shape=(mmax,norder,4,numx+1,numy+1))
+        mC = np.memmap(eof_file, dtype=float64, offset=76, shape=(mmax+1,norder,4,numx+1,numy+1))
+        mS = np.memmap(eof_file, dtype=float64, offset=76+(8*4*(mmax+1)*norder*(numx+1)*(numy+1)), shape=(mmax,norder,4,numx+1,numy+1))
 
     else:
-        mC = np.memmap(eof_file, dtype=np.float64, offset=76, shape=(mmax+1,norder,3,numx+1,numy+1))
-        mS = np.memmap(eof_file, dtype=np.float64, offset=76+(8*3*(mmax+1)*norder*(numx+1)*(numy+1)), shape=(mmax,norder,3,numx+1,numy+1))
+        mC = np.memmap(eof_file, dtype=float64, offset=76, shape=(mmax+1,norder,3,numx+1,numy+1))
+        mS = np.memmap(eof_file, dtype=float64, offset=76+(8*3*(mmax+1)*norder*(numx+1)*(numy+1)), shape=(mmax,norder,3,numx+1,numy+1))
 
     # note that mS and potS are different sized owing to m orders: need to homogenize for full usefulness?
 
@@ -2258,7 +2258,7 @@ def read_binary_eof_coefficients(coeffile):
     # return to beginning
     f.seek(0)
 
-    [time0] = np.fromfile(f, dtype=np.float,count=1)
+    [time0] = np.fromfile(f, dtype=float,count=1)
     [mmax,nmax] = np.fromfile(f, dtype=np.uint32,count=2)
 
     # hard-coded to match specifications.
@@ -2274,17 +2274,17 @@ def read_binary_eof_coefficients(coeffile):
 
     for tt in range(0,n_outputs):
 
-        [time0] = np.fromfile(f, dtype=np.float,count=1)
+        [time0] = np.fromfile(f, dtype=float,count=1)
         [dummym,dummyn] = np.fromfile(f, dtype=np.uint32,count=2)
 
         times[tt] = time0
 
         for mm in range(0,mmax+1):
 
-            coef_array[tt,0,mm,:] = np.fromfile(f, dtype=np.float,count=nmax)
+            coef_array[tt,0,mm,:] = np.fromfile(f, dtype=float,count=nmax)
 
             if mm > 0:
-                coef_array[tt,1,mm,:] = np.fromfile(f, dtype=np.float,count=nmax)
+                coef_array[tt,1,mm,:] = np.fromfile(f, dtype=float,count=nmax)
 
 
     return times,coef_array
@@ -2320,7 +2320,7 @@ def read_binary_eof_coefficients_dict(coeffile):
     # return to beginning
     f.seek(0)
 
-    [time0] = np.fromfile(f, dtype=np.float,count=1)
+    [time0] = np.fromfile(f, dtype=float,count=1)
     [mmax,nmax] = np.fromfile(f, dtype=np.uint32,count=2)
 
     # hard-coded to match specifications.
@@ -2336,7 +2336,7 @@ def read_binary_eof_coefficients_dict(coeffile):
 
         EOF_Obj = EOF_Object()
 
-        [EOF_Obj.time] = np.fromfile(f, dtype=np.float,count=1)
+        [EOF_Obj.time] = np.fromfile(f, dtype=float,count=1)
         [EOF_Obj.mmax,EOF_Obj.nmax] = np.fromfile(f, dtype=np.uint32,count=2)
 
         # fill in dummy values
@@ -2350,10 +2350,10 @@ def read_binary_eof_coefficients_dict(coeffile):
 
         for mm in range(0,mmax+1):
 
-            EOF_Obj.cos[mm,:] = np.fromfile(f, dtype=np.float,count=nmax)
+            EOF_Obj.cos[mm,:] = np.fromfile(f, dtype=float,count=nmax)
 
             if mm > 0:
-                EOF_Obj.sin[mm,:] = np.fromfile(f, dtype=np.float,count=nmax)
+                EOF_Obj.sin[mm,:] = np.fromfile(f, dtype=float,count=nmax)
 
         EOF_Dict[EOF_Obj.time] = EOF_Obj
 
diff --git a/exptool/observables/__init__.py b/exptool/observables/__init__.py
index 2c28ec4..67451ac 100644
--- a/exptool/observables/__init__.py
+++ b/exptool/observables/__init__.py
@@ -1,4 +1,13 @@
 
 
-# this is EXPtool
+__all__ = ['bulge',
+           'fourier',
+           'gaiatransform',
+           'rotationcurve',
+           'scaleheightcheck',
+           'transform',
+           'velocity',
+           'visualize']
+
+
 
diff --git a/exptool/observables/fourier.py b/exptool/observables/fourier.py
index 293458f..150df73 100644
--- a/exptool/observables/fourier.py
+++ b/exptool/observables/fourier.py
@@ -100,7 +100,7 @@ def within_radius(self, radius):
         r = np.sqrt(self.xpos**2 + self.ypos**2 + self.zpos**2)
         return r < radius
     
-    def compute_fourier(self, harmonic=2):
+    def compute_fourier(self, mask=None, harmonic=2, normalize=True):
         """
         Compute Fourier moments for the particle system.
 
@@ -111,14 +111,36 @@ def compute_fourier(self, harmonic=2):
         - mod: The modulus of the Fourier moment.
         - angle: The phase angle of the Fourier moment.
         """
-        phi = np.arctan2(self.y, self.x)
-        A = np.nansum(self.mass * np.cos(harmonic * phi)) / np.nansum(self.mass)
-        B = np.nansum(self.mass * np.sin(harmonic * phi)) / np.nansum(self.mass)
-        mod = np.sqrt(A * A + B * B)
-        angle = np.arctan2(B, A) / 2.
+
+        norm = 1.0
+
+        if mask is None:
+            phi = np.arctan2(self.y, self.x)
+
+            if normalize:
+                norm = np.nansum(self.mass)
+
+            A = np.nansum(self.mass * np.cos(harmonic * phi)) / norm
+            B = np.nansum(self.mass * np.sin(harmonic * phi)) / norm
+
+            mod = np.sqrt(A * A + B * B)
+            angle = np.arctan2(B, A) / 2.
+
+        else:
+            phi = np.arctan2(self.y[mask], self.x[mask])
+
+            if normalize:
+                norm = np.nansum(self.mass[mask])
+
+            A = np.nansum(self.mass[mask] * np.cos(harmonic * phi)) / norm
+            B = np.nansum(self.mass[mask] * np.sin(harmonic * phi)) / norm
+
+            mod = np.sqrt(A * A + B * B)
+            angle = np.arctan2(B, A) / 2.
+
         return mod, angle
     
-    def compute_fourier_vel(self, harmonic=2):
+    def compute_fourier_vel(self, mask=None, harmonic=2):
         """
         Compute velocity-weighted Fourier moments.
 
@@ -128,15 +150,22 @@ def compute_fourier_vel(self, harmonic=2):
         Returns:
         - mod: The modulus of the velocity-weighted Fourier moment.
         """
-        phi = np.arctan2(self.y, self.x)
-        # Compute the tangential velocity
-        vel = (self.x * self.vy - self.y * self.vx) / np.sqrt(self.x**2 + self.y**2)
-        A = np.nansum(self.mass * vel * np.cos(harmonic * phi)) / np.nansum(self.mass)
-        B = np.nansum(self.mass * vel * np.sin(harmonic * phi)) / np.nansum(self.mass)
-        mod = np.sqrt(A * A + B * B)
+        if mask is None:
+            phi = np.arctan2(self.y, self.x)
+            # Compute the tangential velocity
+            vel = (self.x * self.vy - self.y * self.vx) / np.sqrt(self.x**2 + self.y**2)
+            A = np.nansum(self.mass * vel * np.cos(harmonic * phi)) / np.nansum(self.mass)
+            B = np.nansum(self.mass * vel * np.sin(harmonic * phi)) / np.nansum(self.mass)
+            mod = np.sqrt(A * A + B * B)
+        else:
+            phi = np.arctan2(self.y[mask], self.x[mask])
+            vel = (self.x[mask] * self.vy[mask] - self.y[mask] * self.vx[mask]) / np.sqrt(self.x[mask]**2 + self.y[mask]**2)
+            A = np.nansum(self.mass[mask] * vel * np.cos(harmonic * phi)) / np.nansum(self.mass[mask])
+            B = np.nansum(self.mass[mask] * vel * np.sin(harmonic * phi)) / np.nansum(self.mass[mask])
+            mod = np.sqrt(A * A + B * B)
         return mod
     
-    def fourier_tabulate(self, bins=75):
+    def fourier_tabulate(self, bins=75,mmax=2, r_max=None,normalize=True):
         """
         Compute the Fourier moments and velocity-weighted Fourier moments in radial bins.
 
@@ -150,19 +179,25 @@ def fourier_tabulate(self, bins=75):
         - fvpower: The velocity-weighted Fourier moments for harmonics 0 to 4.
         """
         rval = np.sqrt(self.x**2 + self.y**2)
-        rtest = np.linspace(0., np.nanpercentile(rval, 75), bins)
+
+        if r_max is None:
+            rtest = np.linspace(0., np.nanpercentile(rval, 75), bins)
+        else:
+            rtest = np.linspace(0., r_max, bins)
+
+
         dr = rtest[1] - rtest[0]
         rbin = (np.floor(rval / dr)).astype(int)
         
-        fpower = np.zeros((5, rtest.size))
-        fangle = np.zeros((5, rtest.size))
-        fvpower = np.zeros((5, rtest.size))
+        fpower = np.zeros((mmax+1, rtest.size))
+        fangle = np.zeros((mmax+1, rtest.size))
+        fvpower = np.zeros((mmax+1, rtest.size))
         
         for ir, rv in enumerate(rtest):
             w = np.where(rbin == ir)[0]
-            for h in range(5):
-                fpower[h, ir], fangle[h, ir] = self.compute_fourier(harmonic=h)
-                fvpower[h, ir] = self.compute_fourier_vel(harmonic=h)
+            for h in range(mmax+1):
+                fpower[h, ir], fangle[h, ir] = self.compute_fourier(harmonic=h,mask=w,normalize=normalize)
+                fvpower[h, ir] = self.compute_fourier_vel(harmonic=h,mask=w)
         
         return rtest, fpower, fangle, fvpower
 
diff --git a/exptool/observables/visualize.py b/exptool/observables/visualize.py
index 27f2ec6..090470f 100644
--- a/exptool/observables/visualize.py
+++ b/exptool/observables/visualize.py
@@ -302,13 +302,13 @@ def show_dump(infile,comp,nout=-1,type='pos',transform=True,\
     if np.array(comp).size == 1:
 
         if nout < 0:
-            PSPDump = particle.Input(infile,comp=comp)
+            PSPDump = particle.Input(infile,comp=comp,legacy=True)
         else:
-            PSPDump = particle.Input(infile,comp=comp,nout=nout)
+            PSPDump = particle.Input(infile,comp=comp,nout=nout,legacy=True)
 
     else:
         # allow for multiple components to be mixed together
-        PartArray = [particle.Input(infile,comp=cc) for cc in comp]
+        PartArray = [particle.Input(infile,comp=cc,legacy=True) for cc in comp]
 
         PSPDump = particle.mix_particles(PartArray)
 
@@ -484,8 +484,8 @@ def compare_dumps(infile1,infile2,comp,type='pos',transform=True,\
 
     # read in files
 
-    PSPDump1 = particle.Input(infile1,comp=comp)
-    PSPDump2 = particle.Input(infile2,comp=comp)
+    PSPDump1 = particle.Input(infile1,comp=comp,legacy=True)
+    PSPDump2 = particle.Input(infile2,comp=comp,legacy=True)
 
     if transform:
         PSPDump1 = pattern.BarTransform(PSPDump1)
diff --git a/exptool/observables/movie.py b/exptool/tests/movie.py
similarity index 96%
rename from exptool/observables/movie.py
rename to exptool/tests/movie.py
index ab52d7f..d16d094 100644
--- a/exptool/observables/movie.py
+++ b/exptool/tests/movie.py
@@ -22,10 +22,10 @@
 
 
 # exptool imports
-from .io import particle
-from .utils import kde_3d
-from .observables import transform
-from .analysis import pattern
+from ..io import particle
+from ..utils import kde_3d
+from ..observables import transform
+from ..analysis import pattern
 
 
 # location of the files to become a movie: satellite edition
diff --git a/exptool/utils/utils.py b/exptool/utils/utils.py
index 0a15082..65de9e6 100644
--- a/exptool/utils/utils.py
+++ b/exptool/utils/utils.py
@@ -442,8 +442,8 @@ def savitzky_golay(y, window_size, order, deriv=0, rate=1):
 
     
     try:
-        window_size = np.abs(np.int(window_size))
-        order = np.abs(np.int(order))
+        window_size = np.abs(int(window_size))
+        order = np.abs(int(order))
     except:# ValueError, msg:
         raise ValueError("window_size and order have to be of type int")
     
diff --git a/test_requirements.txt b/test_requirements.txt
index 7b500b3..75afa50 100644
--- a/test_requirements.txt
+++ b/test_requirements.txt
@@ -1,3 +1,5 @@
 numpy==1.26.4
 scipy==1.11.4
+scikit-image==0.22.0
+pyyaml==6.0.2
 -e .
\ No newline at end of file
diff --git a/tests/importtest.py b/tests/importtest.py
new file mode 100644
index 0000000..02ce3c9
--- /dev/null
+++ b/tests/importtest.py
@@ -0,0 +1,4 @@
+
+
+from exptool.analysis import *
+from exptool.observables import *
\ No newline at end of file