diff --git a/arlpy/uwa.py b/arlpy/uwa.py
index 6e310d0..14667bc 100644
--- a/arlpy/uwa.py
+++ b/arlpy/uwa.py
@@ -167,20 +167,25 @@ def doppler(speed, frequency, c=soundspeed()):
     """
     return (1+speed/float(c))*frequency
 
-def bubble_resonance(radius, depth=0):
-    """Get the resonant frequency of a bubble of a given radius.
-
-    The bubble resonance is computed based on Medwin & Clay (1998).
-
-    :param radius: radius of the bubble in m
-    :param depth: depth in m
-    :returns: resonant frequency of the bubble in Hz
+def bubble_resonance(radius, depth=0.0, gamma = 1.4, p0 = 1.013e5, rho_water = 1022.476):
+    """Compute resonance frequency of a freely oscillating has bubble in water,
+    using implementation based on Medwin & Clay (1998). This ignores surface-tension,
+    thermal, viscous and acoustic damping effects, and the pressure-volume relationship
+    is taken to be adiabatic. Parameters:
+
+    :radius: bubble `radius` in meters
+    :depth: depth of bubble in water in meters
+    :gamma: gas ratio of specific heats. Default 1.4 for air
+    :p0: atmospheric pressure. Default 1.013e5 Pa
+    :rho_water: Density of water. Default 1022.476 kg/m³
 
     >>> import arlpy
     >>> arlpy.uwa.bubble_resonance(100e-6)
-    32500.0
+    32465.56
     """
-    return 3.25/radius * _np.sqrt(1+0.1*depth)
+    g = 9.80665 #acceleration due to gravity
+    p_air = p0 + rho_water*g*depth
+    return 1/(2*_np.pi*radius)*_np.sqrt(3*gamma*p_air/rho_water)
 
 def bubble_surface_loss(windspeed, frequency, angle):
     """Get the surface loss due to bubbles.
@@ -226,21 +231,21 @@ def bubble_soundspeed(void_fraction, c=soundspeed(), c_gas=340, relative_density
 
 def pressure(x, sensitivity, gain, volt_params=None):
     """Convert the real signal x to an acoustic pressure signal in micropascal.
-    
+
     :param x: real signal in voltage or bit depth (number of bits)
     :param sensitivity: receiving sensitivity in dB re 1V per micropascal
     :param gain: preamplifier gain in dB
-    :param volt_params: (nbits, v_ref) is used to convert the number of bits 
-        to voltage where nbits is the number of bits of each sample and v_ref 
-        is the reference voltage, default to None  
+    :param volt_params: (nbits, v_ref) is used to convert the number of bits
+        to voltage where nbits is the number of bits of each sample and v_ref
+        is the reference voltage, default to None
     :returns: acoustic pressure signal in micropascal
-    
+
     If `volt_params` is provided, the sample unit of x is in number of bits,
-    else is in voltage.  
+    else is in voltage.
 
-    >>> import arlpy 
+    >>> import arlpy
     >>> nbits = 16
-    >>> V_ref = 1.0 
+    >>> V_ref = 1.0
     >>> x_volt = V_ref*signal.cw(64, 1, 512)
     >>> x_bit = x_volt*(2**(nbits-1))
     >>> sensitivity = 0
@@ -251,17 +256,17 @@ def pressure(x, sensitivity, gain, volt_params=None):
     nu = 10**(sensitivity/20)
     G = 10**(gain/20)
     if volt_params is not None:
-        nbits, v_ref = volt_params   
+        nbits, v_ref = volt_params
         x = x*v_ref/(2**(nbits-1))
     return x/(nu*G)
 
 def spl(x, ref=1):
     """Get Sound Pressure Level (SPL) of the acoustic pressure signal x.
-    
+
     :param x: acoustic pressure signal in micropascal
     :param ref: reference acoustic pressure in micropascal, default to 1
     :returns: average SPL in dB re micropascal
-    
+
     In water, the common reference is 1 micropascal. In air, the common
     reference is 20 micropascal.
 
diff --git a/arlpy/uwapm.py b/arlpy/uwapm.py
index fbc3f55..d1be189 100644
--- a/arlpy/uwapm.py
+++ b/arlpy/uwapm.py
@@ -108,7 +108,8 @@ def create_env2d(**kv):
         'depth_interp': linear,         # curvilinear/linear
         'min_angle': -80,               # deg
         'max_angle': 80,                # deg
-        'nbeams': 0                     # number of beams (0 = auto)
+        'nbeams': 0,                     # number of beams (0 = auto)
+        'stepsize': 0                     # raytrace step size (0 = auto)
     }
     for k, v in kv.items():
         if k not in env.keys():
@@ -695,7 +696,7 @@ def _create_env_file(self, env, taskcode):
             self._create_sbp_file(fname_base+'.sbp', env['tx_directionality'])
         self._print(fh, "%d" % (env['nbeams']))
         self._print(fh, "%0.6f %0.6f /" % (env['min_angle'], env['max_angle']))
-        self._print(fh, "0.0 %0.6f %0.6f" % (1.01*max_depth, 1.01*_np.max(env['rx_range'])/1000))
+        self._print(fh, "%0.6f %0.6f %0.6f" % (env['stepsize'], 1.01*max_depth, 1.01*_np.max(env['rx_range'])/1000))
         _os.close(fh)
         return fname_base
 
diff --git a/tests/test_basic.py b/tests/test_basic.py
index 9f7f997..9d7f739 100644
--- a/tests/test_basic.py
+++ b/tests/test_basic.py
@@ -130,9 +130,9 @@ def test_doppler(self):
         self.assertApproxEqual(uwa.doppler(-10, 50000), 49675)
 
     def test_bubble_resonance(self):
-        self.assertApproxEqual(uwa.bubble_resonance(100e-6), 32500)
-        self.assertApproxEqual(uwa.bubble_resonance(32e-6), 101562)
-        self.assertApproxEqual(uwa.bubble_resonance(100e-6, depth=10), 45962)
+        self.assertApproxEqual(uwa.bubble_resonance(100e-6), 32465.562964469198)
+        self.assertApproxEqual(uwa.bubble_resonance(32e-6), 101454.88426396625)
+        self.assertApproxEqual(uwa.bubble_resonance(100e-6, depth=10), 45796.45437634176)
 
     def test_bubble_surface_loss(self):
         self.assertApproxEqual(utils.mag2db(uwa.bubble_surface_loss(3, 10000, 0)), -1.44, precision=2)
@@ -277,7 +277,7 @@ def test_detect_impulses(self):
         x += np.random.normal(0, 0.1, nsamp)
         ind_imp, _ = signal.detect_impulses(x, fs=100000, k=10, tdist=1e-3)
         self.assertArrayEqual(true_ind_imp, ind_imp)
-        
+
 class CommsTestSuite(MyTestCase):
 
     def test_random_data(self):