diff --git a/Project.toml b/Project.toml
index 403aa08..91ab9aa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,4 +1,4 @@
-name = "BEM"
+name = "MarineHydro"
 uuid = "8d7a907f-d253-4c14-a318-41dd3f310a02"
 authors = ["KapilKhanal <kk733@cornell.edu>"]
 version = "0.1.0"
diff --git a/README.md b/README.md
index 14a8495..028447e 100644
--- a/README.md
+++ b/README.md
@@ -42,7 +42,7 @@ Fully-differentiable boundary element solver for marine hydrodynamics. This new
   Includes plots and data generated for the paper.
 
 - **📜 src**  
-  Source code files for the `BEM.jl` package, including the main functionality.
+  Source code files for the `MarineHydro.jl` package, including the main functionality.
 
 - **🧪 test**  
   Contains test files and resources to verify the functionality of the source code.
@@ -74,10 +74,10 @@ Fully-differentiable boundary element solver for marine hydrodynamics. This new
    ENV["PYTHON"] = "/path/to/capytaine_env/bin/python" # run 'which python' in your terminal for this
    ```
   
-6. **Using the BEM Module**  
+6. **Using the MarineHydro Module**  
 Load the module in your Julia session:  
 ```julia
-using BEM
+using MarineHydro
 using PyCall
 # import your capytaine mesh
 cpt = pyimport("capytaine")
diff --git a/paper/ArraySensitivities.jl b/paper/ArraySensitivities.jl
index 9ae33fc..fcd6d35 100644
--- a/paper/ArraySensitivities.jl
+++ b/paper/ArraySensitivities.jl
@@ -1,7 +1,7 @@
 using Plots
 using CSV
 using DataFrames
-using BEM
+using MarineHydro
 include("/home/cornell/BEMJulia/BEM.jl/paper/meshGradients_pair.jl")
  #takes a while depending on this - some faces_max_radius gives weird answer
 #check meshes.jl to change faces_max_radius.
@@ -23,7 +23,7 @@ function added_mass_off_diagonal(radius,omega ,dx1)
     sphere_2_heave_normal = [element_is_in_sphere_2(j) ? mesh.normals[j,:]' * [0.0, 0.0, 1.0] : 0.0 for j in 1:total_nfaces]
     k = omega^2 / 9.81  # Wave number
     S, D = assemble_matrices((Rankine(), RankineReflected(), GFWu()), mesh, k)
-    potential = BEM.solve(D, S, -1im * omega * sphere_1_heave_normal)
+    potential = MarineHydro.solve(D, S, -1im * omega * sphere_1_heave_normal)
     pressure = 1im * 1000 * omega * potential
     # force = -sum(pressure .* sphere_2_heave_normal .* mesh.areas)
     # A12 = real(force) / omega^2
@@ -41,7 +41,7 @@ function damping_off_diagonal(radius,omega ,dx1)
     sphere_2_heave_normal = [element_is_in_sphere_2(j) ? mesh.normals[j,:]' * [0.0, 0.0, 1.0] : 0.0 for j in 1:total_nfaces]
     k = omega^2 / 9.81  # Wave number
     S, D = assemble_matrices((Rankine(), RankineReflected(), GFWu()), mesh, k) # Assemble matrices tuple error -- use default
-    potential = BEM.solve(D, S, -1im * omega * sphere_1_heave_normal)
+    potential = MarineHydro.solve(D, S, -1im * omega * sphere_1_heave_normal)
     pressure = 1im * 1000 * omega * potential
     force = -sum(pressure .* sphere_1_heave_normal .* mesh.areas)
     B11 = imag(force) / omega
diff --git a/paper/CoeficientsGradients.jl b/paper/CoeficientsGradients.jl
index 7332720..65a96db 100644
--- a/paper/CoeficientsGradients.jl
+++ b/paper/CoeficientsGradients.jl
@@ -1,4 +1,4 @@
-using BEM
+using MarineHydro
 using Zygote
 using FiniteDifferences
 using Test
diff --git a/paper/ComparisonPlots.jl b/paper/ComparisonPlots.jl
index d7ac5d8..e134463 100644
--- a/paper/ComparisonPlots.jl
+++ b/paper/ComparisonPlots.jl
@@ -1,6 +1,6 @@
 ## COEFFICIENTS CHECKS - DO NOT DELETE even if commented out; MAKE TEST
 using LinearAlgebra,Statistics,Plots
-using BEM
+using MarineHydro
 using ColorTypes
 include("MeshGradients_singlebody.jl")
 include("HulmeData.jl")
@@ -184,4 +184,4 @@ _grad_fine =  [6552.7934287808785, 26211.17371512438, 58975.14085905893, 104844.
 #plot(radii, _grad_coarse, xlabel="r [m]", ylabel="∂A_∂ω", label="AD (coarse mesh)",marker = "*")
 plot(radii, _gradients , xlabel="r [m]", ylabel="∂A_∂r", label="AD",marker = "*",color = vermillion)
 plot!(radii,analy_grady  , xlabel="r [m]", label="analytical",marker = "*" , color = bluishgreen)
-savefig("/home/cornell/BEMJulia/MarineHydro.jl/paper/Plots/analy_ad_AMass_surge.pdf")
\ No newline at end of file
+savefig("/home/cornell/BEMJulia/MarineHydro.jl/paper/Plots/analy_ad_AMass_surge.pdf")
diff --git a/paper/DelhommeauVsWu.jl b/paper/DelhommeauVsWu.jl
index 8cd6c01..550fd51 100644
--- a/paper/DelhommeauVsWu.jl
+++ b/paper/DelhommeauVsWu.jl
@@ -1,4 +1,4 @@
-using PyCall, BEM, Zygote, Plots, ColorTypes
+using PyCall, MarineHydro, Zygote, Plots, ColorTypes
 
 orange = RGB(230/255,159/255,0/255)  
 vermillion = RGB(213/255, 94/255, 0/255) 
@@ -63,7 +63,7 @@ xlabel!("k")
 ylabel!("Gradient")
 title!("Comparison of Gradients: Wu vs Delhommeau")
 legend()
-savefig("/home/cornell/BEMJulia/BEM.jl/paper/Plots/added_mass_del_wu.pdf")
+savefig("/home/cornell/BEMJulia/MarineHydro.jl/paper/Plots/added_mass_del_wu.pdf")
 
 plot(k_range, A_w_gradwu, label="A_w_gradwu", lw=2, linestyle=:dash, color=vermillion)
 plot!(k_range, A_w_graddel, label="A_w_graddel", lw=2, linestyle=:solid, color=bluishgreen)
@@ -71,4 +71,4 @@ xlabel!("k")
 ylabel!("Gradient")
 title!("Comparison of Gradients: Wu vs Delhommeau")
 legend()
-savefig("/home/cornell/BEMJulia/BEM.jl/paper/Plots/gradw_added_mass_del_wu.pdf")
\ No newline at end of file
+savefig("/home/cornell/BEMJulia/MarineHydro.jl/paper/Plots/gradw_added_mass_del_wu.pdf")
diff --git a/paper/HydrostaticsRule.jl b/paper/HydrostaticsRule.jl
index cf23691..66582e6 100644
--- a/paper/HydrostaticsRule.jl
+++ b/paper/HydrostaticsRule.jl
@@ -1,6 +1,6 @@
 using Zygote
 using ChainRulesCore
-using BEM 
+using MarineHydro 
 cpt = pyimport("capytaine")
 
 ∂J_r_fd_hs(f,r; h=1e-5) = (f(r+h) .- f(r-h)) ./ (2*h)
diff --git a/paper/Multibody.jl b/paper/Multibody.jl
index eedde4e..a49fdcf 100644
--- a/paper/Multibody.jl
+++ b/paper/Multibody.jl
@@ -1,4 +1,4 @@
-using BEM
+using MarineHydro
 using PyCall
 using LinearAlgebra
 
@@ -33,7 +33,7 @@ Am = zeros((2, 2))
 Bm = zeros((2, 2))
 
 # radiation of first sphere
-potential = BEM.solve(D, S, -1im * omega * sphere_1_heave_normal)
+potential = MarineHydro.solve(D, S, -1im * omega * sphere_1_heave_normal)
 pressure = 1im * 1000 * omega * potential 
 force_on_sphere_1 = -sum(pressure .* sphere_1_heave_normal .* mesh.areas)
 A11 = real(force_on_sphere_1)/omega^2
@@ -44,7 +44,7 @@ B12 = imag(force_on_sphere_2)/omega
 Am = [A11 A12; A12 A11]
 Bm = [B11 B12; B12 B11]
 # # radiation of second sphere (very symmetric here, but that is just for demonstration)
-# potential = BEM.solve(D, S, -1im * omega * sphere_2_heave_normal)
+# potential = MarineHydro.solve(D, S, -1im * omega * sphere_2_heave_normal)
 # pressure = 1im * 1000 * omega * potential 
 # force_on_sphere_1 = -sum(pressure .* sphere_1_heave_normal .* mesh.areas)
 # Am[2, 1] = real(force_on_sphere_1)/omega^2
diff --git a/paper/Power.jl b/paper/Power.jl
index b2be2c6..2c438eb 100644
--- a/paper/Power.jl
+++ b/paper/Power.jl
@@ -4,7 +4,7 @@ using ForwardDiff
 using Zygote, ChainRulesCore 
 ForwardDiff.can_dual(::Type{Any}) = true
 using LinearAlgebra
-using BEM
+using MarineHydro
 using ImplicitAD
 include("/home/cornell/BEMJulia/MarineHydro.jl/paper/meshGradients_pair.jl")
 include("/home/cornell/BEMJulia/MarineHydro.jl/paper/HydrostaticsRule.jl")
@@ -38,7 +38,7 @@ function power(r1,dx1)
     # sphere_2_center = hcat(sphere_2_center'...)
 
     # radiation of first sphere
-    potential = BEM.solve(D, S, -1im * omega * sphere_1_heave_normal)
+    potential = MarineHydro.solve(D, S, -1im * omega * sphere_1_heave_normal)
     pressure = 1im * 1000 * omega * potential 
     force_on_sphere_1 = -sum(pressure .* sphere_1_heave_normal .* mesh.areas)
     A11 = real(force_on_sphere_1)/omega^2
diff --git a/paper/direct_vs_indirect.jl b/paper/direct_vs_indirect.jl
index 43c9da6..7eb3b36 100644
--- a/paper/direct_vs_indirect.jl
+++ b/paper/direct_vs_indirect.jl
@@ -1,4 +1,4 @@
-using BEM
+using MarineHydro
 using Zygote
 using FiniteDifferences
 using Test
diff --git a/paper/meshGradients_pair.jl b/paper/meshGradients_pair.jl
index 027327c..ddaf132 100644
--- a/paper/meshGradients_pair.jl
+++ b/paper/meshGradients_pair.jl
@@ -209,7 +209,7 @@ function ChainRulesCore.rrule(::typeof(sphere_faces), r,dx1)
     return Mesh(vertices,faces,centers,normals,areas,radii,nvertices,nfaces) #
 end
 
-using BEM
+using MarineHydro
   #test differentiability
 function func(r1,dx1)
     mesh = differentiableMeshPairs(r1,dx1)
diff --git a/paper/profileBEM.jl b/paper/profileBEM.jl
index a6ef484..3b5f776 100644
--- a/paper/profileBEM.jl
+++ b/paper/profileBEM.jl
@@ -1,4 +1,4 @@
-using BEM
+using MarineHydro
 using BenchmarkTools
 include("BemProgram.jl")
 include("MeshGradients_singlebody.jl")
@@ -24,21 +24,21 @@ element_2 = (
 
 function integral_and_gradient(GF, element_1, element_2, wavenumber=nothing)
     return (
-        BEM.integral(GF, element_1, element_2, wavenumber),
-        BEM.integral_gradient(GF, element_1, element_2, wavenumber)
+        MarineHydro.integral(GF, element_1, element_2, wavenumber),
+        MarineHydro.integral_gradient(GF, element_1, element_2, wavenumber)
     )
 end
 omega = 1.03
 
 println("Rankine")
-@btime BEM.integral($(Rankine()), $element_1, $element_2)
-@btime BEM.integral_gradient($(Rankine()), $element_1, $element_2)
+@btime MarineHydro.integral($(Rankine()), $element_1, $element_2)
+@btime MarineHydro.integral_gradient($(Rankine()), $element_1, $element_2)
 @btime integral_and_gradient($(Rankine()), $element_1, $element_2)
 
 wavenumber = 1.0
 println("GFWu")
-@btime BEM.integral($(GFWu()), $element_1, $element_2, $wavenumber)
-@btime BEM.integral_gradient($(GFWu()), $element_1, $element_2, $wavenumber)
+@btime MarineHydro.integral($(GFWu()), $element_1, $element_2, $wavenumber)
+@btime MarineHydro.integral_gradient($(GFWu()), $element_1, $element_2, $wavenumber)
 @btime integral_and_gradient($(GFWu()), $element_1, $element_2, $wavenumber)
 
 
@@ -52,4 +52,4 @@ B(radius) = damping_program(radius,omega,dof)
 @show A(1.0)
 @btime Zygote.gradient(A, 1.0)
 
-nothing
\ No newline at end of file
+nothing
diff --git a/src/BEM.jl b/src/MarineHydro.jl
similarity index 99%
rename from src/BEM.jl
rename to src/MarineHydro.jl
index 62f4daf..0d4bbd2 100644
--- a/src/BEM.jl
+++ b/src/MarineHydro.jl
@@ -1,5 +1,5 @@
 
-module BEM
+module MarineHydro
 
 export greens, gradient_greens, integral, integral_gradient, solve,calculate_radiation_forces
 export Mesh, element, integrate_pressure, DiffractionForce,assemble_matrix_wu
diff --git a/test/analytical.jl b/test/analytical.jl
index e94384f..9a001c0 100644
--- a/test/analytical.jl
+++ b/test/analytical.jl
@@ -1,4 +1,4 @@
-using BEM
+using MarineHydro
 using Test
 using PyCall
 
diff --git a/test/benchmarks_test.jl b/test/benchmarks_test.jl
index 43e0656..f41e38f 100644
--- a/test/benchmarks_test.jl
+++ b/test/benchmarks_test.jl
@@ -1,5 +1,5 @@
 using BenchmarkTools
-using BEM 
+using MarineHydro 
 using PyCall
 using Zygote
 
@@ -8,7 +8,7 @@ radius = 1.0
 resolution = (6, 6) #smaller data just to test
 cptmesh = cpt.mesh_sphere(name="sphere", radius=radius, center=(0, 0, 0), resolution=resolution)
 cptmesh.keep_immersed_part(inplace=true)
-mesh = BEM.Mesh(cptmesh)
+mesh = MarineHydro.Mesh(cptmesh)
 npanels = mesh.nfaces
 println( "Mesh information $resolution with  panels $npanels")
 omega = 1.12
diff --git a/test/coeff_gradient_tests.jl b/test/coeff_gradient_tests.jl
index bcb8d20..8f87d2b 100644
--- a/test/coeff_gradient_tests.jl
+++ b/test/coeff_gradient_tests.jl
@@ -1,6 +1,6 @@
 using Test
 using Zygote
-using BEM
+using MarineHydro
 using PyCall
 using LinearAlgebra
 using ImplicitAD
diff --git a/test/differentiability_tests.jl b/test/differentiability_tests.jl
index 889f2ed..71da5bb 100644
--- a/test/differentiability_tests.jl
+++ b/test/differentiability_tests.jl
@@ -1,6 +1,6 @@
 using Test
 using Zygote
-using BEM
+using MarineHydro
 using PyCall
 using LinearAlgebra
 using ImplicitAD
diff --git a/test/green_function.jl b/test/green_function.jl
index a0f5bf9..d827873 100644
--- a/test/green_function.jl
+++ b/test/green_function.jl
@@ -1,5 +1,5 @@
 using Test
-using BEM
+using MarineHydro
 
 @testset "Green function" begin
 
@@ -31,7 +31,7 @@ using BEM
         end
 
         @testset "Alternative definitions of the reflected Rankine" begin
-            import BEM.free_surface_symmetry
+            import MarineHydro.free_surface_symmetry
 
             @test greens(Rankine(), free_surface_symmetry(e1), e2) ≈ greens(Rankine(), e1, free_surface_symmetry(e2))
             @test gradient_greens(Rankine(), free_surface_symmetry(e1), e2, with_respect_to_first_variable=true) ≈ vertical_reflection(gradient_greens(Rankine(), e1, free_surface_symmetry(e2), with_respect_to_first_variable=true))
@@ -54,7 +54,7 @@ using BEM
 
     @testset "Guével-Delhommeau" begin
 
-        import BEM._dimless_wave_term, BEM._d_dimless_wave_term_dr
+        import MarineHydro._dimless_wave_term, MarineHydro._d_dimless_wave_term_dr
         # From cpt.Delhommeau().fortran_core.delhommeau_integrals.numerical_integration(1.0, -1.0, 1000)
         capy_ref = [-1.84003542, 1.76871985, -0.59117747, -1.01715699]
         @test _dimless_wave_term(1.0, -1.0) ≈ capy_ref[1] + im*capy_ref[2] atol=1e-4
diff --git a/test/with_capytaine.jl b/test/with_capytaine.jl
index 9cd6759..bf6833d 100644
--- a/test/with_capytaine.jl
+++ b/test/with_capytaine.jl
@@ -1,4 +1,4 @@
-using BEM
+using MarineHydro
 using Test
 using PyCall