diff --git a/examples/2D_ibm_ellipse/case.py b/examples/2D_ibm_ellipse/case.py
index 1c5d8052f7..ceb9fda4b6 100644
--- a/examples/2D_ibm_ellipse/case.py
+++ b/examples/2D_ibm_ellipse/case.py
@@ -81,7 +81,7 @@
             "patch_icpp(1)%pres": 1.0e00,
             "patch_icpp(1)%alpha_rho(1)": 1.0e00,
             "patch_icpp(1)%alpha(1)": 1.0e00,
-            # Patch: Cylinder Immersed Boundary
+            # Patch: Ellipse Immersed Boundary
             "patch_ib(1)%geometry": 6,
             "patch_ib(1)%x_centroid": 1.5e-03,
             "patch_ib(1)%y_centroid": 3.0e-03,
diff --git a/examples/2D_kelvin_helmholtz/case.py b/examples/2D_kelvin_helmholtz/case.py
index 3a12034696..e6873f8145 100644
--- a/examples/2D_kelvin_helmholtz/case.py
+++ b/examples/2D_kelvin_helmholtz/case.py
@@ -2,7 +2,6 @@
 import json
 import math
 
-
 eps = 1e-6
 time_end = 1.0
 time_save = time_end / 100.0
diff --git a/examples/2D_richtmyer_meshkov/case.py b/examples/2D_richtmyer_meshkov/case.py
index efb30f17f4..fed6529c55 100644
--- a/examples/2D_richtmyer_meshkov/case.py
+++ b/examples/2D_richtmyer_meshkov/case.py
@@ -2,7 +2,6 @@
 import json
 import math
 
-
 mu = 1.0e-4
 lambd = 1.0
 time_end = 15.0
diff --git a/examples/2D_viscous_shock_tube/case.py b/examples/2D_viscous_shock_tube/case.py
index e61f34cb97..7c760e696b 100644
--- a/examples/2D_viscous_shock_tube/case.py
+++ b/examples/2D_viscous_shock_tube/case.py
@@ -2,7 +2,6 @@
 import json
 import math
 
-
 mu = 1.0e-3
 time_end = 1.0
 time_save = time_end / 100.0
diff --git a/examples/nD_perfect_reactor/analyze.py b/examples/nD_perfect_reactor/analyze.py
index ab0c6c0a94..b2bd4e1fa8 100644
--- a/examples/nD_perfect_reactor/analyze.py
+++ b/examples/nD_perfect_reactor/analyze.py
@@ -6,7 +6,6 @@
 import mfc.viz
 from case import dt, Tend, SAVE_COUNT, sol
 
-
 case = mfc.viz.Case(".", dt)
 
 sns.set_theme(style=mfc.viz.generate_cpg_style())
diff --git a/examples/nD_perfect_reactor/export.py b/examples/nD_perfect_reactor/export.py
index f60f55e3fe..112622a32f 100644
--- a/examples/nD_perfect_reactor/export.py
+++ b/examples/nD_perfect_reactor/export.py
@@ -5,7 +5,6 @@
 import mfc.viz
 from case import dt, NS, Tend, SAVE_COUNT, sol
 
-
 case = mfc.viz.Case(".", dt)
 
 for name in tqdm(sol.species_names, desc="Loading Variables"):
diff --git a/src/common/include/case.fpp b/src/common/include/case.fpp
index 84d1d2cc14..aa0e0637b9 100644
--- a/src/common/include/case.fpp
+++ b/src/common/include/case.fpp
@@ -6,3 +6,8 @@
 #:def analytical()
 
 #:enddef
+
+! For moving immersed boundaries in simulation
+#:def mib_analytical()
+
+#:enddef
diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp
index c7a3df21cf..da2b29956c 100644
--- a/src/common/m_compute_levelset.fpp
+++ b/src/common/m_compute_levelset.fpp
@@ -94,6 +94,7 @@ contains
             do j = 0, n
                 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
                 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate
+                xy_local = xy_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
 
                 if (xy_local(2) >= 0._wp) then
                     ! finds the location on the airfoil grid with the minimum distance (closest)
@@ -189,6 +190,7 @@ contains
 
                     xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
                     xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
+                    xyz_local = xyz_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
 
                     if (xyz_local(2) >= center(2)) then
                         do k = 1, Np
diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
index f3ef91e831..5c36a7ae11 100644
--- a/src/common/m_derived_types.fpp
+++ b/src/common/m_derived_types.fpp
@@ -296,6 +296,7 @@ module m_derived_types
         !! is specified through its x-, y- and z-coordinates, respectively.
         real(wp) :: step_x_centroid, step_y_centroid, step_z_centroid !<
         !! Centroid locations of intermediate steps in the time_stepper module
+        real(wp), dimension(1:3) :: centroid_offset ! offset of center of mass from computed cell center for odd-shaped IBs
 
         real(wp), dimension(1:3) :: angles
         real(wp), dimension(1:3) :: step_angles
diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp
index 1bb00c44a3..3db86ceab5 100644
--- a/src/common/m_ib_patches.fpp
+++ b/src/common/m_ib_patches.fpp
@@ -279,6 +279,7 @@ contains
             do i = 0, m
                 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
                 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates
+                xy_local = xy_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
 
                 if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then
                     xa = xy_local(1)/ca_in
@@ -433,6 +434,7 @@ contains
                 do i = 0, m
                     xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
                     xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
+                    xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
 
                     if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then
 
diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp
index ea679ba6f5..c95e1a2946 100644
--- a/src/pre_process/m_global_parameters.fpp
+++ b/src/pre_process/m_global_parameters.fpp
@@ -572,6 +572,7 @@ contains
             patch_ib(i)%angular_vel(:) = 0._wp
             patch_ib(i)%mass = dflt_real
             patch_ib(i)%moment = dflt_real
+            patch_ib(i)%centroid_offset(:) = 0._wp
 
             ! sets values of a rotation matrix which can be used when calculating rotations
             patch_ib(i)%rotation_matrix = 0._wp
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
index 0b8ec1ab81..bccb194d42 100644
--- a/src/simulation/m_ibm.fpp
+++ b/src/simulation/m_ibm.fpp
@@ -95,6 +95,7 @@ contains
         integer :: i, j, k
         integer :: max_num_gps, max_num_inner_gps
 
+        ! do all set up for moving immersed boundaries
         moving_immersed_boundary_flag = .false.
         do i = 1, num_ibs
             if (patch_ib(i)%moving_ibm /= 0) then
@@ -102,6 +103,7 @@ contains
                 moving_immersed_boundary_flag = .true.
             end if
             call s_update_ib_rotation_matrix(i)
+            call s_compute_centroid_offset(i)
         end do
         $:GPU_ENTER_DATA(copyin='[patch_ib]')
 
@@ -1012,24 +1014,24 @@ contains
     ! compute the surface integrals of the IB via a volume integraion method described in
     ! "A coupled IBM/Euler-Lagrange framework for simulating shock-induced particle size segregation"
     ! by Archana Sridhar and Jesse Capecelatro
-    subroutine s_compute_ib_forces(q_prim_vf, dynamic_viscosity)
+    subroutine s_compute_ib_forces(q_prim_vf, fluid_pp)
 
         ! real(wp), dimension(idwbuff(1)%beg:idwbuff(1)%end, &
         !             idwbuff(2)%beg:idwbuff(2)%end, &
         !             idwbuff(3)%beg:idwbuff(3)%end), intent(in) :: pressure
         type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf
-        real(wp), intent(in) :: dynamic_viscosity
+        type(physical_parameters), dimension(1:num_fluids), intent(in) :: fluid_pp
 
-        integer :: gp_id, i, j, k, l, q, ib_idx
+        integer :: gp_id, i, j, k, l, q, ib_idx, fluid_idx
         real(wp), dimension(num_ibs, 3) :: forces, torques
         real(wp), dimension(1:3, 1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2 ! viscous stress tensor with temp vectors to hold divergence calculations
         real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution, vel
-        real(wp) :: cell_volume, dx, dy, dz
+        real(wp) :: cell_volume, dx, dy, dz, dynamic_viscosity
 
         forces = 0._wp
         torques = 0._wp
 
-        $:GPU_PARALLEL_LOOP(private='[ib_idx,radial_vector,local_force_contribution,cell_volume,local_torque_contribution, viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2, dx, dy, dz]', copy='[forces,torques]', copyin='[ib_markers,patch_ib,dynamic_viscosity]', collapse=3)
+        $:GPU_PARALLEL_LOOP(private='[ib_idx,fluid_idx, radial_vector,local_force_contribution,cell_volume,local_torque_contribution, dynamic_viscosity, viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2, dx, dy, dz]', copy='[forces,torques]', copyin='[ib_markers,patch_ib,fluid_pp]', collapse=3)
         do i = 0, m
             do j = 0, n
                 do k = 0, p
@@ -1044,26 +1046,33 @@ contains
                         dx = x_cc(i + 1) - x_cc(i)
                         dy = y_cc(j + 1) - y_cc(j)
 
-                        ! Get the pressure contribution to force via a finite difference to compute the 2D components of the gradient of the pressure and cell volume
-                        local_force_contribution(1) = -1._wp*(q_prim_vf(E_idx)%sf(i + 1, j, k) - q_prim_vf(E_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient
-                        local_force_contribution(2) = -1._wp*(q_prim_vf(E_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx)%sf(i, j - 1, k))/(2._wp*dy)
-                        cell_volume = abs(dx*dy)
-                        ! add the 3D component of the pressure gradient, if we are working in 3 dimensions
-                        if (num_dims == 3) then
-                            dz = z_cc(k + 1) - z_cc(k)
-                            local_force_contribution(3) = -1._wp*(q_prim_vf(E_idx)%sf(i, j, k + 1) - q_prim_vf(E_idx)%sf(i, j, k - 1))/(2._wp*dz)
-                            cell_volume = abs(cell_volume*dz)
-                        else
-                            local_force_contribution(3) = 0._wp
-                        end if
+                        local_force_contribution(:) = 0._wp
+                        do fluid_idx = 0, num_fluids - 1
+                            ! Get the pressure contribution to force via a finite difference to compute the 2D components of the gradient of the pressure and cell volume
+                            local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i + 1, j, k) - q_prim_vf(E_idx + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient
+                            local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
+                            cell_volume = abs(dx*dy)
+                            ! add the 3D component of the pressure gradient, if we are working in 3 dimensions
+                            if (num_dims == 3) then
+                                dz = z_cc(k + 1) - z_cc(k)
+                                local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j, k + 1) - q_prim_vf(E_idx + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz)
+                                cell_volume = abs(cell_volume*dz)
+                            end if
+                        end do
 
                         ! Update the force values atomically to prevent race conditions
                         call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)
 
                         ! get the viscous stress and add its contribution if that is considered
                         ! TODO :: This is really bad code
-                        ! if (.false.) then
                         if (viscous) then
+                            ! compute the volume-weighted local dynamic viscosity
+                            dynamic_viscosity = 0._wp
+                            do fluid_idx = 1, num_fluids
+                                ! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid
+                                if (fluid_pp(fluid_idx)%Re(1) /= 0._wp) dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + advxb - 1)%sf(i, j, k)/fluid_pp(fluid_idx)%Re(1))
+                            end do
+
                             ! get the linear force component first
                             call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, j, k)
                             call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, j, k)
@@ -1150,6 +1159,64 @@ contains
 
     end subroutine s_finalize_ibm_module
 
+    !> Computes the center of mass for IB patch types where we are unable to determine their center of mass analytically.
+    !> These patches include things like NACA airfoils and STL models
+    subroutine s_compute_centroid_offset(ib_marker)
+
+        integer, intent(in) :: ib_marker
+
+        integer :: i, j, k, num_cells, num_cells_local
+        real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local
+
+        ! Offset only needs to be computes for specific geometries
+        if (patch_ib(ib_marker)%geometry == 4 .or. &
+            patch_ib(ib_marker)%geometry == 5 .or. &
+            patch_ib(ib_marker)%geometry == 11 .or. &
+            patch_ib(ib_marker)%geometry == 12) then
+
+            center_of_mass_local = [0._wp, 0._wp, 0._wp]
+            num_cells_local = 0
+
+            ! get the summed mass distribution and number of cells to divide by
+            do i = 0, m
+                do j = 0, n
+                    do k = 0, p
+                        if (ib_markers%sf(i, j, k) == ib_marker) then
+                            num_cells_local = num_cells_local + 1
+                            center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp]
+                            if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k)
+                        end if
+                    end do
+                end do
+            end do
+
+            ! reduce the mass contribution over all MPI ranks and compute COM
+            call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
+            if (num_cells /= 0) then
+                call s_mpi_allreduce_sum(center_of_mass_local(1), center_of_mass(1))
+                call s_mpi_allreduce_sum(center_of_mass_local(2), center_of_mass(2))
+                call s_mpi_allreduce_sum(center_of_mass_local(3), center_of_mass(3))
+                center_of_mass = center_of_mass/real(num_cells, wp)
+            else
+                patch_ib(ib_marker)%centroid_offset = [0._wp, 0._wp, 0._wp]
+                return
+            end if
+
+            ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the current centroid
+            patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid] &
+                                                  - center_of_mass
+            patch_ib(ib_marker)%x_centroid = center_of_mass(1)
+            patch_ib(ib_marker)%y_centroid = center_of_mass(2)
+            patch_ib(ib_marker)%z_centroid = center_of_mass(3)
+
+            ! rotate the centroid offset back into the local coords of the IB
+            patch_ib(ib_marker)%centroid_offset = matmul(patch_ib(ib_marker)%rotation_matrix_inverse, patch_ib(ib_marker)%centroid_offset)
+        else
+            patch_ib(ib_marker)%centroid_offset(:) = [0._wp, 0._wp, 0._wp]
+        end if
+
+    end subroutine s_compute_centroid_offset
+
     subroutine s_compute_moment_of_inertia(ib_marker, axis)
 
         real(wp), dimension(3), intent(in) :: axis !< the axis about which we compute the moment. Only required in 3D.
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index b196b3afa3..3502326e85 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -3,6 +3,7 @@
 !! @brief Contains module m_time_steppers
 
 #:include 'macros.fpp'
+#:include 'case.fpp'
 
 !> @brief The following module features a variety of time-stepping schemes.
 !!              Currently, it includes the following Runge-Kutta (RK) algorithms:
@@ -611,44 +612,10 @@ contains
             if (ib) then
                 ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms
                 if (moving_immersed_boundary_flag) then
-                    do i = 1, num_ibs
-                        if (s == 1) then
-                            patch_ib(i)%step_vel = patch_ib(i)%vel
-                            patch_ib(i)%step_angular_vel = patch_ib(i)%angular_vel
-                            patch_ib(i)%step_angles = patch_ib(i)%angles
-                            patch_ib(i)%step_x_centroid = patch_ib(i)%x_centroid
-                            patch_ib(i)%step_y_centroid = patch_ib(i)%y_centroid
-                            patch_ib(i)%step_z_centroid = patch_ib(i)%z_centroid
-                        end if
-
-                        if (patch_ib(i)%moving_ibm > 0) then
-                            patch_ib(i)%vel = (rk_coef(s, 1)*patch_ib(i)%step_vel + rk_coef(s, 2)*patch_ib(i)%vel)/rk_coef(s, 4)
-                            patch_ib(i)%angular_vel = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel + rk_coef(s, 2)*patch_ib(i)%angular_vel)/rk_coef(s, 4)
-
-                            if (patch_ib(i)%moving_ibm == 2) then ! if we are using two-way coupling, apply force and torque
-                                ! compute the force and torque on the IB from the fluid
-                                call s_compute_ib_forces(q_prim_vf, 1._wp/fluid_pp(1)%Re(1))
-
-                                ! update the velocity from the force value
-                                patch_ib(i)%vel = patch_ib(i)%vel + rk_coef(s, 3)*dt*(patch_ib(i)%force/patch_ib(i)%mass)/rk_coef(s, 4)
-
-                                ! update the angular velocity with the torque value
-                                patch_ib(i)%angular_vel = (patch_ib(i)%angular_vel*patch_ib(i)%moment) + (rk_coef(s, 3)*dt*patch_ib(i)%torque/rk_coef(s, 4)) ! add the torque to the angular momentum
-                                call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel) ! update the moment of inertia to be based on the direction of the angular momentum
-                                patch_ib(i)%angular_vel = patch_ib(i)%angular_vel/patch_ib(i)%moment ! convert back to angular velocity with the new moment of inertia
-                            end if
-
-                            ! Update the angle of the IB
-                            patch_ib(i)%angles = (rk_coef(s, 1)*patch_ib(i)%step_angles + rk_coef(s, 2)*patch_ib(i)%angles + rk_coef(s, 3)*patch_ib(i)%angular_vel*dt)/rk_coef(s, 4)
-
-                            ! Update the position of the IB
-                            patch_ib(i)%x_centroid = (rk_coef(s, 1)*patch_ib(i)%step_x_centroid + rk_coef(s, 2)*patch_ib(i)%x_centroid + rk_coef(s, 3)*patch_ib(i)%vel(1)*dt)/rk_coef(s, 4)
-                            patch_ib(i)%y_centroid = (rk_coef(s, 1)*patch_ib(i)%step_y_centroid + rk_coef(s, 2)*patch_ib(i)%y_centroid + rk_coef(s, 3)*patch_ib(i)%vel(2)*dt)/rk_coef(s, 4)
-                            patch_ib(i)%z_centroid = (rk_coef(s, 1)*patch_ib(i)%step_z_centroid + rk_coef(s, 2)*patch_ib(i)%z_centroid + rk_coef(s, 3)*patch_ib(i)%vel(3)*dt)/rk_coef(s, 4)
-                        end if
-                    end do
-                    call s_update_mib(num_ibs, levelset, levelset_norm)
+                    call s_propagate_immersed_boundaries(s)
                 end if
+
+                ! update the ghost fluid properties point values based on IB state
                 if (qbmm .and. .not. polytropic) then
                     call s_ibm_correct_state(q_cons_ts(1)%vf, q_prim_vf, pb_ts(1)%sf, mv_ts(1)%sf)
                 else
@@ -802,6 +769,61 @@ contains
 
     end subroutine s_apply_bodyforces
 
+    subroutine s_propagate_immersed_boundaries(s)
+
+        integer, intent(in) :: s
+        integer :: i
+        logical :: forces_computed
+
+        forces_computed = .false.
+
+        do i = 1, num_ibs
+            if (s == 1) then
+                patch_ib(i)%step_vel = patch_ib(i)%vel
+                patch_ib(i)%step_angular_vel = patch_ib(i)%angular_vel
+                patch_ib(i)%step_angles = patch_ib(i)%angles
+                patch_ib(i)%step_x_centroid = patch_ib(i)%x_centroid
+                patch_ib(i)%step_y_centroid = patch_ib(i)%y_centroid
+                patch_ib(i)%step_z_centroid = patch_ib(i)%z_centroid
+            end if
+
+            if (patch_ib(i)%moving_ibm > 0) then
+                patch_ib(i)%vel = (rk_coef(s, 1)*patch_ib(i)%step_vel + rk_coef(s, 2)*patch_ib(i)%vel)/rk_coef(s, 4)
+                patch_ib(i)%angular_vel = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel + rk_coef(s, 2)*patch_ib(i)%angular_vel)/rk_coef(s, 4)
+
+                if (patch_ib(i)%moving_ibm == 1) then
+                    ! plug in analytic velocities for 1-way coupling, if it exists
+                    @:mib_analytical()
+                else if (patch_ib(i)%moving_ibm == 2) then ! if we are using two-way coupling, apply force and torque
+                    ! compute the force and torque on the IB from the fluid
+                    if (.not. forces_computed) then
+                        call s_compute_ib_forces(q_prim_vf, fluid_pp)
+                        forces_computed = .true.
+                    end if
+
+                    ! update the velocity from the force value
+                    patch_ib(i)%vel = patch_ib(i)%vel + rk_coef(s, 3)*dt*(patch_ib(i)%force/patch_ib(i)%mass)/rk_coef(s, 4)
+
+                    ! update the angular velocity with the torque value
+                    patch_ib(i)%angular_vel = (patch_ib(i)%angular_vel*patch_ib(i)%moment) + (rk_coef(s, 3)*dt*patch_ib(i)%torque/rk_coef(s, 4)) ! add the torque to the angular momentum
+                    call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel) ! update the moment of inertia to be based on the direction of the angular momentum
+                    patch_ib(i)%angular_vel = patch_ib(i)%angular_vel/patch_ib(i)%moment ! convert back to angular velocity with the new moment of inertia
+                end if
+
+                ! Update the angle of the IB
+                patch_ib(i)%angles = (rk_coef(s, 1)*patch_ib(i)%step_angles + rk_coef(s, 2)*patch_ib(i)%angles + rk_coef(s, 3)*patch_ib(i)%angular_vel*dt)/rk_coef(s, 4)
+
+                ! Update the position of the IB
+                patch_ib(i)%x_centroid = (rk_coef(s, 1)*patch_ib(i)%step_x_centroid + rk_coef(s, 2)*patch_ib(i)%x_centroid + rk_coef(s, 3)*patch_ib(i)%vel(1)*dt)/rk_coef(s, 4)
+                patch_ib(i)%y_centroid = (rk_coef(s, 1)*patch_ib(i)%step_y_centroid + rk_coef(s, 2)*patch_ib(i)%y_centroid + rk_coef(s, 3)*patch_ib(i)%vel(2)*dt)/rk_coef(s, 4)
+                patch_ib(i)%z_centroid = (rk_coef(s, 1)*patch_ib(i)%step_z_centroid + rk_coef(s, 2)*patch_ib(i)%z_centroid + rk_coef(s, 3)*patch_ib(i)%vel(3)*dt)/rk_coef(s, 4)
+            end if
+        end do
+
+        call s_update_mib(num_ibs, levelset, levelset_norm)
+
+    end subroutine s_propagate_immersed_boundaries
+
     !> This subroutine saves the temporary q_prim_vf vector
         !!      into the q_prim_ts vector that is then used in p_main
         !! @param t_step current time-step
diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py
index 269c8438bf..3125490dd5 100644
--- a/toolchain/mfc/case.py
+++ b/toolchain/mfc/case.py
@@ -7,12 +7,15 @@
 from .state import ARG
 from .run   import case_dicts
 
-
 QPVF_IDX_VARS = {
     'alpha_rho': 'contxb', 'vel'  : 'momxb',         'pres': 'E_idx', 
     'alpha':     'advxb',  'tau_e': 'stress_idx%beg', 'Y':   'chemxb',
     'cf_val': 'c_idx', 'Bx': 'B_idx%beg', 'By': 'B_idx%end-1', 'Bz': 'B_idx%end',
 }
+
+MIBM_ANALYTIC_VARS = [
+    'vel(1)', 'vel(2)', 'vel(3)', 'angular_vel(1)', 'angular_vel(2)', 'angular_vel(3)'
+]
 # "B_idx%end - 1" not "B_idx%beg + 1" must be used because 1D does not have Bx
 
 @dataclasses.dataclass(init=False)
@@ -48,7 +51,7 @@ def get_inp(self, _target) -> str:
         dict_str = ""
         for key, val in self.params.items():
             if key in MASTER_KEYS and key not in case_dicts.IGNORE:
-                if self.__is_ic_analytical(key, val):
+                if self.__is_ic_analytical(key, val) or self.__is_mib_analytical(key, val):
                     dict_str += f"{key} = 0d0\n"
                     ignored.append(key)
                     continue
@@ -95,8 +98,20 @@ def __is_ic_analytical(self, key: str, val: str) -> bool:
 
         return False
 
+    def __is_mib_analytical(self, key: str, val: str) -> bool:
+        '''Is this initial condition analytical?
+        More precisely, is this an arbitrary expression or a string representing a number?'''
+        if common.is_number(val) or not isinstance(val, str):
+            return False
+
+        for variable in MIBM_ANALYTIC_VARS:
+            if re.match(fr'^patch_ib\([0-9]+\)%{re.escape(variable)}', key):
+                return True
+
+        return False
+
     # pylint: disable=too-many-locals
-    def __get_pre_fpp(self, print: bool) -> str:
+    def __get_analytic_ic_fpp(self, print: bool) -> str:
         # generates the content of an FFP file that will hold the functions for
         # some initial condition
         DATA = {
@@ -181,6 +196,72 @@ def rhs_replace(match):
 #:def analytical()
 {f'{chr(10)}{chr(10)}'.join(srcs)}
 #:enddef
+"""
+        return content
+
+    # gets the analytic description of a moving IB's velocity and rotation rate
+    def __get_analytic_mib_fpp(self, print: bool) -> str:
+        # iterates over the parameters and checks if they are defined as an
+        # analytical function. If so, append it to the `patches`` object
+        ib_patches = {}
+
+        for key, val in self.params.items():
+            if not self.__is_mib_analytical(key, val):
+                continue
+
+            patch_id = re.search(r'[0-9]+', key).group(0)
+
+            if patch_id not in ib_patches:
+                ib_patches[patch_id] = []
+
+            ib_patches[patch_id].append((key, val))
+
+        srcs = []
+
+        # for each analytical patch that is required to be added, generate
+        # the string that contains that function.
+        for pid, items in ib_patches.items():
+
+            # function that defines how we will replace variable names with
+            # values from the case file
+            def rhs_replace(match):
+                return {
+                    'x': 'x_cc(i)', 'y': 'y_cc(j)', 'z': 'z_cc(k)',
+                    't': 'mytime',
+
+                    'r':     f'patch_ib({pid})%radius',
+
+                    'e' : f'{math.e}',
+                }.get(match.group(), match.group())
+
+            lines = []
+            # perform the replacement of strings for each analytic function
+            # to generate some fortran string representing the code passed in
+            for attribute, expr in items:
+                if print:
+                    cons.print(f"* Codegen: {attribute} = {expr}")
+
+                lhs = attribute
+                rhs = re.sub(r"[a-zA-Z]+", rhs_replace, expr)
+
+                lines.append(f"        {lhs} = {rhs}")
+
+            # concatenates all of the analytic lines into a single string with
+            # each element separated by new line characters. Then write those
+            # new lines as a fully concatenated string with fortran syntax
+            srcs.append(f"""\
+    if (i == {pid}) then
+{f'{chr(10)}'.join(lines)}
+    end if\
+""")
+
+        content = f"""\
+! This file was generated by MFC. It is only used when we analytically
+! parameterize the velocity and rotation rate of a moving IB.
+
+#:def mib_analytical()
+{f'{chr(10)}{chr(10)}'.join(srcs)}
+#:enddef
 """
         return content
 
@@ -266,9 +347,16 @@ def __get_sim_fpp(self, print: bool) -> str:
         else:
             out = ""
 
+        out = out + self.__get_analytic_mib_fpp(print)
+
         # We need to also include the pre_processing includes so that common subroutines have access to the @:analytical function
         return out + f"\n{self.__get_pre_fpp(print)}"
 
+
+    def __get_pre_fpp(self, print: bool) -> str:
+        out = self.__get_analytic_ic_fpp(print)
+        return out
+
     def get_fpp(self, target, print = True) -> str:
         def _prepend() -> str:
             return f"""\
diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py
index 9240201527..dfdc13b38d 100644
--- a/toolchain/mfc/run/case_dicts.py
+++ b/toolchain/mfc/run/case_dicts.py
@@ -121,9 +121,7 @@ def analytic(self):
         PRE_PROCESS[f"patch_ib({ib_id})%{real_attr}"] = ty
 
     for dir_id in range(1, 4):
-        PRE_PROCESS[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL
         PRE_PROCESS[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL
-        PRE_PROCESS[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL
 
     for cmp_id, cmp in enumerate(["x", "y", "z"]):
         cmp_id += 1
@@ -372,9 +370,9 @@ def analytic(self):
         SIMULATION[f"patch_ib({ib_id})%{real_attr}"] = ty
 
     for dir_id in range(1, 4):
-        SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL
+        SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
         SIMULATION[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL
-        SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL
+        SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL.analytic()
 
     for cmp_id, cmp in enumerate(["x", "y", "z"]):
         cmp_id += 1