diff --git a/examples/data/bem.nc b/examples/data/bem.nc
index 33d70816..bca7066d 100644
Binary files a/examples/data/bem.nc and b/examples/data/bem.nc differ
diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb
index 34c0cd99..01dcc0b9 100644
--- a/examples/tutorial_1_WaveBot.ipynb
+++ b/examples/tutorial_1_WaveBot.ipynb
@@ -97,25 +97,6 @@
     "ndof = fb.nb_dofs"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Next we will add the mass and hydrostatic stiffness properties. \n",
-    "If these values are known they can be added directly.\n",
-    "Here we will use the fact that the WaveBot is free floating and assume constant density to calculate these properties, which Capytaine can natively perform with the `FloatingBody` created above. For convenience, this functionality has been wrapped in `wecopttool.hydrostatics`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "stiffness = wot.hydrostatics.stiffness_matrix(fb).values\n",
-    "mass = wot.hydrostatics.inertia_matrix(fb).values"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -145,7 +126,9 @@
     "\n",
     "We will analyze 50 frequencies with a spacing of 0.05 Hz. These frequencies will be used for the Fourier representation of both the wave and the desired PTO force in the pseudo-spectral problem. See the Theory section of the Documentation for more details on the pseudo-spectral problem formulation.\n",
     "\n",
-    "If you would like to save our BEM data to a NetCDF file for future use, see the `wecopttool.write_netcdf` function."
+    "If you would like to save our BEM data to a NetCDF file for future use, see the `wecopttool.write_netcdf` function.\n",
+    "\n",
+    "The run_bem() function now calculates the hydrostatics and stores them as part of bem_data. The inertia and stiffness can still be defined manually as part of the `FloatingBody` if desired."
    ]
   },
   {
@@ -236,8 +219,6 @@
    "source": [
     "wec = wot.WEC.from_bem(\n",
     "    bem_data,\n",
-    "    inertia_matrix=mass,\n",
-    "    hydrostatic_stiffness=stiffness,\n",
     "    constraints=constraints,\n",
     "    friction=None,\n",
     "    f_add=f_add,\n",
@@ -512,8 +493,6 @@
     "# Update WEC\n",
     "\n",
     "wec_2 = wot.WEC.from_bem(bem_data,\n",
-    "                         inertia_matrix=mass,\n",
-    "                         hydrostatic_stiffness=stiffness,\n",
     "                         constraints=constraints_2,\n",
     "                         friction=None,\n",
     "                         f_add=f_add_2\n",
@@ -558,8 +537,6 @@
    "source": [
     "wec_2_nocon = wot.WEC.from_bem(\n",
     "    bem_data,\n",
-    "    inertia_matrix=mass,\n",
-    "    hydrostatic_stiffness=stiffness,\n",
     "    constraints=None,\n",
     "    friction=None,\n",
     "    f_add=f_add_2)\n",
@@ -712,10 +689,6 @@
     "    nfreq = 50\n",
     "    bem_data = wot.run_bem(fb, freq)\n",
     "\n",
-    "    # Mass & hydrostatic stiffness\n",
-    "    stiffness_3 = wot.hydrostatics.stiffness_matrix(fb).values\n",
-    "    mass_3 = wot.hydrostatics.inertia_matrix(fb).values\n",
-    "\n",
     "    # Impedance definition\n",
     "    omega = bem_data.omega.values\n",
     "    gear_ratio = 12.0\n",
@@ -764,8 +737,6 @@
     "\n",
     "    # Create WEC\n",
     "    wec = wot.WEC.from_bem(bem_data,\n",
-    "            inertia_matrix=mass_3,\n",
-    "            hydrostatic_stiffness=stiffness_3,\n",
     "            constraints=constraints,\n",
     "            friction=None, \n",
     "            f_add=f_add,\n",
@@ -874,7 +845,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.10"
+   "version": "3.10.9"
   },
   "vscode": {
    "interpreter": {
diff --git a/examples/tutorial_2_AquaHarmonics.ipynb b/examples/tutorial_2_AquaHarmonics.ipynb
index dc0c1425..a1548f14 100644
--- a/examples/tutorial_2_AquaHarmonics.ipynb
+++ b/examples/tutorial_2_AquaHarmonics.ipynb
@@ -1,7 +1,6 @@
 {
  "cells": [
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -29,12 +28,15 @@
     "import matplotlib.pyplot as plt\n",
     "from matplotlib import cm\n",
     "from scipy.optimize import brute\n",
+    "import xarray as xr\n",
+    "import logging\n",
+    "logging.basicConfig()\n",
+    "logging.getLogger().setLevel(logging.DEBUG)\n",
     "\n",
     "import wecopttool as wot"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -62,7 +64,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -83,14 +84,12 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "#### Hydrostatics and mass\n",
     "The AquaHarmonics device is positively buoyant (i.e., the buoyancy force is greater than the force due to gravity).\n",
-    "Therefore, we'll calculate the hydrostatic stiffness from the mesh, but set the rigid-body mass manually. \n",
-    "We will also calculate the displaced volume and mass from the mesh, which we will need later for defining forces and constraints."
+    "Therefore, we'll set the rigid-body mass manually, but allow the hydrostatic stiffness to be set automatically by run_bem() based on the mesh. We will also calculate the displaced volume and mass from the mesh (before manually defining the mass of the FloatingBody), which we will need later for defining forces and constraints."
    ]
   },
   {
@@ -99,17 +98,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "mass = np.atleast_2d(5e3) # [kg]\n",
-    "stiffness = wot.hydrostatics.stiffness_matrix(fb).values\n",
-    "\n",
     "g = 9.81\n",
     "rho = 1025\n",
-    "displaced_mass = wot.hydrostatics.inertia_matrix(fb, rho).values  # [kg]\n",
-    "displacement = displaced_mass/rho # [m^3] "
+    "fb.center_of_mass = [0, 0, 0]\n",
+    "fb.rotation_center = fb.center_of_mass\n",
+    "displaced_mass = fb.compute_rigid_body_inertia(rho=rho).values  # [kg]\n",
+    "displacement = displaced_mass/rho # [m^3] \n",
+    "\n",
+    "fb.mass = np.atleast_2d(5e3) # [kg]"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -134,7 +133,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -167,7 +165,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -233,7 +230,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -257,7 +253,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -281,7 +276,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -321,7 +315,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -349,7 +342,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -410,7 +402,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -474,15 +465,13 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "#### WEC object\n",
     "Finally, we can use all the different components we've developed thus far to construct a `WEC` object:\n",
     "\n",
-    " * **BEM data** - defines the hydrodynamics of the hull\n",
-    " * **Inertia and hydrostatics** - rigid-body inertia and hydrostatic stiffness\n",
+    " * **BEM data** - defines the hydrodynamics of the hull and includes hydrostatics\n",
     " * **Constraints** - limitations on the hardware (max power, max torque, etc.) and the constraint to prevent the tether from going slack\n",
     " * **Additional forces** - this captures all of the forces on the WEC that are not due to the interaction between the hull and water (PTO, etc.)"
    ]
@@ -490,20 +479,19 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": false
+   },
    "outputs": [],
    "source": [
     "wec = wot.WEC.from_bem(\n",
     "    bem_data,\n",
-    "    inertia_matrix = mass,\n",
-    "    hydrostatic_stiffness = stiffness,\n",
     "    constraints = constraints,\n",
     "    f_add = f_add,\n",
     ")"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -526,7 +514,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -544,7 +531,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -558,7 +544,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "scale_x_wec = 1e1\n",
@@ -583,7 +571,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -603,7 +590,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -621,7 +607,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -719,7 +704,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -747,7 +731,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -783,12 +766,11 @@
     "    # Unpack geometry variables\n",
     "    mass_ratio = x[0]\n",
     "    mass_var = mass_ratio * max_mass\n",
-    "\n",
+    "    bem_data['inertia_matrix'] = mass_var\n",
+    "    \n",
     "    # update WEC \n",
     "    wec_mass = wot.WEC.from_bem(\n",
     "        bem_data,\n",
-    "        inertia_matrix = mass_var,\n",
-    "        hydrostatic_stiffness = stiffness,\n",
     "        constraints = constraints,\n",
     "        friction = None,\n",
     "        f_add = f_add,\n",
@@ -819,7 +801,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -853,7 +834,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -912,6 +892,13 @@
     "ax3.set_yticks([1, 5, 10, 15, 20, 25])\n",
     "ax3.legend(loc=9)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb
index 52a89035..a57602dc 100644
--- a/examples/tutorial_3_LUPA.ipynb
+++ b/examples/tutorial_3_LUPA.ipynb
@@ -232,7 +232,7 @@
    "source": [
     "#### Combined `FloatingBody`\n",
     "\n",
-    "With both WEC bodies defined separately, we can now create a union of the bodies and define the properties for the overall LUPA device. At the equilibrium position the float is neutrally buoyant while the spar is positively buoyant and requires mooring pre-tension. The combined center of mass and moment of inertia can be found by using the given values weighted by the mass (via the parallel axis theorem for the moment of inertia). This is also when we can specify the surge and pitch degrees of freedom.\n",
+    "With both WEC bodies defined separately, we should first define the respective centers of mass and rotation centers for the `FloatingBody`s. Then, we can create a union of the bodies and define the properties for the overall LUPA device. At the equilibrium position the float is neutrally buoyant while the spar is positively buoyant and requires mooring pre-tension. The combined center of mass and moment of inertia can be found by using the given values weighted by the mass (via the parallel axis theorem for the moment of inertia). This is also when we can specify the surge and pitch degrees of freedom.\n",
     "\n",
     "We are using the density of *fresh* water, $\\rho = 1000 kg/m^3$, since we are modeling LUPA in a wave flume."
    ]
@@ -246,19 +246,23 @@
     "# density of fresh water\n",
     "rho = 1000\n",
     "\n",
-    "# floating body\n",
-    "lupa_fb = float_fb + spar_fb\n",
-    "lupa_fb.name = 'LUPA'\n",
-    "\n",
     "# mass properties float\n",
     "mass_float = float_mass_properties['mass'] \n",
     "cm_float = np.array(float_mass_properties['CG'])\n",
     "pitch_inertia_float = float_mass_properties['MOI'][1] \n",
+    "float_fb.center_of_mass = cm_float\n",
+    "float_fb.rotation_center = float_fb.center_of_mass\n",
     "\n",
     "# mass properties spar\n",
     "mass_spar = spar_mass_properties['mass']  \n",
     "cm_spar = np.array(spar_mass_properties['CG'])\n",
     "pitch_inertia_spar = spar_mass_properties['MOI'][1] \n",
+    "spar_fb.center_of_mass = cm_spar\n",
+    "spar_fb.rotation_center = spar_fb.center_of_mass\n",
+    "\n",
+    "# floating body\n",
+    "lupa_fb = float_fb + spar_fb\n",
+    "lupa_fb.name = 'LUPA'\n",
     "\n",
     " # mass properties LUPA\n",
     "lupa_fb.center_of_mass = ((mass_float*cm_float + mass_spar*cm_spar)\n",
@@ -279,6 +283,36 @@
     "lupa_fb.add_rotation_dof(name='Pitch')"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Define Hydrostatics Manually\n",
+    "\n",
+    "When combining multiple bodies, we need to be careful to set up the hydrostatics correctly. The bodies move separately in heave but move together in surge and pitch. Therefore, we should define the individual heave inertia values for each body and define the total surge and pitch inertia values. The inertia then needs to be reformatted as a DataArray for Capytaine. Lastly, the hydrostatic stiffness can be calculated for the total immersed body."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# reorganize inertia values into DataArray for Capytaine\n",
+    "rigid_inertia_matrix_xr = xr.DataArray(data=np.asarray((inertia)),\n",
+    "                            dims=['influenced_dof', 'radiating_dof'],\n",
+    "                            coords={'influenced_dof': list(lupa_fb.dofs),\n",
+    "                                    'radiating_dof': list(lupa_fb.dofs)},\n",
+    "                            name=\"inertia_matrix\")\n",
+    "\n",
+    "# Set FloatingBody inertia matrix\n",
+    "lupa_fb.inertia_matrix = rigid_inertia_matrix_xr\n",
+    "\n",
+    "# Calculate hydrostatic stiffness after keeping immersed value\n",
+    "lupa_fb = lupa_fb.keep_immersed_part()\n",
+    "lupa_fb.hydrostatic_stiffness = lupa_fb.compute_hydrostatic_stiffness(rho=rho)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -329,12 +363,8 @@
     "try:\n",
     "    bem_data = wot.read_netcdf(filename)\n",
     "except:\n",
-    "    bem_data = wot.run_bem(lupa_fb, freq)\n",
-    "    wot.write_netcdf(filename, bem_data)\n",
-    "\n",
-    "# hydrostatics\n",
-    "_ = lupa_fb.compute_hydrostatics(rho=rho)\n",
-    "hydrostatic_stiffness = lupa_fb.hydrostatic_stiffness"
+    "    bem_data = wot.run_bem(lupa_fb, freq, rho=rho)\n",
+    "    wot.write_netcdf(filename, bem_data)"
    ]
   },
   {
@@ -910,8 +940,6 @@
     "\n",
     "# WEC\n",
     "wec = wot.WEC.from_bem(bem_data,\n",
-    "                       inertia_matrix=inertia,\n",
-    "                       hydrostatic_stiffness=hydrostatic_stiffness,\n",
     "                       constraints=constraints,\n",
     "                       friction=friction,\n",
     "                       f_add=f_add,\n",
@@ -1296,8 +1324,6 @@
     "\n",
     "    # create WEC object\n",
     "    wec = wot.WEC.from_bem(bem_data,\n",
-    "                           inertia_matrix=inertia,\n",
-    "                           hydrostatic_stiffness=hydrostatic_stiffness,\n",
     "                           constraints=constraints,\n",
     "                           friction=friction,\n",
     "                           f_add=f_add,\n",
diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb
index d0a8fde6..a87d36b6 100644
--- a/examples/tutorial_4_Pioneer.ipynb
+++ b/examples/tutorial_4_Pioneer.ipynb
@@ -1,7 +1,6 @@
 {
  "cells": [
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -34,6 +33,7 @@
     "import autograd.numpy as np\n",
     "import matplotlib.pyplot as plt\n",
     "from scipy.linalg import block_diag\n",
+    "import xarray as xr\n",
     "\n",
     "import wecopttool as wot\n",
     "\n",
@@ -42,7 +42,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -85,7 +84,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -108,7 +106,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -137,7 +134,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -178,7 +174,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -197,12 +192,11 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "#### Hydrodynamics and hydrostatics\n",
-    "As mentioned above, the `FloatingBody` object in Capytaine only needs to model the buoy, since no other components are being excited by the waves."
+    "As mentioned above, the `FloatingBody` object in Capytaine only needs to model the buoy, since no other components are being excited by the waves. The inertia matrix needs to be defined manually here since it is not based on the displaced mass."
    ]
   },
   {
@@ -216,29 +210,36 @@
     "pnr_fb.center_of_mass = np.array([0., 0., buoy_props['CG']])\n",
     "pnr_fb.rotation_center = pnr_fb.center_of_mass\n",
     "ndof = pnr_fb.nb_dofs\n",
-    "pnr_fb.show_matplotlib()"
+    "pnr_fb.show_matplotlib()\n",
+    "\n",
+    "pnr_fb.inertia_matrix = xr.DataArray([[buoy_props['MOI']]],\n",
+    "                            dims=['influenced_dof', 'radiating_dof'],\n",
+    "                            coords={'influenced_dof' : ['Pitch'],\n",
+    "                            'radiating_dof' : ['Pitch']},\n",
+    "                            name=\"inertia_matrix\"\n",
+    "                            )"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "rho = 1025. # kg/m^3\n",
     "freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n",
-    "bem_data = wot.run_bem(pnr_fb, freq)\n",
+    "bem_data = wot.run_bem(pnr_fb, freq, rho=rho)\n",
     "omega = bem_data.omega.values\n",
     "\n",
     "pnr_fb.keep_immersed_part()\n",
-    "k_buoy = pnr_fb.compute_hydrostatic_stiffness(rho=rho).values.squeeze()\n",
     "k_spring = spring_props['Max torque'] / spring_props['Max displacement']\n",
-    "print(f'Hydrostatic stiffness from Capytaine: {k_buoy} N-m/rad')\n",
+    "print(f\"Hydrostatic stiffness from Capytaine: {bem_data['hydrostatic_stiffness'].values.squeeze()} N-m/rad\")\n",
     "print('Hydrostatic stiffness from experiment: 37204 N-m/rad')"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -280,7 +281,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -292,7 +292,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -316,7 +315,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -351,7 +349,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -372,7 +369,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -397,7 +393,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -439,7 +434,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -484,7 +478,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -506,7 +499,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -537,7 +529,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -569,7 +560,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -595,7 +585,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -610,8 +599,6 @@
    "outputs": [],
    "source": [
     "wec = wot.WEC.from_bem(bem_data,\n",
-    "                       inertia_matrix=np.array([[buoy_props['MOI']]]),\n",
-    "                       hydrostatic_stiffness=k_buoy,\n",
     "                       f_add=f_add,\n",
     "                       constraints=constraints,\n",
     "                       uniform_shift=False,\n",
@@ -619,7 +606,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -647,7 +633,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -678,7 +663,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -745,7 +729,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -753,7 +736,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -767,7 +749,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hydro_data = wot.linear_hydrodynamics(bem_data, np.array([[buoy_props['MOI']]]), k_buoy, friction=None)    \n",
+    "hydro_data = wot.add_linear_friction(bem_data, friction=None)    \n",
     "hydro_data = wot.check_linear_damping(hydro_data, uniform_shift=False)    \n",
     "\n",
     "Zi = wot.hydrodynamic_impedance(hydro_data)\n",
@@ -784,7 +766,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -856,7 +837,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -886,7 +866,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -968,11 +947,18 @@
     "    axi.label_outer()\n",
     "    axi.autoscale(axis='x', tight=True)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "conda-wot",
+   "display_name": "wot_dev",
    "language": "python",
    "name": "python3"
   },
@@ -986,7 +972,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.10"
+   "version": "3.10.9"
   }
  },
  "nbformat": 4,
diff --git a/tests/test_core.py b/tests/test_core.py
index 6695357d..f372d760 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -55,17 +55,22 @@ def bem_data(f1, nfreq):
     ndof = 2; ndir = 3;
     radiation_dims = ['radiating_dof', 'influenced_dof', 'omega']
     excitation_dims = ['influenced_dof', 'wave_direction', 'omega']
+    hydrostatics_dims = ['radiating_dof', 'influenced_dof']
 
     added_mass = np.ones([ndof, ndof, nfreq])
     radiation_damping = np.ones([ndof, ndof, nfreq])
     diffraction_force = np.ones([ndof, ndir, nfreq], dtype=complex) + 1j
     Froude_Krylov_force = np.ones([ndof, ndir, nfreq], dtype=complex) + 1j
+    inertia_matrix = np.ones([ndof, ndof])
+    hydrostatic_stiffness = np.ones([ndof, ndof])
 
     data_vars = {
         'added_mass': (radiation_dims, added_mass),
         'radiation_damping': (radiation_dims, radiation_damping),
         'diffraction_force': (excitation_dims, diffraction_force),
-        'Froude_Krylov_force': (excitation_dims, Froude_Krylov_force)
+        'Froude_Krylov_force': (excitation_dims, Froude_Krylov_force),
+        'inertia_matrix': (hydrostatics_dims, inertia_matrix),
+        'hydrostatic_stiffness': (hydrostatics_dims, hydrostatic_stiffness)
     }
     return xr.Dataset(data_vars=data_vars, coords=coords)
 
@@ -75,11 +80,9 @@ def hydro_data(bem_data):
     """Synthetic hydro-data containing inertia, stiffness, and friction
     in addition to the coefficients in `bem_data`."""
     ndof = len(bem_data.influenced_dof)
-    inertia_matrix = np.ones([ndof, ndof])
-    stiffness = np.ones([ndof, ndof])
     friction = np.ones([ndof, ndof])
-    data = wot.linear_hydrodynamics(
-        bem_data, inertia_matrix, stiffness, friction
+    data = wot.add_linear_friction(
+        bem_data, friction
     )
     return data
 
@@ -1300,12 +1303,12 @@ def test_round_trip(self, bem_data, data_back):
 
 
 class TestLinearHydrodynamics:
-    """Test function :python:`linear_hydrodynamics`."""
+    """Test function :python:`add_linear_friction`."""
 
     def test_values(self, bem_data, hydro_data):
         """Test the function returns expected values."""
         mat = np.array([[1, 1], [1, 1]])
-        calculated = wot.linear_hydrodynamics(bem_data, mat, mat, mat)
+        calculated = wot.add_linear_friction(bem_data, mat)
         xr.testing.assert_allclose(calculated, hydro_data)
 
 
diff --git a/tests/test_hydrostatics.py b/tests/test_hydrostatics.py
deleted file mode 100644
index d4178cf3..00000000
--- a/tests/test_hydrostatics.py
+++ /dev/null
@@ -1,204 +0,0 @@
-""" Unit tests for functions in the :python:`hydrostatics.py` module.
-"""
-
-import pytest
-import capytaine as cpy
-import numpy as np
-
-import wecopttool as wot
-from wecopttool.core import _default_parameters
-
-
-# setup simple constant density rectangular barge
-@pytest.fixture()
-def rho():
-    """Water density [kg/m^3]."""
-    return _default_parameters['rho']
-
-
-@pytest.fixture()
-def g():
-    """Gravitational acceleration [m/s^2]."""
-    return _default_parameters['g']
-
-
-@pytest.fixture()
-def lx():
-    """Barge length [m]."""
-    return 4.0
-
-
-@pytest.fixture()
-def ly():
-    """Barge width [m]."""
-    return 5.0
-
-
-@pytest.fixture()
-def lz():
-    """Barge height [m]."""
-    return 2.0
-
-
-@pytest.fixture()
-def cog():
-    """Center of gravity of the barge."""
-    return (0.0, 0.0, 0.0)
-
-
-@pytest.fixture()
-def mass(lx, ly, lz, rho):
-    """Mass of the barge."""
-    return lx*ly*lz/2 * rho
-
-
-@pytest.fixture()
-def fb(lx, ly, lz):
-    """Simple constant density rectangular barge."""
-    rect = cpy.RectangularParallelepiped(
-        (lx, ly, lz), resolution=(100, 100, 10), center=(0.0, 0.0, 0.0,))
-    rect.add_all_rigid_body_dofs()
-    rect.center_of_mass = None
-    fb.mass = None
-    return rect
-
-
-@pytest.fixture()
-def fb_mass(fb, cog, mass):
-    """Simple constant density rectangular barge with the mass
-    properties specified.
-    """
-    fbc = fb.copy()
-    fbc.center_of_mass = cog
-    fbc.mass = mass
-    return fbc
-
-
-@pytest.fixture()
-def stiffness(lx, ly, lz, rho, g):
-    """True (theoretical/calculated) value of the stiffness matrix for
-    the barge.
-    """
-    stiffness = np.zeros([6, 6])
-    stiffness[2, 2] = rho*g * lx*ly
-    stiffness[2, 3] = 0.0
-    stiffness[2, 4] = 0.0
-    stiffness[3, 3] = (rho*g * lx*ly**3 /12) + rho*g*(lx*ly*lz/2)*(-lz/4)
-    stiffness[3, 4] = 0.0
-    stiffness[3, 5] = 0.0
-    stiffness[4, 4] = (rho*g * ly*lx**3 /12) + rho*g*(lx*ly*lz/2)*(-lz/4)
-    stiffness[4, 5] = 0.0
-    stiffness[3, 2] = stiffness[2, 3]
-    stiffness[4, 2] = stiffness[2, 4]
-    stiffness[4, 3] = stiffness[3, 4]
-    return stiffness
-
-
-@pytest.fixture()
-def inertia(lx, ly, lz, mass):
-    """True (theoretical/calculated) value of the inertia matrix for
-    the barge.
-    """
-    inertia = np.zeros([6, 6])
-    inertia[0, 0] = mass
-    inertia[1, 1] = mass
-    inertia[2, 2] = mass
-    inertia[3, 3] = mass/12 * (ly**2 + lz**2)
-    inertia[3, 4] = 0.0
-    inertia[3, 5] = 0.0
-    inertia[4, 3] = 0.0
-    inertia[4, 4] = mass/12 * (lx**2 + lz**2)
-    inertia[4, 5] = 0.0
-    inertia[5, 3] = 0.0
-    inertia[5, 4] = 0.0
-    inertia[5, 5] = mass/12 * (lx**2 + ly**2)
-    return inertia
-
-
-class TestStiffnessMatrix:
-    """Test function :python:`hydrostatics.stiffness_matrix`."""
-
-    def test_inferred_cog(self, fb, rho, g, stiffness):
-        """Test the function with the center of gravity not provided,
-        but inferred.
-        """
-        stiffness_calc = wot.hydrostatics.stiffness_matrix(fb, rho, g)
-        assert np.allclose(stiffness, stiffness_calc, rtol=0.01)
-
-
-    def test_given_cog(self, fb, rho, g, cog, stiffness):
-        """Test the function with the center of gravity provided as a
-        function input and not in the floating body.
-        """
-        stiffness_calc = wot.hydrostatics.stiffness_matrix(fb, rho, g, cog)
-        assert np.allclose(stiffness, stiffness_calc, rtol=0.01)
-
-
-    def test_given_cog_fb(self, fb_mass, rho, g, stiffness):
-        """Test the function with the center of gravity provided in the
-        floating body and not as a function input.
-        """
-        stiffness_calc = wot.hydrostatics.stiffness_matrix(fb_mass, rho, g)
-        assert np.allclose(stiffness, stiffness_calc, rtol=0.01)
-
-
-    def test_given_cog_redundant(self, fb_mass, rho, g, cog, stiffness):
-        """Test the function with the center of gravity provided in both
-        the floating body and the function inputs.
-        """
-        stiffness_calc = wot.hydrostatics.stiffness_matrix(fb_mass, rho, g, cog)
-        assert np.allclose(stiffness, stiffness_calc, rtol=0.01)
-
-
-    def test_cog_mismatch(self, fb_mass, rho, g):
-        """Test that the function fails if the center of gravity is
-        provided in both the floating body and the function inputs but
-        have different values.
-        """
-        cog_wrong = (0, 0, -0.1)
-        with pytest.raises(ValueError):
-            wot.hydrostatics.stiffness_matrix(fb_mass, rho, g, cog_wrong)
-
-
-class TestInertiaMatrix:
-    """Test function :python:`hydrostatics.test_inertia_matrix`."""
-
-    def test_inferred_mass(self, fb, rho, inertia):
-        """Test the function with the mass not provided, but inferred.
-        """
-        inertia_calc = wot.hydrostatics.inertia_matrix(fb, rho)
-        assert np.allclose(inertia, inertia_calc, rtol=0.01)
-
-
-    def test_given_mass(self, fb, rho, cog, mass, inertia):
-        """Test the function with the mass provided as a function input
-        and not in the floating body.
-        """
-        inertia_calc = wot.hydrostatics.inertia_matrix(fb, rho, cog, mass)
-        assert np.allclose(inertia, inertia_calc, rtol=0.01)
-
-
-    def test_given_mass_fb(self, fb_mass, rho, cog, mass, inertia):
-        """Test the function with the mass provided in the floating body
-        and not as a function input.
-        """
-        inertia_calc = wot.hydrostatics.inertia_matrix(fb_mass, rho)
-        assert np.allclose(inertia, inertia_calc, rtol=0.01)
-
-
-    def test_given_mass_redundant(self, fb_mass, rho, cog, mass, inertia):
-        """Test the function with the mass provided in both the floating
-        body and the function inputs.
-        """
-        inertia_calc = wot.hydrostatics.inertia_matrix(fb_mass, rho, cog, mass)
-        assert np.allclose(inertia, inertia_calc, rtol=0.01)
-
-
-    def test_mass_mismatch(self, fb_mass, rho, cog, mass):
-        """Test that the function fails if the mass is provided in both
-        the floating body and the function inputs but have different
-        values.
-        """
-        mass_wrong = 0.9*mass
-        with pytest.raises(ValueError):
-            wot.hydrostatics.inertia_matrix(fb_mass, rho, cog, mass_wrong)
diff --git a/tests/test_integration.py b/tests/test_integration.py
index 46ed235e..a93c178e 100644
--- a/tests/test_integration.py
+++ b/tests/test_integration.py
@@ -104,20 +104,16 @@ def irregular_wave(f1, nfreq):
 @pytest.fixture(scope='module')
 def wec_from_bem(f1, nfreq, bem, fb, pto):
     """Simple WEC: 1 DOF, no constraints."""
-    mass = wot.hydrostatics.inertia_matrix(fb).values
-    hstiff = wot.hydrostatics.stiffness_matrix(fb).values
     f_add = {"PTO": pto.force_on_wec}
-    wec = wot.WEC.from_bem(bem, mass, hstiff, f_add=f_add)
+    wec = wot.WEC.from_bem(bem, f_add=f_add)
     return wec
 
 
 @pytest.fixture(scope='module')
 def wec_from_floatingbody(f1, nfreq, fb, pto):
     """Simple WEC: 1 DOF, no constraints."""
-    mass = wot.hydrostatics.inertia_matrix(fb).values
-    hstiff = wot.hydrostatics.stiffness_matrix(fb).values
     f_add = {"PTO": pto.force_on_wec}
-    wec = wot.WEC.from_floating_body(fb, f1, nfreq, mass, hstiff, f_add=f_add)
+    wec = wot.WEC.from_floating_body(fb, f1, nfreq, f_add=f_add)
     return wec
 
 
@@ -130,10 +126,13 @@ def wec_from_impedance(bem, pto, fb):
     w = np.expand_dims(omega, [0, 1])
     A = bemc['added_mass'].values
     B = bemc['radiation_damping'].values
-    mass = wot.hydrostatics.inertia_matrix(fb).values
-    hstiff = wot.hydrostatics.stiffness_matrix(fb).values
+    fb.center_of_mass = [0, 0, 0]
+    fb.rotation_center = fb.center_of_mass
+    fb = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part()
+    mass = bemc['inertia_matrix'].values
+    hstiff = bemc['hydrostatic_stiffness'].values
     K = np.expand_dims(hstiff, 2)
-
+    
     freqs = omega / (2 * np.pi)
     impedance = (A + mass)*(1j*w) + B + K/(1j*w)
     exc_coeff = bem['Froude_Krylov_force'] + bem['diffraction_force']
@@ -146,9 +145,7 @@ def wec_from_impedance(bem, pto, fb):
 @pytest.fixture(scope='module')
 def resonant_wave(f1, nfreq, fb, bem):
     """Regular wave at natural frequency of the WEC"""
-    mass = wot.hydrostatics.inertia_matrix(fb).values
-    hstiff = wot.hydrostatics.stiffness_matrix(fb).values
-    hd = wot.linear_hydrodynamics(bem, mass, hstiff)
+    hd = wot.add_linear_friction(bem)
     Zi = wot.hydrodynamic_impedance(hd)
     wn = Zi['omega'][np.abs(Zi).argmin(dim='omega')].item()
     waves = wot.waves.regular_wave(f1, nfreq, freq=wn/2/np.pi, amplitude=0.1)
@@ -215,7 +212,7 @@ def test_same_wec_init(wec_from_bem,
     bem_res = wec_from_bem.residual(x_wec_0, x_opt_0, waves)
     fb_res = wec_from_floatingbody.residual(x_wec_0, x_opt_0, waves)
     imp_res = wec_from_impedance.residual(x_wec_0, x_opt_0, waves)
-
+    
     assert fb_res == approx(bem_res, rel=0.01)
     assert imp_res == approx(bem_res, rel=0.01)
 
@@ -227,17 +224,17 @@ class TestTheoreticalPowerLimits:
     @pytest.fixture(scope='class')
     def mass(self, fb):
         """Rigid-body mass"""
-        return wot.hydrostatics.inertia_matrix(fb).values
+        return fb.compute_rigid_body_inertia()
 
     @pytest.fixture(scope='class')
     def hstiff(self, fb):
         """Hydrostatic stiffness"""
-        return wot.hydrostatics.stiffness_matrix(fb).values
+        return fb.compute_hydrostatic_stiffness()
 
     @pytest.fixture(scope='class')
-    def hydro_impedance(self, bem, mass, hstiff):
+    def hydro_impedance(self, bem):
         """Intrinsic hydrodynamic impedance"""
-        hd = wot.linear_hydrodynamics(bem, mass, hstiff)
+        hd = wot.add_linear_friction(bem)
         hd = wot.check_linear_damping(hd)
         Zi = wot.hydrodynamic_impedance(hd)
         return Zi
@@ -246,14 +243,12 @@ def test_p_controller_resonant_wave(self,
                                         bem,
                                         resonant_wave,
                                         p_controller_pto,
-                                        mass,
-                                        hstiff,
                                         hydro_impedance):
         """Proportional controller should match optimum for natural resonant
         wave"""
-
+        
         f_add = {"PTO": p_controller_pto.force_on_wec}
-        wec = wot.WEC.from_bem(bem, mass, hstiff, f_add=f_add)
+        wec = wot.WEC.from_bem(bem, f_add=f_add)
 
         res = wec.solve(waves=resonant_wave,
                         obj_fun=p_controller_pto.average_power,
@@ -275,19 +270,17 @@ def test_p_controller_resonant_wave(self,
         power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())
                          ).squeeze().sum('omega').item()
 
-        assert power_sol == approx(power_optimal, rel=0.02)
+        assert power_sol == approx(power_optimal, rel=0.03)
 
     def test_pi_controller_regular_wave(self,
                                         bem,
                                         regular_wave,
                                         pi_controller_pto,
-                                        mass,
-                                        hstiff,
                                         hydro_impedance):
         """PI controller matches optimal for any regular wave"""
 
         f_add = {"PTO": pi_controller_pto.force_on_wec}
-        wec = wot.WEC.from_bem(bem, mass, hstiff, f_add=f_add)
+        wec = wot.WEC.from_bem(bem, f_add=f_add)
 
         res = wec.solve(waves=regular_wave,
                         obj_fun=pi_controller_pto.average_power,
@@ -317,14 +310,12 @@ def test_unstructured_controller_irregular_wave(self,
                                                     regular_wave,
                                                     pto,
                                                     nfreq,
-                                                    mass,
-                                                    hstiff,
                                                     hydro_impedance):
         """Unstructured (numerical optimal) controller matches optimal for any
         irregular wave when unconstrained"""
 
         f_add = {"PTO": pto.force_on_wec}
-        wec = wot.WEC.from_bem(bem, mass, hstiff, f_add=f_add)
+        wec = wot.WEC.from_bem(bem, f_add=f_add)
 
         res = wec.solve(waves=regular_wave,
                         obj_fun=pto.average_power,
@@ -343,4 +334,4 @@ def test_unstructured_controller_irregular_wave(self,
         power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())
                          ).squeeze().sum('omega').item()
 
-        assert power_sol == approx(power_optimal, rel=1e-3)
+        assert power_sol == approx(power_optimal, rel=1e-2)
diff --git a/wecopttool/__init__.py b/wecopttool/__init__.py
index d36607f3..9b3d6b66 100644
--- a/wecopttool/__init__.py
+++ b/wecopttool/__init__.py
@@ -31,7 +31,6 @@
 
 from wecopttool.core import *
 from wecopttool import waves
-from wecopttool import hydrostatics
 from wecopttool import pto
 
 try:
diff --git a/wecopttool/core.py b/wecopttool/core.py
index e60f14a2..1ea44943 100644
--- a/wecopttool/core.py
+++ b/wecopttool/core.py
@@ -40,7 +40,7 @@
     "standard_forces",
     "run_bem",
     "change_bem_convention",
-    "linear_hydrodynamics",
+    "add_linear_friction",
     "wave_excitation",
     "hydrodynamic_impedance",
     "atleast_2d",
@@ -50,6 +50,7 @@
     "decompose_state",
     "frequency_parameters",
     "time_results",
+    "set_fb_centers",
 ]
 
 
@@ -284,8 +285,6 @@ def __repr__(self) -> str:
     @staticmethod
     def from_bem(
         bem_data: Union[Dataset, Union[str, Path]],
-        inertia_matrix: Optional[ndarray] = None,
-        hydrostatic_stiffness: Optional[ndarray] = None,
         friction: Optional[ndarray] = None,
         f_add: Optional[TIForceDict] = None,
         constraints: Optional[Iterable[Mapping]] = None,
@@ -312,26 +311,16 @@ def from_bem(
         :py:func:`wecopttool.run_bem`,
         rather than running Capytaine directly, which outputs the
         results in the correct convention. The results can be saved
-        using :py:func:`wecopttool.write_netcdf`.
-
-        In addition to the Capytaine results, if the dataset contains
-        the :python:`inertia_matrix`, :python:`hydrostatic_stiffness`,
-        or :python:`friction` these do not need to be provided
-        separately.
+        using :py:func:`wecopttool.write_netcdf`. 
+        :py:func:`wecopttool.run_bem` also computes the inertia and 
+        hydrostatic stiffness which should be included in bem_data.
 
         Parameters
         ----------
         bem_data
             Linear hydrodynamic coefficients obtained using the boundary
             element method (BEM) code Capytaine, with sign convention
-            corrected.
-        inertia_matrix
-           Inertia matrix of size :python:`(ndof, ndof)`.
-           :python:`None` if included in :python:`bem_data`.
-        hydrostatic_stiffness
-            Linear hydrostatic restoring coefficient of size
-            :python:`(nodf, ndof)`.
-            :python:`None` if included in :python:`bem_data`.
+            corrected. Also includes inertia and hydrostatic stiffness.
         friction
             Linear friction, in addition to radiation damping, of size
             :python:`(nodf, ndof)`.
@@ -362,32 +351,16 @@ def from_bem(
             If :python:`None` the names
             :python:`['DOF_0', ..., 'DOF_N']` are used.
 
-        Raises
-        ------
-        ValueError
-            If either :python:`inertia_matrix` or
-            :python:`hydrostatic_stiffness` are :python:`None` and is
-            not included in :python:`bem_data`.
-            See :py:func:`wecopttool.linear_hydrodynamics`.
-        ValueError
-            If any of :python:`inertia_matrix`,
-            :python:`hydrostatic_stiffness`, or :python:`stiffness` are
-            both provided and included in :python:`bem_data` but have
-            different values.
-            See :py:func:`wecopttool.linear_hydrodynamics`.
-
         See Also
         --------
-        run_bem, linear_hydrodynamics, change_bem_convention,
+        run_bem, add_linear_friction, change_bem_convention,
         write_netcdf, check_linear_damping
         """
         if isinstance(bem_data, (str, Path)):
             bem_data = read_netcdf(bem_data)
-        # add inertia_matrix, hydrostatic stiffness, and friction
-        hydro_data = linear_hydrodynamics(
-            bem_data, inertia_matrix, hydrostatic_stiffness, friction)
-        if inertia_matrix is None:
-            inertia_matrix = hydro_data['inertia_matrix'].values
+        # add friction
+        hydro_data = add_linear_friction(bem_data, friction)
+        inertia_matrix = hydro_data['inertia_matrix'].values
 
         # frequency array
         f1, nfreq = frequency_parameters(
@@ -411,8 +384,6 @@ def from_floating_body(
         fb: cpy.FloatingBody,
         f1: float,
         nfreq: int,
-        inertia_matrix: ndarray,
-        hydrostatic_stiffness: ndarray,
         friction: Optional[ndarray] = None,
         f_add: Optional[TIForceDict] = None,
         constraints: Optional[Iterable[Mapping]] = None,
@@ -454,11 +425,6 @@ def from_floating_body(
         nfreq
             Number of frequencies (not including zero frequency),
             i.e., :python:`freqs = [0, f1, 2*f1, ..., nfreq*f1]`.
-        inertia_matrix
-           Inertia matrix of size :python:`(ndof, ndof)`.
-        hydrostatic_stiffness
-            Linear hydrostatic restoring coefficient of size
-            :python:`(ndof, ndof)`.
         friction
             Linear friction, in addition to radiation damping, of size
             :python:`(ndof, ndof)`.
@@ -504,7 +470,7 @@ def from_floating_body(
         bem_data = run_bem(
             fb, freq, wave_directions, rho=rho, g=g, depth=depth)
         wec = WEC.from_bem(
-            bem_data, inertia_matrix, hydrostatic_stiffness, friction, f_add,
+            bem_data, friction, f_add,
             constraints, min_damping=min_damping)
         return wec
 
@@ -2157,12 +2123,17 @@ def run_bem(
         # radiation only problem, no diffraction or excitation
         test_matrix = test_matrix.drop_vars('wave_direction')
     if write_info is None:
-        write_info = {'hydrostatics': False,
+        write_info = {'hydrostatics': True,
                       'mesh': False,
                       'wavelength': False,
                       'wavenumber': False,
                      }
     wec_im = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part()
+    wec_im = set_fb_centers(wec_im, rho=rho)
+    if not hasattr(wec_im, 'inertia_matrix'):
+        wec_im.inertia_matrix = wec_im.compute_rigid_body_inertia(rho=rho)
+    if not hasattr(wec_im, 'hydrostatic_stiffness'):
+        wec_im.hydrostatic_stiffness = wec_im.compute_hydrostatic_stiffness(rho=rho, g=g)
     bem_data = solver.fill_dataset(
         test_matrix, wec_im, n_jobs=njobs, **write_info)
     return change_bem_convention(bem_data)
@@ -2189,14 +2160,11 @@ def change_bem_convention(bem_data: Dataset) -> Dataset:
     return bem_data
 
 
-def linear_hydrodynamics(
+def add_linear_friction(
     bem_data: Dataset,
-    inertia_matrix: Optional[ArrayLike] = None,
-    hydrostatic_stiffness: Optional[ArrayLike] = None,
     friction: Optional[ArrayLike] = None
 ) -> Dataset:
-    """Add rigid body inertia_matrix, hydrostatic stiffness, and linear
-    friction to BEM data.
+    """Add linear friction to BEM data.
 
     Returns the Dataset with the additional information added.
 
@@ -2205,62 +2173,32 @@ def linear_hydrodynamics(
     bem_data
         Linear hydrodynamic coefficients obtained using the boundary
         element method (BEM) code Capytaine, with sign convention
-        corrected.
-    inertia_matrix
-        Inertia matrix of size :python:`(ndof, ndof)`.
-        :python:`None` if included in :python:`bem_data`.
-    hydrostatic_stiffness
-        Linear hydrostatic restoring coefficient of size
-        :python:`(nodf, ndof)`.
-        :python:`None` if included in :python:`bem_data`.
+        corrected. Also includes inertia and hydrostatic stiffness.
     friction
         Linear friction, in addition to radiation damping, of size
         :python:`(nodf, ndof)`.
         :python:`None` if included in :python:`bem_data` or to set to zero.
-
-    Raises
-    ------
-    ValueError
-        If either :python:`inertia_matrix` or
-        :python:`hydrostatic_stiffness` are :python:`None` and is not
-        included in :python:`bem_data`.
-    ValueError
-        If any of :python:`inertia_matrix`,
-        :python:`hydrostatic_stiffness`, or :python:`friction` are both
-        provided and included in :python:`bem_data` but have different
-        values.
     """
-    vars = {'inertia_matrix': inertia_matrix, 'friction': friction,
-            'hydrostatic_stiffness': hydrostatic_stiffness}
-
     dims = ['radiating_dof', 'influenced_dof']
-
     hydro_data = bem_data.copy(deep=True)
-
-    for name, data in vars.items():
-        org = name in hydro_data.variables.keys()
-        new = data is not None
-        if new and org:
+    
+    if friction is not None:
+        if 'friction' in hydro_data.variables.keys():
             if not np.allclose(data, hydro_data.variables[name]):
                 raise ValueError(
-                    f'BEM data already has variable "{name}" ' +
-                    'with diferent values')
-            else :
+                        f'Variable "friction" is already in BEM data ' +
+                        f'with different values.')
+            else:
                 _log.warning(
                     f'Variable "{name}" is already in BEM data ' +
                     'with same value.')
-        elif (not new) and (not org):
-            if name=='friction':
-                ndof = len(hydro_data["influenced_dof"])
-                hydro_data[name] = (dims, np.zeros([ndof, ndof]))
-            else:
-                raise ValueError(
-                    f'Variable "{name}" is not in BEM data and ' +
-                    'was not provided.')
-        elif new:
-            data = atleast_2d(data)
-            hydro_data[name] = (dims, data)
-
+        else:
+            data = atleast_2d(friction)
+            hydro_data['friction'] = (dims, friction)
+    elif friction is None:
+        ndof = len(hydro_data["influenced_dof"])
+        hydro_data['friction'] = (dims, np.zeros([ndof, ndof]))
+            
     return hydro_data
 
 
@@ -2320,7 +2258,7 @@ def hydrodynamic_impedance(hydro_data: Dataset) -> Dataset:
     ----------
     hydro_data
         Dataset with linear hydrodynamic coefficients produced by
-        :py:func:`wecopttool.linear_hydrodynamics`.
+        :py:func:`wecopttool.add_linear_friction`.
     """
 
     Zi = (hydro_data['inertia_matrix'] \
@@ -2569,3 +2507,36 @@ def time_results(fd: DataArray, time: DataArray) -> ndarray:
         out = out + \
             np.real(mag)*np.cos(w*time) - np.imag(mag)*np.sin(w*time)
     return out
+
+
+def set_fb_centers(
+    fb: FloatingBody,
+    rho: float = _default_parameters["rho"],
+) -> FloatingBody:
+    """Sets default properties if not provided by the user:
+        - `center_of_mass` is set to the geometric centroid
+        - `rotation_center` is set to the center of mass
+    """
+    valid_properties = ['center_of_mass', 'rotation_center']
+    
+    for property in valid_properties:
+        if not hasattr(fb, property):
+            setattr(fb, property, None)
+        if getattr(fb, property) is None:
+            if property == 'center_of_mass':
+                def_val = fb.center_of_buoyancy
+                log_str = (
+                    "Using the geometric centroid as the center of gravity (COG).")
+            elif property == 'rotation_center':
+                def_val = fb.center_of_mass
+                log_str = (
+                    "Using the center of gravity (COG) as the rotation center " +
+                    "for hydrostatics.")
+            setattr(fb, property, def_val)
+            _log.warning(log_str)
+        elif getattr(fb, property) is not None:
+            _log.warning(
+                f'{property} already defined as {getattr(fb, property)}.')
+            
+    return fb
+
diff --git a/wecopttool/hydrostatics.py b/wecopttool/hydrostatics.py
deleted file mode 100644
index 29454ddd..00000000
--- a/wecopttool/hydrostatics.py
+++ /dev/null
@@ -1,169 +0,0 @@
-"""Functions for calculating hydrostatic and mass properties for
-floating bodies.
-"""
-
-
-from __future__ import annotations
-
-
-__all__ = [
-    "stiffness_matrix",
-    "inertia_matrix"
-]
-
-
-from typing import Iterable, Optional
-import logging
-
-import numpy as np
-from capytaine import FloatingBody
-from xarray import DataArray
-
-from wecopttool.core import _default_parameters
-
-
-# logger
-_log = logging.getLogger(__name__)
-
-def stiffness_matrix(
-    fb: FloatingBody,
-    rho: float = _default_parameters['rho'],
-    g: float = _default_parameters['g'],
-    center_of_mass: Optional[Iterable[float]] = None,
-    rotation_center: Optional[Iterable[float]] = None
-) -> DataArray:
-    """Compute the hydrostatic stiffness of a Capytaine floating body.
-
-    .. note:: Only works for rigid body DOFs which must be named
-              according to the Capytaine convention (e.g.,
-              :python:`"Heave"`).
-
-    Uses
-    :py:meth:`capytaine.bodies.bodies.FloatingBody.compute_hydrostatic_stiffness`
-    on the immersed part of the mesh.
-
-    Parameters
-    ----------
-    fb
-        A capytaine floating body (mesh + DOFs, and optionally center of
-        mass).
-    rho
-        Water density in :math:`kg/m^3`.
-    g
-        Gravitational acceleration in :math:`m/s^2`.
-    center_of_mass
-        Center of gravity/mass :python:`(cx, cy, cz)`.
-    rotation_center
-        Center of rotation :python:`(rx, ry, rz)` for hydrostatics
-        calculations.
-
-    Raises
-    ------
-    ValueError
-        If :python:`fb.center_of_mass is not None` and
-        :python:`center_of_mass` is provided with a different value.
-    """
-    fb = _set_property(fb, 'center_of_mass', center_of_mass)
-    fb = _set_property(fb, 'rotation_center', rotation_center)
-    fb_im = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part()
-    return fb_im.compute_hydrostatic_stiffness(rho=rho, g=g)
-
-
-def inertia_matrix(
-    fb: FloatingBody,
-    rho: Optional[float] = _default_parameters['rho'],
-    center_of_mass: Optional[Iterable[float]] = None,
-    mass: Optional[float] = None,
-    rotation_center: Optional[Iterable[float]] = None,
-) -> DataArray:
-    """Compute the inertia (mass) matrix assuming a constant density for
-    the WEC.
-
-    .. note:: This function assumes a constant density WEC.
-
-    Uses
-    :py:meth:`capytaine.bodies.bodies.FloatingBody.compute_rigid_body_inertia`
-    on the full mesh.
-
-    Parameters
-    ----------
-    fb
-        A capytaine floating body (mesh + DOFs, and optionally center of
-        mass and mass).
-    rho
-        Water density in :math:`kg/m^3`.
-    center_of_mass
-        Center of gravity/mass.
-    mass
-        Rigid body mass.
-    rotation_center
-        Center of rotation :python:`(rx, ry, rz)` for hydrostatics
-        calculations.
-
-    Raises
-    ------
-    ValueError
-        If :python:`fb.center_of_mass is not None` and
-        :python:`center_of_mass` is provided with a different value.
-    ValueError
-        If :python:`fb.mass is not None` and :python:`mass` is provided
-        with a different value.
-    """
-    fb = _set_property(fb, 'center_of_mass', center_of_mass)
-    fb = _set_property(fb, 'rotation_center', rotation_center)
-    fb = _set_property(fb, 'mass', mass, rho)
-    return fb.compute_rigid_body_inertia(rho=rho)
-
-
-def _set_property(
-    fb: FloatingBody,
-    property: str,
-    value: Optional[Iterable[float]],
-    rho: float = _default_parameters["rho"],
-) -> FloatingBody:
-    """Sets default properties if not provided by the user:
-        - `center_of_mass` is set to the geometric centroid
-        - `mass` is set to the displaced mass
-        - `rotation_center` is set to the center of mass
-    """
-    valid_properties = ['mass', 'center_of_mass', 'rotation_center']
-    if property not in valid_properties:
-        raise ValueError(
-            "`property` is not a recognized property. Valid properties are " +
-           f"{valid_properties}."
-        )
-    if not hasattr(fb, property):
-        setattr(fb, property, None)
-    prop_org = getattr(fb, property) is not None
-    prop_new = value is not None
-
-    if not prop_org and not prop_new:
-        if property == 'mass':
-            vol = fb.copy(
-                name=f"{fb.name}_immersed"
-                ).keep_immersed_part().volume
-            def_val = rho * vol
-            log_str = (
-                "Setting the mass to the displaced mass.")
-        elif property == 'center_of_mass':
-            def_val = fb.center_of_buoyancy
-            log_str = (
-                "Using the geometric centroid as the center of gravity (COG).")
-        elif property == 'rotation_center':
-            def_val = fb.center_of_mass
-            log_str = (
-                "Using the center of gravity (COG) as the rotation center " +
-                "for hydrostatics.")
-        setattr(fb, property, def_val)
-        _log.info(log_str)
-    elif prop_org and prop_new:
-        if not np.allclose(getattr(fb, property), value):
-            raise ValueError(
-               f"Both :python:`fb.{property}` and " +
-               f":python:`{property}` were provided but have " +
-                "different values."
-            )
-    elif prop_new:
-        setattr(fb, property, value)
-
-    return fb