diff --git a/setup.cfg b/setup.cfg
index aa1d709fbe2d258183ac9990729b8ab6bc509f98..4e37be177faadd4cdc12d8d9df833414ff7367b1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -14,6 +14,7 @@ classifiers = Development Status :: 4 - Beta,
               Programming Language :: Python :: 3,
               Programming Language :: Python :: 3.3,
               Programming Language :: Python :: 3.4,
+			  Programming Language :: Python :: 3.5,
               Environment :: Console,
               Intended Audience :: Education,
               Intended Audience :: Science/Research,
diff --git a/wetb/__init__.py b/wetb/__init__.py
index 24f6acf059a1ef5281ab652238ad343fc3f9a518..a4c36dbed94f3e105321a29086a19d6f89cd7a84 100644
--- a/wetb/__init__.py
+++ b/wetb/__init__.py
@@ -4,9 +4,9 @@ from __future__ import division
 from __future__ import absolute_import
 from future import standard_library
 standard_library.install_aliases()
-import pkg_resources
-test = "TEST" 
+test = "TEST"
 try:
-    __version__ = pkg_resources.get_distribution(__name__).version
+    import pkg_resources
+    __version__ = pkg_resources.safe_version(pkg_resources.get_distribution(__name__).version)
 except:
     __version__ = 'unknown'
diff --git a/wetb/hawc2/simulation.py b/wetb/hawc2/simulation.py
index 0c6f584a0a4e35cba74b332571be6ead43679c78..589fc44b997b11ffe4f5d9f7d7312736b8ecdfb0 100755
--- a/wetb/hawc2/simulation.py
+++ b/wetb/hawc2/simulation.py
@@ -85,7 +85,7 @@ class Simulation(object):
         if not os.path.isabs(htcfilename):
             htcfilename = os.path.join(modelpath, htcfilename)
         self.filename = os.path.basename(htcfilename)
-        self.htcFile = H2aeroHTCFile(htcfilename)
+        self.htcFile = HTCFile(htcfilename)
         self.time_stop = self.htcFile.simulation.time_stop[0]
         self.hawc2exe = hawc2exe
         self.copy_turbulence = copy_turbulence
diff --git a/wetb/hawc2/st_file.py b/wetb/hawc2/st_file.py
index 1b22366757a8bcdee47c97f4b1f6932b3e9ea7ba..f3c36d4470cc1434cc550e3ae441461438009efe 100644
--- a/wetb/hawc2/st_file.py
+++ b/wetb/hawc2/st_file.py
@@ -22,8 +22,8 @@ class StFile(object):
 
     - r : curved length distance from main_body node 1 [m]
     - m : mass per unit length [kg/m]
-    - xm : xc2-coordinate from C1/2 to mass center [m]
-    - ym : yc2-coordinate from C1/2 to mass center [m]
+    - x_cg : xc2-coordinate from C1/2 to mass center [m]
+    - y_cg : yc2-coordinate from C1/2 to mass center [m]
     - rix : radius of gyration related to elastic center. Corresponds to rotation about principal bending xe axis [m]
     - riy : radius of gyration related to elastic center. Corresponds to rotation about principal bending ye axis [m]
     - xs : xc2-coordinate from C1/2 to shear center [m]. The shear center is the point where external forces only contributes to pure bending and no torsion.
diff --git a/wetb/utils/caching.py b/wetb/utils/caching.py
index 6f5df273165e5ffd56b539017b239165b8265796..84a38a69b5a171244a2aadc4519cbacb3415646a 100644
--- a/wetb/utils/caching.py
+++ b/wetb/utils/caching.py
@@ -8,6 +8,7 @@ from __future__ import print_function
 from __future__ import division
 from __future__ import absolute_import
 from future import standard_library
+import sys
 standard_library.install_aliases()
 import inspect
 
@@ -80,7 +81,11 @@ def cache_function(f):
             self.cache_attr_lst.add(name)
 
         return getattr(self, name)
-    if 'reload' in inspect.getargspec(f)[0]:
+    version = sys.version_info
+    if version >= (3,3):
+        if 'reload' in inspect.signature(f).parameters.values():
+            raise AttributeError("Functions decorated with cache_function are not allowed to take a parameter called 'reload'")
+    elif 'reload' in inspect.getargspec(f)[0]:
         raise AttributeError("Functions decorated with cache_function are not allowed to take a parameter called 'reload'")
     return wrap
 
diff --git a/wetb/utils/cython_compile/cython_compile.py b/wetb/utils/cython_compile/cython_compile.py
index 2d09be147d2454350552da9f3ea63a3c055b3ba3..1e17e16f68208962ccf5a20e72efe528e29c178b 100644
--- a/wetb/utils/cython_compile/cython_compile.py
+++ b/wetb/utils/cython_compile/cython_compile.py
@@ -151,7 +151,7 @@ def cython_import(import_module_name, compiler=None):
         fid.close()
 
         # compile, import compiled module and delete temporary files
-        module_relname = os.path.relpath(eval(module_name).__file__, str(os.getcwd())).replace(os.path.sep, ".")[:-3]
+        module_relname = eval(module_name).__package__ + "." + module_name  #"os.path.relpath(eval(module_name).__file__, str(os.getcwd())).replace(os.path.sep, ".")[:-3]
         return compile_and_cleanup(import_module_name, pyx_filename, module_relname, compiler)
     return eval(module_name)
 
diff --git a/wetb/utils/geometry.py b/wetb/utils/geometry.py
index 05c7369e4b3a09ecf79e2ea86082933b2e700ae6..e82594ab3f77f7c8e43963006ec834553e223c84 100644
--- a/wetb/utils/geometry.py
+++ b/wetb/utils/geometry.py
@@ -162,7 +162,7 @@ def rpm2rads(rpm):
     return rpm * 2 * np.pi / 60
 
 
-def abvrel2xyz(alpha, beta, vrel):
+def abvrel2xyz_old(alpha, beta, vrel):
     """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
 
     Parameters
@@ -203,3 +203,57 @@ def abvrel2xyz(alpha, beta, vrel):
     z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 + 1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2))
 
     return x, y, z
+
+def abvrel2xyz(alpha, beta, vrel):
+    """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
+
+    x : parallel to pitot tube, direction pitot tube root to tip, i.e. normal flow gives negative x\n
+    y : component in alpha plane
+    z : component in beta plane
+
+    For typical usage where pitot tube is mounted on leading edge:\n
+    x: Opposite rotational direction\n
+    y: Direction of mean wind\n
+    z: From blade root to tip\n
+
+
+    Parameters
+    ----------
+    alpha : array_like
+        Pitot tube angle of attack [rad]. Zero for flow towards pitot tube. Positive around z-axis. I.e.
+        negative alpha (normal flow) gives positive y component
+    beta : array_like
+        Pitot tube side slip angle [rad]. Zero for flow towards pitot tube. Positive around y-axis. I.e.
+        Positive beta (normal flow due to expansion and position in front of blade) gives positive z
+    vrel : array_like
+        Pitot tube relative velocity. Positive: flow towards pitot tube
+
+    Returns
+    -------
+    x : array_like
+        Wind component away from pitot tube (positive for postive vrel and -90<beta<90)
+    y : array_like
+        Wind component in alpha plane (positive for positive alpha)
+    z : array_like
+        Wind component in beta plane (positive for negative beta)
+    """
+    alpha = np.array(alpha, dtype=np.float)
+    beta = np.array(beta, dtype=np.float)
+    vrel = np.array(vrel, dtype=np.float)
+
+    sign_vsx = ((np.abs(beta) > np.pi / 2) * 2 - 1)  # -1 for |beta| < 90, +1 for |beta|>90
+    sign_vsy = -np.sign(alpha)  #- for alpha > 0
+    sign_vsz = np.sign(beta)  # for beta>0
+
+
+    x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
+
+    m = alpha != 0
+    y = np.zeros_like(alpha)
+    y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 + 1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2))
+
+    m = beta != 0
+    z = np.zeros_like(alpha)
+    z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 + 1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2))
+
+    return np.array([x, y, z]).T
diff --git a/wetb/utils/tests/test_geometry.py b/wetb/utils/tests/test_geometry.py
index 291864c108a2819e279be9c787e92899305241f2..e2fbf1f1faa5d951d611b4386cf0a95ca800ab42 100644
--- a/wetb/utils/tests/test_geometry.py
+++ b/wetb/utils/tests/test_geometry.py
@@ -14,7 +14,7 @@ import unittest
 import wetb.gtsdf
 import numpy as np
 from wetb.utils.geometry import rad, deg, mean_deg, sind, cosd, std_deg, xyz2uvw, \
-    wsp_dir2uv, wsp_dir_tilt2uvw, tand
+    wsp_dir2uv, wsp_dir_tilt2uvw, tand, abvrel2xyz
 import os
 
 
@@ -117,6 +117,15 @@ class TestGeometry(unittest.TestCase):
         np.testing.assert_array_almost_equal(wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal=True, dir_ref=180), np.array([x, -y, z]), 5)
         np.testing.assert_array_almost_equal(xyz2uvw(*wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal=True), left_handed=False), xyz2uvw(x, y, z))
 
+    def test_abvrel2xyz(self):
+        abvrel = np.array([(0., 0, 1), (30, 0, 1), (-30, 0, 1), (0, 30, 1), (0, -30, 1), (30, 30, 1), (30, 30, 2)])
+        abvrel[:, :2] = rad(abvrel[:, :2])
+        for (x, y, z), (a, b, vrel) in zip(abvrel2xyz(*abvrel.T), abvrel):
+            #print (deg(a), deg(b), vrel, x, y, z)
+            self.assertAlmostEqual(a, np.arctan(y / x))
+            self.assertAlmostEqual(b, np.arctan(-z / x))
+            self.assertAlmostEqual(vrel, np.sqrt(x ** 2 + y ** 2 + z ** 2))
+
 
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.test_rad']
diff --git a/wetb/wind/shear.py b/wetb/wind/shear.py
index 4916a86ce087f06f6b0c892a36fb4b799496ad17..dd6122470e3a5a460660b7673371fb9042a147dd 100644
--- a/wetb/wind/shear.py
+++ b/wetb/wind/shear.py
@@ -19,7 +19,7 @@ def _z_u(z_u_lst):
     u = np.array([np.mean(np.array([u])[:]) for _, u in z_u_lst])
     return z, u
 
-def power_shear(alpha, z_ref, u_ref, z):
+def power_shear(alpha, z_ref, u_ref):
     """Power shear
 
     Parameters
@@ -28,22 +28,20 @@ def power_shear(alpha, z_ref, u_ref, z):
         The alpha shear parameter
     z_ref : int or float
         The reference height
-    z : int, float or array_like
-        The heights of interest
     u_ref : int or float
         The wind speed of the reference height
 
     Returns
     -------
-    power_shear : array-like
-        Wind speeds at height z
+    power_shear : function
+        Function returning for wsp at input heights: f(height) -> wsp
 
     Example
     --------
-    >>> power_shear(.5, 70, [20,50,70], 9)
+    >>> power_shear(.5, 70, 9)([20,50,70])
     [ 4.81070235  7.60638829  9.        ]
     """
-    return u_ref * (np.array(z) / z_ref) ** alpha
+    return lambda z : u_ref * (np.array(z) / z_ref) ** alpha
 
 
 def fit_power_shear(z_u_lst):
@@ -72,7 +70,7 @@ def fit_power_shear(z_u_lst):
     alpha, _ = np.polyfit(np.log(z / z_hub), np.log(u / u_hub), 1)
     return alpha
 
-def fit_power_shear_ref(z_u_lst, z_ref):
+def fit_power_shear_ref(z_u_lst, z_ref, plt=None):
     """Estimate power shear parameter, alpha, from two or more specific reference heights using polynomial fit.
 
     Parameters
@@ -84,6 +82,8 @@ def fit_power_shear_ref(z_u_lst, z_ref):
         - u_z2: Wind speeds or mean wind speeds at z2
     z_ref : float or int
         Reference height (hub height)
+    plt : matplotlib.pyplot (or similar) or None
+        Used to plot result if not None
 
     Returns
     -------
@@ -101,7 +101,13 @@ def fit_power_shear_ref(z_u_lst, z_ref):
         alpha, u_ref = x
         return np.sum([(u - u_ref * (z / z_ref) ** alpha) ** 2 for z, u in z_u_lst])
     z_u_lst = [(z, np.mean(u)) for z, u in z_u_lst]
-    return fmin(shear_error, (.1, 10), (z_u_lst, z_ref), disp=False)
+    alpha, u_ref = fmin(shear_error, (.1, 10), (z_u_lst, z_ref), disp=False)
+    if plt:
+        z, u = list(zip(*z_u_lst))
+        plt.plot(u, z, '.')
+        z = np.linspace(min(z), max(z), 100)
+        plt.plot(power_shear(alpha, z_ref, u_ref)(z), z)
+    return alpha, u_ref
 
 
 
@@ -188,7 +194,7 @@ if __name__ == '__main__':
     from matplotlib.pyplot import plot, show
     z = np.arange(0, 211)
     for alpha, c in zip([0.00001, 1, 2], ['r', 'b', 'g']):
-        u = power_shear(alpha, 120, 10, z)
+        u = power_shear(alpha, 120, 10)(z)
         plot(u, z, c)
         plot(u.mean(), 120, c + '.')
 
diff --git a/wetb/wind/tests/test_Shear.py b/wetb/wind/tests/test_Shear.py
index af1f6f0a539395b97053cddd7b11f94ce86aef71..318d1c32dcddfad62398f345ffea1a65743f0ff9 100644
--- a/wetb/wind/tests/test_Shear.py
+++ b/wetb/wind/tests/test_Shear.py
@@ -98,8 +98,10 @@ class TestShear(unittest.TestCase):
         wsp53 = data[:, 2]
         wsp21 = data[:, 4]
 
-
-        alpha, u_ref = fit_power_shear_ref([(85, wsp85), (21, wsp21)], 87.13333)
+        import matplotlib.pyplot as plt
+        alpha, u_ref = fit_power_shear_ref([(85, wsp85), (21, wsp21)], 87.13333, plt)
+        if 0:
+            plt.show()
         self.assertAlmostEqual(alpha, .5, delta=.001)
         self.assertAlmostEqual(u_ref, 9, delta=.01)