diff --git a/examples/scripts/plot_phase_tensor_map.py b/examples/scripts/plot_phase_tensor_map.py
index 4617a740b..7e7b03ad9 100644
--- a/examples/scripts/plot_phase_tensor_map.py
+++ b/examples/scripts/plot_phase_tensor_map.py
@@ -15,11 +15,12 @@
 import mtpy.imaging.phase_tensor_maps as pptmaps
 
 # directory containing edis
+
 #edipath = r'C:\mtpywin\mtpy\examples\data\edi2'
 edipath = '/mtpy/examples/data/ET_edi'
 
 # whether or not to save the figure to file
-save = True
+save = False
 
 # full path to file to save to
 savepath = '/tmp'
@@ -28,8 +29,9 @@
 plot_freq = 1e-2
 
 # value to color ellipses by, options are phimin,phimax,skew
-colorby = 'skew'
+colorby = 'strike'
 ellipse_range = [-9, 9]
+rotation_angle = 0
 
 image_fn = 'phase_tensor_map%1is_' % (int(1. / plot_freq)) + colorby + '.png'
 
@@ -61,7 +63,7 @@
                                 background_image='/mtpy/examples/data/gravity/ET_gravity.tif'
  #                               bimg_band=1, #  Optional, set to None by default
                                 #bimg_cmap='viridis'  # Optional, set to 'viridis' by default
-                                )
+                                rotation_angle=rotation_angle)
 
 
 if save:
diff --git a/examples/scripts/plot_phase_tensor_section.py b/examples/scripts/plot_phase_tensor_section.py
index 979632cb7..3efc9025a 100644
--- a/examples/scripts/plot_phase_tensor_section.py
+++ b/examples/scripts/plot_phase_tensor_section.py
@@ -28,6 +28,7 @@
                                  plot_tipper = 'yri', # plot tipper ('y') + 'ri' means real+imag
                                  font_size=5,
                                  lw=0.5,
+                                 rotation_angle=20,
                                  ellipse_dict = {'ellipse_colorby':'skew_seg',# option to colour by phimin, phimax, skew, skew_seg
                                                  'ellipse_range':[-12, 12, 3]} # set color limits - default 0,90 for phimin or max,
                                                                          # [-12,12] for skew. If plotting skew_seg need to provide
@@ -37,8 +38,7 @@
 # update ellipse size (tweak for your dataset)
 plotObj.ellipse_size = 2.5
 
-
 plotObj.plot()
 
 #plotObj.save_figure(save_fn = op.join(savepath,'PhaseTensorSection.png'),
-#                    fig_dpi=400) # change to your preferred file resolution
\ No newline at end of file
+#                    fig_dpi=400) # change to your preferred file resolution
diff --git a/examples/scripts/plot_strike_roseplot.py b/examples/scripts/plot_strike_roseplot.py
index 4316fd230..574b6a6b5 100644
--- a/examples/scripts/plot_strike_roseplot.py
+++ b/examples/scripts/plot_strike_roseplot.py
@@ -11,12 +11,11 @@
 import os
 os.chdir(r'C:\mtpywin\mtpy') # change to path where mtpy is installed
 
-
 from mtpy.imaging.plotstrike import PlotStrike
 
 
 # directory containing edis
-edipath = r'C:\mtpywin\mtpy\examples\data\edi_files_2'
+edipath = r'C:\mtpywin\mtpy\data\edifiles2'
 
 
 # full path to file to save to
@@ -27,13 +26,50 @@
 elst = [op.join(edipath,f) for f in os.listdir(edipath) if (f.endswith('.edi'))]# and f.startswith('GL')
 
 
-strikeplot = PlotStrike(fn_list=elst,
-                        fold=False,
-                        show_ptphimin=False,
-                        plot_type=2 # 1 means divide into separate plots for different decades
-                                    # 2 means combine all data into one rose plot
-                                    )
-#strikeplot.save_plot(savepath,
-#                     file_format='png',
-#                     fig_dpi=400 
-#                     )
\ No newline at end of file
+### this will plot the estimated strike duplicated across quadrants
+#  Plot type 2 will plot all estimates of  strike into one ploe 
+strike_plot = PlotStrike(fn_list=elst,
+                         fold=False,
+                         plot_type=2)
+
+# If you want to plot the orthogonal estimation
+strike_plot.plot_orthogonal = True
+strike_plot.fig_num = 2
+strike_plot.plot()
+
+# if you want to plot Tipper strike estimates as well
+strike_plot.plot_orthogonal = False
+strike_plot.plot_tipper = True
+strike_plot.fig_num = 3
+strike_plot.plot()
+
+# if you want to rotate the data
+# note this will rotate the data N30W
+strike_plot.rotation_angle = -30
+strike_plot.fig_num = 4
+strike_plot.plot()
+
+# to rotate back
+strike_plot.rotation_angle = 30
+
+# if you want to plot per decade plots
+strike_plot.plot_type = 1
+strike_plot.fig_num = 5
+strike_plot.text_pad = 1.4
+strike_plot.plot()
+
+# if you want to only look at a few period ranges
+# not the range is given in log10 of the period
+strike_plot.plot_range = [-2, 0] 
+strike_plot.fig_num = 6
+strike_plot.plot()
+
+# if you want a vertical orientation instead of horizontal
+strike_plot.plot_orientation = 'v'
+strike_plot.plot_range = 'data'
+strike_plot.fig_num = 7
+strike_plot.plot()
+
+
+
+
diff --git a/mtpy/analysis/pt.py b/mtpy/analysis/pt.py
index c2ccf9c71..61a002836 100644
--- a/mtpy/analysis/pt.py
+++ b/mtpy/analysis/pt.py
@@ -108,7 +108,12 @@ def __init__(self, pt_array=None, pt_err_array=None, z_array=None,
     #  define get/set functions and properties
     # ==========================================================================
     # ---phase tensor array----------------------------------------
-    def _set_pt(self, pt_array):
+    @property
+    def pt(self):
+        return self._pt
+   
+    @pt.setter
+    def pt(self, pt_array):
         """
             Set the attribute 'pt'.
 
@@ -176,13 +181,13 @@ def _set_pt(self, pt_array):
         else:
             pass
 
-    def _get_pt(self):
+    # ---phase tensor Error-----------------------------------------------------
+    @property
+    def pt_err(self):
         return self._pt
 
-    pt = property(_get_pt, _set_pt, doc="Phase tensor array")
-
-    # ---phase tensor Error-----------------------------------------------------
-    def _set_pt_err(self, pt_err_array):
+    @pt_err.setter
+    def pt_err(self, pt_err_array):
         """
             Set the attribute 'pt_err'.
 
@@ -224,14 +229,13 @@ def _set_pt_err(self, pt_err_array):
         else:
             pass
 
-    def _get_pt_err(self):
-        return self._pt_err
-
-    pt_err = property(_get_pt_err, _set_pt_err,
-                      doc='Phase tensor error array, must be same shape as pt')
-
     # ---freq------------------------------------------------------------
-    def _set_freq(self, lo_freq):
+    @property
+    def freq(self):
+        return self._freq
+    
+    @freq.setter
+    def freq(self, lo_freq):
         """
             Set array of freq.
 
@@ -253,11 +257,6 @@ def _set_freq(self, lo_freq):
         except:
             self._freq = None
 
-    def _get_freq(self):
-        return self._freq
-
-    freq = property(_get_freq, _set_freq, doc="freq array")
-
     # ---z_object---------------------------------------------------------------
 
     def set_z_object(self, z_object):
@@ -278,12 +277,7 @@ def set_z_object(self, z_object):
                     self._pt[idx_f], self._pt_err[idx_f] = z2pt(self._z[idx_f],
                                                              self._z_err[idx_f])
                 except MTex.MTpyError_PT:
-                    try:
-                        print('Singular Matrix at {0:.5g} Hz'.format(
-                            self._freq[idx_f]))
-                    except AttributeError:
-                        print('Computed singular matrix')
-                        print('  --> pt[{0}]=np.zeros((2,2))'.format(idx_f))
+                    pass
 
         # --> if there is not error to the impedance tensor
         else:
@@ -291,30 +285,19 @@ def set_z_object(self, z_object):
                 try:
                     self._pt[idx_f] = z2pt(self._z[idx_f])[0]
                 except MTex.MTpyError_PT:
-                    try:
-                        print('Singular Matrix at {0:.5g}'.format(
-                            self._freq[idx_f]))
-                    except AttributeError:
-                        print('Computed singular matrix')
-                        print('  --> pt[{0}]=np.zeros((2,2))'.format(idx_f))
+                    pass
 
         self.rotation_angle = z_object.rotation_angle
 
-    # def _get_z_object(self):
-    #     z_object = MTz.Z(z_array=self._z, z_err_array=self._z_err)
-    #     z_object.freq = self._freq
-    #     z_object.rotation_angle = self.rotation_angle
-
-    #     return z_object
-
-    # _z_object = property(_get_z_object, _set_z_object, 
-    #                     doc="class mtpy.core.z.Z")
-
-
     # ---z array---------------------------------------------------------------
-    def _set_z(self, z_array):
+    @property
+    def z(self):
+        return self._z
+    
+    @z.setter
+    def z(self, z_array):
         """
-            Set  Z array as PhaseTensor object attribute.
+        Set  Z array as PhaseTensor object attribute.
         """
 
         self._z = z_array
@@ -327,12 +310,7 @@ def _set_z(self, z_array):
                     self._pt[idx_f], self._pt_err[idx_f] = z2pt(self._z[idx_f],
                                                              self._z_err[idx_f])
                 except MTex.MTpyError_PT:
-                    try:
-                        print('Singular Matrix at {0:.5g} Hz'.format(
-                            self._freq[idx_f]))
-                    except AttributeError:
-                        print('Computed singular matrix')
-                        print('  --> pt[{0}]=np.zeros((2,2))'.format(idx_f))
+                    pass
 
         # --> if there is not error to the impedance tensor
         elif self._z is not None:
@@ -340,20 +318,15 @@ def _set_z(self, z_array):
                 try:
                     self._pt[idx_f] = z2pt(self._z[idx_f])[0]
                 except MTex.MTpyError_PT:
-                    try:
-                        print('Singular Matrix at {0:.5g}'.format(
-                            self._freq[idx_f]))
-                    except AttributeError:
-                        print('Computed singular matrix')
-                        print('  --> pt[{0}]=np.zeros((2,2))'.format(idx_f))
-
-    # def _get_z(self):
-    #     return self._z
+                    pass
 
-    # z = property(_get_z, _set_z, 
-    #              doc="impedance tensor numpy.array((nf, 2, 2))")
 
     # ---Z Error array---------------------------------------------------------------
+    @property
+    def z_err(self):
+        return self._z_err
+    
+    @z_err.setter
     def _set_z_err(self, z_err_array):
         """
             Set  Z-error array as PhaseTensor object attribute.
@@ -373,12 +346,7 @@ def _set_z_err(self, z_err_array):
                     self.pt[idx_f], self.pt_err[idx_f] = z2pt(self._z[idx_f],
                                                              self._z_err[idx_f])
                 except MTex.MTpyError_PT:
-                    try:
-                        print('Singular Matrix at {0:.5g} Hz'.format(
-                            self._freq[idx_f]))
-                    except AttributeError:
-                        print('Computed singular matrix')
-                        print('  --> pt[{0}]=np.zeros((2,2))'.format(idx_f))
+                    pass
 
             # --> if there is not error to the impedance tensor
             else:
@@ -386,20 +354,7 @@ def _set_z_err(self, z_err_array):
                     try:
                         self._pt[idx_f] = z2pt(self._z[idx_f])[0]
                     except MTex.MTpyError_PT:
-                        try:
-                            print('Singular Matrix at {0:.5g}'.format(
-                                self._freq[idx_f]))
-                        except AttributeError:
-                            print('Computed singular matrix')
-                            print('  --> pt[{0}]=np.zeros((2,2))'.format(idx_f))
-
-    # def _get_z_err(self):
-    #     return self._z_err
-
-    # z_err = property(_get_z_err, _set_z_err, 
-    #                  doc="impedance tensor numpy.array((nf, 2, 2))")
-
-
+                        pass
 
     # ==========================================================================
     #  define get methods for read only properties
@@ -570,7 +525,7 @@ def azimuth(self):
     @property
     def azimuth_err(self):
         if self.pt_err is not None:
-            az_err = np.sqrt(self.alpha+self.beta)
+            az_err = np.sqrt(self.alpha + self.beta)
         else:
             az_err = None
 
@@ -740,14 +695,13 @@ def phimax(self):
         if self.pt is None:
             return None
 
-#        return self._pi2()[0] + self._pi1()[0]
         return np.degrees(np.arctan(self._pi2()[0] + self._pi1()[0]))
 
     @property
     def phimax_err(self):
         phimaxerr = None
         if self.pt_err is not None:
-            phimaxerr = np.sqrt(self._pi2()[1]**2+self._pi1()[1]**2)
+            phimaxerr = np.sqrt(self._pi2()[1]**2 + self._pi1()[1]**2)
  
             return np.degrees(np.arctan(phimaxerr))
         else:
@@ -755,14 +709,12 @@ def phimax_err(self):
 
     def rotate(self, alpha):
         """
-            Rotate PT array. Change the rotation angles attribute respectively.
-
-            Rotation angle must be given in degrees. All angles are referenced to 
-			geographic North, positive in clockwise direction. 
-			(Mathematically negative!)
+        Rotate PT array. Change the rotation angles attribute respectively.
 
-            In non-rotated state, X refs to North and Y to East direction.
+        Rotation angle must be given in degrees. All angles are referenced to 
+        North, positive in clockwise direction. (Mathematically negative!)
 
+        In non-rotated state, X refs to North and Y to East direction.
 
         """
 
@@ -819,9 +771,14 @@ def rotate(self, alpha):
                 angle = 0.
 
             if self.pt_err is not None:
-                pt_rot[idx_freq], pt_err_rot[idx_freq] = MTcc.rotatematrix_incl_errors(self.pt[idx_freq,:,:], angle, self.pt_err[idx_freq,:,:])
+                pt_rot[idx_freq], pt_err_rot[idx_freq] = \
+                    MTcc.rotate_matrix_with_errors(self.pt[idx_freq,:,:],
+                                                  -angle,
+                                                  self.pt_err[idx_freq,:,:])
             else:
-                pt_rot[idx_freq], pt_err_rot = MTcc.rotatematrix_incl_errors(self.pt[idx_freq,:,:], angle)
+                pt_rot[idx_freq], pt_err_rot = \
+                    MTcc.rotate_matrix_with_errors(self.pt[idx_freq,:,:],
+                                                  -angle)
 
         # --> set the rotated tensors as the current attributes
         self._pt = pt_rot
@@ -1306,87 +1263,4 @@ def z2pt(z_array, z_err_array=None):
                                         np.abs( -pt_array[idx_f,1,1] * realz[0,0] * z_err_array[1,1]) + \
                                         np.abs(  realz[0,0] * z_err_array[1,1]) + np.abs( -realz[1,0] * z_err_array[0,1]) )
 
-    return pt_array, pt_err_array
-
-
-def z_object2pt(z_object):
-    """
-        Calculate Phase Tensor from Z object (incl. uncertainties)
-
-        Input:
-        - Z-object : instance of the MTpy Z class
-
-
-        Return:
-        - PT object
-    """
-    #     - PT : nx2x2 real valued Numpy array
-    #     - PT-error : nx2x2 real valued Numpy array
-
-    # """
-
-    try:
-        p = PhaseTensor(z_object=z_object)
-    except:
-        raise MTex.MTpyError_Z('Input argument is not a valid instance of the Z class')
-
-    # pt_array = p.pt
-    # pterr_array = p.pterr
-
-    # return pt_array, pterr_array
-    return p
-
-
-def _edi_object2pt(edi_object):
-    """
-        Calculate Phase Tensor from Edi object (incl. uncertainties)
-
-        Input:
-        - Edi-object : instance of the MTpy Edi class
-
-        Return:
-        - PT : nx2x2 real valued Numpy array
-        - PT-error : nx2x2 real valued Numpy array
-
-    """
-
-    if not isinstance(edi_object, MTedi.Edi):
-        raise MTex.MTpyError_EDI('Input argument is not an instance of the Edi class')
-    p = PhaseTensor(edi_object=edi_object)
-
-    pt_array = p.pt
-
-    pterr_array = p.pterr
-
-    return pt_array, pterr_array
-
-
-def edi_file2pt(filename):
-    """
-        Calculate Phase Tensor from Edi-file (incl. uncertainties)
-
-        Input:
-        - Edi-file : full path to the Edi-file
-
-        Return:
-        - PT object
-
-    """
-    # Return:
-    # - PT : nx2x2 real valued Numpy array
-    # - PT-error : nx2x2 real valued Numpy array
-
-    # """
-
-    e = MTedi.Edi()
-    e.readfile(filename)
-
-    p = PhaseTensor(z_object=e.Z)
-
-    # pt_array = p.pt
-
-    # pterr_array = p.pterr
-
-    # return pt_array, pterr_array
-
-    return p
+    return pt_array, pt_err_array
\ No newline at end of file
diff --git a/mtpy/analysis/zinvariants.py b/mtpy/analysis/zinvariants.py
index eace6360f..ef0eb4859 100644
--- a/mtpy/analysis/zinvariants.py
+++ b/mtpy/analysis/zinvariants.py
@@ -151,7 +151,6 @@ def compute_invariants(self):
             **q**          : dependent variable suggesting dimensionality
             
         """
-        print("computing invariants")
         # get the length of z to initialize some empty arrays           
         nz = self.z.shape[0]
 
@@ -184,8 +183,8 @@ def compute_invariants(self):
             ex = x1 * e1 - x2 * e2 - x3 * e3 + x4 * e4
 
             if ex == 0.0:
-                print('Could not compute invariants for {0:5e} Hz'.format(
-                       self.freq[ii]))
+                # print('Could not compute invariants for {0:5e} Hz'.format(
+                #        self.freq[ii]))
                 self.inv1[ii] = np.nan
                 self.inv2[ii] = np.nan
                 self.inv3[ii] = np.nan
diff --git a/mtpy/core/mt.py b/mtpy/core/mt.py
index bde853484..b97b316fa 100644
--- a/mtpy/core/mt.py
+++ b/mtpy/core/mt.py
@@ -302,7 +302,7 @@ def rotation_angle(self, theta_r):
         self._Tipper.rotate(theta_r)
         self.pt.rotate(theta_r)
 
-        print(("Rotated Z, Tipper, Phase Tensor and Zinvariants by"
+        print(("Rotated Z, Tipper, Phase Tensor and Zinvariants by "
                "{0:.3f} degrees".format(self._rotation_angle)))
 
     @Z.setter
diff --git a/mtpy/core/z.py b/mtpy/core/z.py
index 60208c1da..cd22cc8cb 100644
--- a/mtpy/core/z.py
+++ b/mtpy/core/z.py
@@ -111,21 +111,12 @@ def compute_resistivity_phase(self, z_array=None, z_err_array=None,
             for idx_f in range(self.freq.size):
                 for ii in range(2):
                     for jj in range(2):
-#                        r_err, phi_err = MTcc.z_error2r_phi_error(
-#                            np.real(self._z[idx_f, ii, jj]),
-#                            self._z_err[idx_f, ii, jj],
-#                            np.imag(self._z[idx_f, ii, jj]),
-#                            self._z_err[idx_f, ii, jj])
-
                         r_err, phi_err = MTcc.z_error2r_phi_error(
                                 self._z[idx_f, ii, jj].real,
                                 self._z[idx_f, ii, jj].imag,
                                 self._z_err[idx_f, ii, jj])
                         self._resistivity_err[idx_f, ii, jj] = \
                             self._resistivity[idx_f, ii, jj] * r_err
-#                        self._resistivity_err[idx_f, ii, jj] = \
-#                            0.4 * np.abs(self._z[idx_f, ii, jj]) / \
-#                            self.freq[idx_f] * r_err
                         self._phase_err[idx_f, ii, jj] = phi_err
 
     def set_res_phase(self, res_array, phase_array, freq, res_err_array=None,
@@ -613,12 +604,12 @@ def rotate(self, alpha):
 
             if self.z_err is not None:
                 z_rot[idx_freq], z_err_rot[idx_freq] = \
-                    MTcc.rotatematrix_incl_errors(self.z[idx_freq, :, :],
+                    MTcc.rotate_matrix_with_errors(self.z[idx_freq, :, :],
                                                   angle,
                                                   self.z_err[idx_freq, :, :])
             else:
                 z_rot[idx_freq], z_err_rot = \
-                    MTcc.rotatematrix_incl_errors(self.z[idx_freq, :, :],
+                    MTcc.rotate_matrix_with_errors(self.z[idx_freq, :, :],
                                                   angle)
 
         self.z = z_rot
@@ -1521,12 +1512,12 @@ def rotate(self, alpha):
 
             if self.tipper_err is not None:
                 tipper_rot[idx_freq], tipper_err_rot[idx_freq] = \
-                    MTcc.rotatevector_incl_errors(self.tipper[idx_freq, :, :],
+                    MTcc.rotate_vector_with_errors(self.tipper[idx_freq, :, :],
                                                   angle,
                                                   self.tipper_err[idx_freq, :, :])
             else:
                 tipper_rot[idx_freq], tipper_err_rot = \
-                    MTcc.rotatevector_incl_errors(self.tipper[idx_freq, :, :],
+                    MTcc.rotate_vector_with_errors(self.tipper[idx_freq, :, :],
                                                   angle)
 
         self.tipper = tipper_rot
diff --git a/mtpy/imaging/phase_tensor_maps.py b/mtpy/imaging/phase_tensor_maps.py
index b4c6af295..57cfab6b8 100644
--- a/mtpy/imaging/phase_tensor_maps.py
+++ b/mtpy/imaging/phase_tensor_maps.py
@@ -367,6 +367,10 @@ def __init__(self, **kwargs):
         :param kwargs: keyword-value pairs
         """
         super(PlotPhaseTensorMaps, self).__init__(**kwargs)
+        
+        self.fig = None
+        self.ax_ptm = None
+        self.ax_cb = None
 
         fn_list = kwargs.pop('fn_list', None)
         z_object_list = kwargs.pop('z_object_list', None)
@@ -382,233 +386,134 @@ def __init__(self, **kwargs):
                                        mt_object_list=mt_object_list)
 
         # set the freq to plot
-        self.plot_freq = kwargs.pop('plot_freq', 1.0)
-        self.ftol = kwargs.pop('ftol', 0.1)
-        self.interpolate = kwargs.pop('interpolate', True)
+        self.plot_freq =  1.0
+        self.ftol = 0.1
+        self.interpolate = True
         # read in map scale
-        self.mapscale = kwargs.pop('mapscale', 'deg')
+        self.mapscale = 'deg'
         # map background image
-        self.background_image = kwargs.pop('background_image', None)
-        self.bimg_band = kwargs.pop('bimg_band', None)
-        self.bimg_cmap = kwargs.pop('bimg_cmap', 'viridis')
+        self.background_image = None
+        self.bimg_band = None
+        self.bimg_cmap = 'viridis'
 
         # --> set the ellipse properties -------------------
         # set default size to 2
         if self.mapscale == 'deg':
-            self.ellipse_size = kwargs.pop('ellipse_size', .05)
-            self.arrow_size = kwargs.pop('arrow_size', .05)
-            self.arrow_head_length = kwargs.pop('arrow_head_length', .005)
-            self.arrow_head_width = kwargs.pop('arrow_head_width', .005)
-            self.arrow_lw = kwargs.pop('arrow_lw', .0005)
-            self.xpad = kwargs.pop('xpad', .05)
-            self.ypad = kwargs.pop('xpad', .05)
+            self.ellipse_size = .05
+            self.arrow_size = .05
+            self.arrow_head_length = .005
+            self.arrow_head_width = .005
+            self.arrow_lw = .0005
+            self.xpad = .05
+            self.ypad = .05
 
         elif self.mapscale == 'm':
-            self.ellipse_size = kwargs.pop('ellipse_size', 500)
-            self.arrow_size = kwargs.pop('arrow_size', 500)
-            self.arrow_head_length = kwargs.pop('arrow_head_length', 50)
-            self.arrow_head_width = kwargs.pop('arrow_head_width', 50)
-            self.arrow_lw = kwargs.pop('arrow_lw', 5)
-            self.xpad = kwargs.pop('xpad', 500)
-            self.ypad = kwargs.pop('xpad', 500)
+            self.ellipse_size = 500
+            self.arrow_size = 500
+            self.arrow_head_length = 50
+            self.arrow_head_width = 50
+            self.arrow_lw = 5
+            self.xpad = 500
+            self.ypad = 500
 
         elif self.mapscale == 'km':
-            self.ellipse_size = kwargs.pop('ellipse_size', .5)
-            self.arrow_size = kwargs.pop('arrow_size', .5)
-            self.arrow_head_length = kwargs.pop('arrow_head_length', .05)
-            self.arrow_head_width = kwargs.pop('arrow_head_width', .05)
-            self.arrow_lw = kwargs.pop('arrow_lw', .005)
-            self.xpad = kwargs.pop('xpad', .5)
-            self.ypad = kwargs.pop('xpad', .5)
+            self.ellipse_size = .5
+            self.arrow_size = .5
+            self.arrow_head_length = .05
+            self.arrow_head_width = .05
+            self.arrow_lw = .005
+            self.xpad = .5
+            self.ypad = .5
 
-        self.minorticks_on = kwargs.pop('minorticks_on', True)
+        self.minorticks_on = True
 
         # --> set colorbar properties---------------------------------
         # set orientation to horizontal
-        cb_dict = kwargs.pop('cb_dict', {})
-        try:
-            self.cb_orientation = cb_dict['orientation']
-        except KeyError:
-            self.cb_orientation = 'vertical'
-
-        # set the position to middle outside the plot
-        try:
-            self.cb_position = cb_dict['position']
-        except KeyError:
-            self.cb_position = None
+        self.cb_orientation = 'vertical'
+        self.cb_position = None
 
         # --> set plot properties ------------------------------
         # set some of the properties as attributes much to Lars' discontent
-        self.fig_num = kwargs.pop('fig_num', 1)
-        self.plot_num = kwargs.pop('plot_num', 1)
-        self.plot_style = kwargs.pop('plot_style', '1')
-        self.plot_title = kwargs.pop('plot_title', None)
-        self.fig_dpi = kwargs.pop('fig_dpi', 300)
+        self.fig_num = 1
+        self.plot_title = 'Phase Tensor Map for '
+        self.fig_dpi = 300
 
-        self.tscale = kwargs.pop('tscale', 'period')
-        self.fig_size = kwargs.pop('fig_size', [8, 8])
-        self.tickstrfmt = kwargs.pop('tickstrfmt', '%.2f')
+        self.tscale = 'period'
+        self.fig_size = None
+        self.tickstrfmt = '%.2f'
 
-        self.font_size = kwargs.pop('font_size', 7)
+        self.font_size = 7
 
-        self.ref_ax_loc = kwargs.pop('ref_ax_loc', (.85, .1, .1, .1))
+        self.ref_ax_loc = (.85, .1, .1, .1)
         # if rotation angle is an int or float make an array the length of
         # mt_list for plotting purposes
 
-        self._rot_z = kwargs.pop('rot_z', 0)
-        if isinstance(self._rot_z, float) or isinstance(self._rot_z, int):
-            self._rot_z = np.array([self._rot_z] * len(self.mt_list))
-
-        # if the rotation angle is an array for rotation of different
-        # freq than repeat that rotation array to the len(mt_list)
-        elif isinstance(self._rot_z, np.ndarray):
-            if self._rot_z.shape[0] != len(self.mt_list):
-                self._rot_z = np.repeat(self._rot_z, len(self.mt_list))
-
-        else:
-            pass
-
         # --> set induction arrow properties -------------------------------
-        self.plot_tipper = kwargs.pop('plot_tipper', 'n')
+        self.plot_tipper = 'n'
 
         # --> set arrow legend properties -------------------------------
-        arrow_legend_dict = kwargs.pop('arrow_legend_dict', {})
 
-        self.arrow_legend_position = arrow_legend_dict.pop('position', 'lower right')
-
-        # set x-border pad
-        self.arrow_legend_xborderpad = arrow_legend_dict.pop('xborderpad', .2)
-
-        # set y-border pad
-        self.arrow_legend_yborderpad = arrow_legend_dict.pop('yborderpad', .2)
-
-        # set font pad
-        self.arrow_legend_fontpad = arrow_legend_dict.pop('fontpad', .05)
-
-        # set font properties
-        self.arrow_legend_fontdict = arrow_legend_dict.pop('fontdict', {'size': self.font_size,
-                                                                        'weight': 'bold'})
+        self.arrow_legend_position = 'lower right'
+        self.arrow_legend_xborderpad = .2
+        self.arrow_legend_yborderpad = .2
+        self.arrow_legend_fontpad = .05
+        self.arrow_legend_fontdict = {'size': self.font_size, 
+                                      'weight': 'bold'}
 
         # --> set a central reference point
-        self.plot_reference_point = kwargs.pop('reference_point', (0, 0))
+        self.plot_reference_point = (0, 0)
 
         # --> set station name properties
-        station_dict = kwargs.pop('station_dict', None)
-        if station_dict is not None:
-            self.station_id = station_dict.pop('id', (0, 2))
+        self.plot_station = False
+        self.station_id = (0, 2)
+        
 
-            # set spacing of station name and ellipse
-            self.station_pad = station_dict.pop('pad', .0005)
+        # set spacing of station name and ellipse
+        self.station_pad = .0005
 
-            # set font properties of the station label
-            self.station_font_dict = station_dict.pop('font_dict', {'size': self.font_size,
-                                                                    'weight': 'bold'})
+        # set font properties of the station label
+        self.station_font_dict = {'size': self.font_size, 'weight': 'bold'}
 
-        self.plot_yn = kwargs.pop('plot_yn', 'y')
-        self.save_fn = kwargs.pop('save_fn', "/c/tmp")
+        self.plot_yn = 'y'
+        self.save_fn = "/c/tmp"
 
-        # By this stage all keyword arguments meant to be set as class properties will have
-        # been processed. Popping all class properties that still exist in kwargs
-        self.kwargs = kwargs
-        for key in vars(self):
-            self.kwargs.pop(key, None)
+        for key, value in kwargs.items():
+            setattr(self, key, value)
 
         # --> plot if desired ------------------------
         if self.plot_yn == 'y':
             self.plot()
-            # self.save_figure(self.save_fn, file_format='png')
 
-    # ---need to rotate data on setting rotz
-    def _set_rot_z(self, rot_z):
+    #---need to rotate data on setting rotz
+    @property
+    def rotation_angle(self):
+        return self._rotation_angle
+    
+    @rotation_angle.setter
+    def rotation_angle(self, value):
         """
-        need to rotate data when setting z
+        only a single value is allowed
         """
-
-        # if rotation angle is an int or float make an array the length of
-        # mt_list for plotting purposes
-        if isinstance(rot_z, float) or isinstance(rot_z, int):
-            self._rot_z = np.array([rot_z] * len(self.mt_list))
-
-        # if the rotation angle is an array for rotation of different
-        # freq than repeat that rotation array to the len(mt_list)
-        elif isinstance(rot_z, np.ndarray):
-            if rot_z.shape[0] != len(self.mt_list):
-                self._rot_z = np.repeat(rot_z, len(self.mt_list))
-
-        else:
-            pass
-
-        for ii, mt in enumerate(self.mt_list):
-            mt.rot_z = self._rot_z[ii]
-
-    def _get_rot_z(self):
-        return self._rot_z
-
-    rot_z = property(fget=_get_rot_z, fset=_set_rot_z,
-                     doc="""rotation angle(s)""")
-
-    # -----------------------------------------------
-    # The main plot method for this module
-    # -----------------------------------------------
-    def plot(self, fig=None, save_path=None, show=True,
-             raster_dict={'lons': [], 'lats': [],
-                          'vals': [], 'levels': 50, 'cmap': 'rainbow',
-                          'cbar_title': 'Arbitrary units',
-                          'cbar_position': None}):
+        if value == 0:
+            return
+        for ii, mt_obj in enumerate(self.mt_list):
+            mt_obj.rotation_angle = value
+            
+        self._rotation_angle = value
+        
+    def fold_strike(self, strike):
         """
-        Plots the phase tensor map.
-        :param fig: optional figure object
-        :param save_path: path to folder for saving plots
-        :param show: show plots if True
-        :param raster_dict: Plotting of raster data is currently only supported when mapscale='deg'.
-                            This parameter is a dictionary of parameters for plotting raster data,
-                            on top of which phase tensor data are plotted. 'lons', 'lats' and 'vals'
-                            are one dimensional lists (or numpy arrays) for longitudes, latitudes
-                            and corresponding values, respectively. 'levels', 'cmap' and 'cbar_title'
-                            are the number of levels to be used in the colormap, the colormap and its
-                            title, respectively.
+        
         """
-
-        # set position properties for the plot
-        plt.rcParams['font.size'] = self.font_size
-        plt.rcParams['figure.subplot.left'] = .1
-        plt.rcParams['figure.subplot.right'] = .98
-        plt.rcParams['figure.subplot.bottom'] = .1
-        plt.rcParams['figure.subplot.top'] = .93
-        plt.rcParams['figure.subplot.wspace'] = .55
-        plt.rcParams['figure.subplot.hspace'] = .70
-        # FZ: tweaks to make plot positioned better
-        # plt.rcParams['font.size']=self.font_size
-        # plt.rcParams['figure.subplot.left']=.1
-        # plt.rcParams['figure.subplot.right']=.90
-        # plt.rcParams['figure.subplot.bottom']=.2
-        # plt.rcParams['figure.subplot.top']=.90
-        # plt.rcParams['figure.subplot.wspace']=.70
-        # plt.rcParams['figure.subplot.hspace']=.70
-
-        lpfig = None
-        lpax = None
-        lpax2 = None
-        # make figure instance
-        if(fig is None):
-            self.fig = plt.figure(self.fig_num, figsize=self.fig_size, dpi=self.fig_dpi)
-            # self.fig = plt.figure(self.fig_num, dpi=self.fig_dpi)
-
-            # clear the figure if there is already one up
-            plt.clf()
-            lpfig = self.fig
-        else:
-            lpfig = fig
-        # end if
-
-        lpax = lpfig.add_subplot(1, 1, 1, aspect='equal')
-
-        # plt.locator_params(axis='x', nbins=3)  # control number of ticks in axis (nbins ticks)
-        plt.xticks(rotation='vertical')  # FZ: control tick rotation=30 not that good
-
-        # get the reference point
-        refpoint = self.plot_reference_point
-
+        strike = strike % 180
+        strike[np.where(strike > 90)] -= 180
+        
+        return strike
+    
+    def add_raster(self, raster_dict):
+        """
+        """
+        
         # plot raster data if provided and if mapscale is 'deg'
         if(len(raster_dict['lons']) and self.mapscale == 'deg'):
             lons = np.array(raster_dict['lons'])
@@ -620,24 +525,28 @@ def plot(self, fig=None, save_path=None, show=True,
             else:
                 vals = np.array(raster_dict['vals'])
 
-            assert len(lons) == len(lats) == len(vals), 'Lons, Lats and Vals must all have the same length'
+            if not len(lons) == len(lats) == len(vals):
+                msg= 'Raster Lons, Lats and Vals must all have the same length'
+                _logger.error(msg)
 
-            lons -= refpoint[0]
-            lats -= refpoint[1]
+            lons -= self.plot_reference_point[0]
+            lats -= self.plot_reference_point[1]
             levels = raster_dict.pop('levels', 50)
             cmap = raster_dict.pop('cmap', 'rainbow')
             cbar_title = raster_dict.pop('cbar_title', 'Arbitrary Units')
             triangulation = tri.Triangulation(lons, lats)
-            cbinfo = lpax.tricontourf(triangulation, vals,
-                                      levels=np.linspace(vals.min(), vals.max(), levels),
+            cbinfo = self.ax_ptm.tricontourf(triangulation, vals,
+                                      levels=np.linspace(vals.min(), 
+                                                         vals.max(), 
+                                                         levels),
                                       cmap=cmap)
             if raster_dict['cbar_position'] is not None:
                 cbax = self.fig.add_axes(raster_dict['cbar_position'])
             else:
-                cbax, kw = mcb.make_axes(lpax,
+                cbax, kw = mcb.make_axes(self.ax_ptm,
                                          orientation=self.cb_orientation,
                                          shrink=.35)
-            cbar = lpfig.colorbar(cbinfo, cbax)
+            cbar = self.fig.colorbar(cbinfo, cbax)
 
             if(self.cb_orientation == 'horizontal'):
                 cbar.ax.set_xlabel(cbar_title)
@@ -652,25 +561,257 @@ def plot(self, fig=None, save_path=None, show=True,
                 cbar.ax.tick_params(axis='y', direction='in')
 
         # end if
+        
+    def add_tipper(self, tipper_obj, plot_x, plot_y):
+        """
+        
+        :param tipper: DESCRIPTION
+        :type tipper: TYPE
+        :return: DESCRIPTION
+        :rtype: TYPE
+
+        """
+
+        # make some local parameters for easier typing
+        ascale = self.arrow_size
+        adir = self.arrow_direction * np.pi
+
+        # plot real tipper
+        ti = tipper_obj
+        if ti is None:
+            return
+
+        if self.plot_tipper == 'yri' or self.plot_tipper == 'yr':
+            if ti.mag_real[self.jj] <= self.arrow_threshold:
+                txr = ti.mag_real[self.jj] * ascale * \
+                    np.sin((ti.angle_real[self.jj]) * np.pi / 180 + adir)
+                tyr = ti.mag_real[self.jj] * ascale * \
+                    np.cos((ti.angle_real[self.jj]) * np.pi / 180 + adir)
+
+                self.ax_ptm.arrow(plot_x,
+                                  plot_y,
+                                  txr,
+                                  tyr,
+                                  width=self.arrow_lw,
+                                  facecolor=self.arrow_color_real,
+                                  edgecolor=self.arrow_color_real,
+                                  length_includes_head=False,
+                                  head_width=self.arrow_head_width,
+                                  head_length=self.arrow_head_length)
+            else:
+                pass
+
+        # plot imaginary tipper
+        if self.plot_tipper == 'yri' or self.plot_tipper == 'yi':
+            if ti.mag_imag[self.jj] <= self.arrow_threshold:
+                txi = ti.mag_imag[self.jj] * ascale * \
+                    np.sin((ti.angle_imag[self.jj]) * np.pi / 180 + adir)
+                tyi = ti.mag_imag[self.jj] * ascale * \
+                    np.cos((ti.angle_imag[self.jj]) * np.pi / 180 + adir)
+
+                self.ax_ptm.arrow(plot_x,
+                                  plot_y,
+                                  txi,
+                                  tyi,
+                                  width=self.arrow_lw,
+                                  facecolor=self.arrow_color_imag,
+                                  edgecolor=self.arrow_color_imag,
+                                  length_includes_head=False,
+                                  head_width=self.arrow_head_width,
+                                  head_length=self.arrow_head_length)
+    
+    def get_xy(self, mt_obj):
+        """
+        """
+        # if map scale is lat lon set parameters
+        if self.mapscale == 'deg':
+            plot_x = mt_obj.lon - self.plot_reference_point[0]
+            plot_y = mt_obj.lat - self.plot_reference_point[1]
+
+        # if map scale is in meters easting and northing
+        elif self.mapscale in ['km', 'm']:
+            if self.mapscale == 'km':
+                scale = 1./1000
+            else:
+                scale = 1.
+            east, north, zone = gis_tools.project_point_ll2utm(mt_obj.lat,
+                                                               mt_obj.lon)
+
+            # set the first point read in as a refernce other points
+            if not hasattr(self._utm_zone_01):
+                self._utm_zone_01 = zone
+                plot_x = east - self.plot_reference_point[0]
+                plot_y = north - self.plot_reference_point[1]
+                
+            else:
+                # check to make sure the zone is the same this needs
+                # to be more rigorously done
+                if self._utm_zone_01 != zone:
+                    print('Zone change at station ' + mt_obj.station)
+                    if zone1[0:2] == zone[0:2]:
+                        pass
+                    elif int(zone1[0:2]) < int(zone[0:2]):
+                        east += 500000
+                    else:
+                        east -= -500000
+                    plot_x = east - self.plot_reference_point[0]
+                    plot_y = north - self.plot_reference_point[1]
+                else:
+                    plot_x = east - self.plot_reference_point[0]
+                    plot_y = north - self.plot_reference_point[1]
+                
+            plot_x *= scale
+            plot_y *= scale
+
+        else:
+            raise NameError('mapscale not recognized')
+            
+        return plot_x, plot_y
+    
+    def get_pt_ellipse(self, pt, plot_x, plot_y):
+        """
+        """
+
+        # --> set local variables
+        phimin = np.nan_to_num(pt.phimin[self.jj])
+        phimax = np.nan_to_num(pt.phimax[self.jj])
+        eangle = np.nan_to_num(pt.azimuth[self.jj])
+
+        if self.ellipse_cmap == 'mt_seg_bl2wh2rd':
+            bounds = np.arange(self.ellipse_range[0], 
+                               self.ellipse_range[1] + self.ellipse_range[2],
+                               self.ellipse_range[2])
+            nseg = float((self.ellipse_range[1] - self.ellipse_range[0]) / \
+                         (2 * self.ellipse_range[2]))
+
+        # get the properties to color the ellipses by
+        if self.ellipse_colorby == 'phiminang' or \
+                self.ellipse_colorby == 'phimin':
+            colorarray = pt.phimin[self.jj]
+
+        elif self.ellipse_colorby == 'phimax':
+            colorarray = pt.phimax[self.jj]
+
+        elif self.ellipse_colorby == 'phidet':
+            colorarray = np.sqrt(abs(pt.det[self.jj])) * (180 / np.pi)
+
+        elif self.ellipse_colorby == 'skew' or \
+                self.ellipse_colorby == 'skew_seg':
+            colorarray = pt.beta[self.jj]
+
+        elif self.ellipse_colorby == 'normalized_skew' or \
+                self.ellipse_colorby == 'normalized_skew_seg':
+            colorarray = 2 * pt.beta[self.jj]
+
+        elif self.ellipse_colorby == 'ellipticity':
+            colorarray = pt.ellipticity[self.jj]
+            
+        elif self.ellipse_colorby in ['strike', 'azimuth']:
+            colorarray = pt.azimuth[self.jj] % 180
+            if colorarray > 90:
+                colorarray -= 180
+            self.ellipse_range = (-90, 90)
+
+        else:
+            raise NameError(self.ellipse_colorby + ' is not supported')
+
+        # --> get ellipse properties
+        # if the ellipse size is not physically correct make it a dot
+        if phimax == 0 or phimax > 100 or phimin == 0 or phimin > 100:
+            eheight = .0000001 * self.ellipse_size
+            ewidth = .0000001 * self.ellipse_size
+        else:
+            scaling = self.ellipse_size / phimax
+            eheight = phimin * scaling
+            ewidth = phimax * scaling
+
+        # make an ellipse
+        ellipd = patches.Ellipse((plot_x, plot_y),
+                                 width=ewidth,
+                                 height=eheight,
+                                 angle=90 - eangle,
+                                 lw=self.lw)
+
+        # get ellipse color
+        if self.ellipse_cmap.find('seg') > 0:
+            ellipd.set_facecolor(mtcl.get_plot_color(colorarray,
+                                                     self.ellipse_colorby,
+                                                     self.ellipse_cmap,
+                                                     self.ellipse_range[0],
+                                                     self.ellipse_range[1],
+                                                     bounds=bounds))
+        else:
+            ellipd.set_facecolor(mtcl.get_plot_color(colorarray,
+                                                     self.ellipse_colorby,
+                                                     self.ellipse_cmap,
+                                                     self.ellipse_range[0],
+                                                     self.ellipse_range[1]))
+            
+        return ellipd
+        
+
+    # -----------------------------------------------
+    # The main plot method for this module
+    # -----------------------------------------------
+    def plot(self, fig=None, save_path=None, show=True,
+             raster_dict={'lons': [], 'lats': [],
+                          'vals': [], 'levels': 50, 'cmap': 'rainbow',
+                          'cbar_title': 'Arbitrary units',
+                          'cbar_position': None}):
+        """
+        Plots the phase tensor map.
+        :param fig: optional figure object
+        :param save_path: path to folder for saving plots
+        :param show: show plots if True
+        :param raster_dict: Plotting of raster data is currently only 
+                            supported when mapscale='deg'.
+                            This parameter is a dictionary of parameters for 
+                            plotting raster data,
+                            on top of which phase tensor data are plotted.
+                            'lons', 'lats' and 'vals'  are one dimensional 
+                            lists (or numpy arrays) for longitudes, latitudes
+                            and corresponding values, respectively. 'levels', 
+                            'cmap' and 'cbar_title' are the number of levels
+                            to be used in the colormap, the colormap and its
+                            title, respectively.
+        """
+
+        # set position properties for the plot
+        plt.rcParams['font.size'] = self.font_size
+
+        # make figure instance
+        if fig is not None:
+            self.fig = fig
+        else:
+            self.fig = plt.figure(self.fig_num, figsize=self.fig_size,
+                                  dpi=self.fig_dpi)
+        plt.clf()
+
+        self.ax_ptm = self.fig.add_subplot(1, 1, 1, aspect='equal')
+
+        # plt.locator_params(axis='x', nbins=3)  
+        # control number of ticks in axis (nbins ticks)
+        # FZ: control tick rotation=30 not that good
+        plt.xticks(rotation='vertical')  
+
+        self.add_raster(raster_dict)
 
-        # set some local parameters
-        es = float(self.ellipse_size)
-        cmap = self.ellipse_cmap
-        ckmin = float(self.ellipse_range[0])
-        ckmax = float(self.ellipse_range[1])
         try:
-            ckstep = float(self.ellipse_range[2])
+            step = float(self.ellipse_range[2])
         except IndexError:
-            if cmap == 'mt_seg_bl2wh2rd':
+            if self.ellipse_cmap == 'mt_seg_bl2wh2rd':
                 raise ValueError('Need to input range as (min, max, step)')
             else:
-                ckstep = 3
-        nseg = float((ckmax - ckmin) / (2 * ckstep))
+                step = 3
+        nseg = float((self.ellipse_range[1] - self.ellipse_range[0]) / \
+                     (2 * step))
         ck = self.ellipse_colorby
 
         # --> set the bounds on the segmented colormap
-        if cmap == 'mt_seg_bl2wh2rd':
-            bounds = np.arange(ckmin, ckmax + ckstep, ckstep)
+        if self.ellipse_cmap == 'mt_seg_bl2wh2rd':
+            bounds = np.arange(self.ellipse_range[0], 
+                               self.ellipse_range[1] + step,
+                               step)
 
         # set tick parameters depending on the mapscale
         if self.mapscale == 'deg':
@@ -681,302 +822,119 @@ def plot(self, fig=None, save_path=None, show=True,
 
         # make some empty arrays
         elliplist = []
-        latlist = np.zeros(len(self.mt_list))
-        lonlist = np.zeros(len(self.mt_list))
         self.plot_xarr = np.zeros(len(self.mt_list))
         self.plot_yarr = np.zeros(len(self.mt_list))
 
-        for ii, mt in enumerate(self.mt_list):
+        for ii, mt_obj in enumerate(self.mt_list):
 
-            newZ = None
-            newTipper = None
+            new_z_obj = None
+            new_tipper_obj = None
             fidx = 0
-            if(self.interpolate):
-                newZ, newTipper = mt.interpolate([self.plot_freq], bounds_error=False)
+            if self.interpolate:
+                new_z_obj, new_tipper_obj = mt_obj.interpolate([self.plot_freq],
+                                                 bounds_error=False)
             else:
-                fidx = np.argmin(np.fabs(mt.Z.freq - self.plot_freq))
+                fidx = np.argmin(np.fabs(mt_obj.Z.freq - self.plot_freq))
 
-            if((not self.interpolate and np.fabs(mt.Z.freq[fidx] - self.plot_freq) < self.ftol) or
-                    (self.interpolate)):
+            if((not self.interpolate and \
+                np.fabs(mt_obj.Z.freq[fidx] - self.plot_freq) < self.ftol) or \
+               (self.interpolate)):
 
                 self.jj = fidx
                 jj = fidx
 
                 # get phase tensor
                 if(not self.interpolate):
-                    pt = mt.pt
+                    pt = mt_obj.pt
                 else:
-                    newZ.compute_resistivity_phase()
-                    pt = MTpt.PhaseTensor(z_object=newZ)
-
-                # if map scale is lat lon set parameters
-                if self.mapscale == 'deg':
-                    latlist[ii] = mt.lat
-                    lonlist[ii] = mt.lon
-                    plotx = mt.lon - refpoint[0]
-                    ploty = mt.lat - refpoint[1]
-
-                # if map scale is in meters easting and northing
-                elif self.mapscale == 'm':
-                    east, north, zone = gis_tools.project_point_ll2utm(mt.lat,
-                                                                       mt.lon)
-
-                    # set the first point read in as a refernce other points
-                    if ii == 0:
-                        zone1 = zone
-                        plotx = east - refpoint[0]
-                        ploty = north - refpoint[1]
-
-                    # read in all the other point
-                    else:
-                        # check to make sure the zone is the same this needs
-                        # to be more rigorously done
-                        if zone1 != zone:
-                            print('Zone change at station ' + mt.station)
-                            if zone1[0:2] == zone[0:2]:
-                                pass
-                            elif int(zone1[0:2]) < int(zone[0:2]):
-                                east += 500000
-                            else:
-                                east -= -500000
-                            latlist[ii] = north - refpoint[1]
-                            lonlist[ii] = east - refpoint[0]
-                            plotx = east - refpoint[0]
-                            ploty = north - refpoint[1]
-                        else:
-                            latlist[ii] = north - refpoint[1]
-                            lonlist[ii] = east - refpoint[0]
-                            plotx = east - refpoint[0]
-                            ploty = north - refpoint[1]
-
-                # if mapscale is in km easting and northing
-                elif self.mapscale == 'km':
-                    east, north, zone = gis_tools.project_point_ll2utm(mt.lat,
-                                                                       mt.lon)
-                    if ii == 0:
-                        zone1 = zone
-                        plotx = (east - refpoint[0]) / 1000.
-                        ploty = (north - refpoint[1]) / 1000.
-
-                    else:
-                        if zone1 != zone:
-                            print('Zone change at station ' + mt.station)
-                            if zone1[0:2] == zone[0:2]:
-                                pass
-                            elif int(zone1[0:2]) < int(zone[0:2]):
-                                east += 500000
-                            else:
-                                east -= 500000
-                            latlist[ii] = (north - refpoint[1]) / 1000.
-                            lonlist[ii] = (east - refpoint[0]) / 1000.
-                            plotx = (east - refpoint[0]) / 1000.
-                            ploty = (north - refpoint[1]) / 1000.
-                        else:
-                            latlist[ii] = (north - refpoint[1]) / 1000.
-                            lonlist[ii] = (east - refpoint[0]) / 1000.
-                            plotx = (east - refpoint[0]) / 1000.
-                            ploty = (north - refpoint[1]) / 1000.
-                else:
-                    raise NameError('mapscale not recognized')
+                    new_z_obj.compute_resistivity_phase()
+                    pt = MTpt.PhaseTensor(z_object=new_z_obj)
 
+                plot_x, plot_y = self.get_xy(mt_obj)
                 # put the location of each ellipse into an array in x and y
-                self.plot_xarr[ii] = plotx
-                self.plot_yarr[ii] = ploty
-
-                # --> set local variables
-                phimin = np.nan_to_num(pt.phimin[jj])
-                phimax = np.nan_to_num(pt.phimax[jj])
-                eangle = np.nan_to_num(pt.azimuth[jj])
-
-                # output to csv file:
-                # print('OUTCSV', mt.station, plotx, ploty, phimin, phimax, eangle)
-
-                if cmap == 'mt_seg_bl2wh2rd':
-                    bounds = np.arange(ckmin, ckmax + ckstep, ckstep)
-                    nseg = float((ckmax - ckmin) / (2 * ckstep))
-
-                # get the properties to color the ellipses by
-                if self.ellipse_colorby == 'phiminang' or \
-                        self.ellipse_colorby == 'phimin':
-                    colorarray = pt.phimin[jj]
-
-                elif self.ellipse_colorby == 'phimax':
-                    colorarray = pt.phimax[jj]
-
-                elif self.ellipse_colorby == 'phidet':
-                    colorarray = np.sqrt(abs(pt.det[jj])) * (180 / np.pi)
-
-                elif self.ellipse_colorby == 'skew' or \
-                        self.ellipse_colorby == 'skew_seg':
-                    colorarray = pt.beta[jj]
-
-                elif self.ellipse_colorby == 'normalized_skew' or \
-                        self.ellipse_colorby == 'normalized_skew_seg':
-                    colorarray = 2 * pt.beta[jj]
-
-                elif self.ellipse_colorby == 'ellipticity':
-                    colorarray = pt.ellipticity[jj]
-
-                else:
-                    raise NameError(self.ellipse_colorby + ' is not supported')
-
-                # --> get ellipse properties
-                # if the ellipse size is not physically correct make it a dot
-                if phimax == 0 or phimax > 100 or phimin == 0 or phimin > 100:
-                    eheight = .0000001 * es
-                    ewidth = .0000001 * es
-                else:
-                    scaling = es / phimax
-                    eheight = phimin * scaling
-                    ewidth = phimax * scaling
-
-                # make an ellipse
-                ellipd = patches.Ellipse((plotx, ploty),
-                                         width=ewidth,
-                                         height=eheight,
-                                         angle=90 - eangle,
-                                         lw=self.lw,
-                                         **self.kwargs)
-
-                # get ellipse color
-                if cmap.find('seg') > 0:
-                    ellipd.set_facecolor(mtcl.get_plot_color(colorarray,
-                                                             self.ellipse_colorby,
-                                                             cmap,
-                                                             ckmin,
-                                                             ckmax,
-                                                             bounds=bounds))
-                else:
-                    ellipd.set_facecolor(mtcl.get_plot_color(colorarray,
-                                                             self.ellipse_colorby,
-                                                             cmap,
-                                                             ckmin,
-                                                             ckmax))
+                self.plot_xarr[ii] = plot_x
+                self.plot_yarr[ii] = plot_y
+                
+                ellipd = self.get_pt_ellipse(pt, plot_x, plot_y)
 
                 # ==> add ellipse to the plot
                 elliplist.append(ellipd)
-                lpax.add_artist(ellipd)
+                self.ax_ptm.add_artist(ellipd)
 
                 # -----------Plot Induction Arrows---------------------------
                 if self.plot_tipper.find('y') == 0:
-
-                    # get tipper
-
-                    if mt.Tipper.tipper is None:
-                        mt.Tipper.tipper = np.zeros((len(mt.period), 1, 2),
-                                                    dtype='complex')
-
-                    # make some local parameters for easier typing
-                    ascale = self.arrow_size
-                    adir = self.arrow_direction * np.pi
-
-                    # plot real tipper
-                    ti = None
-                    if(not self.interpolate):
-                        ti = mt.Tipper
+                    # get phase tensor
+                    if not self.interpolate:
+                        new_tipper_obj = mt_obj.Tipper
                     else:
-                        ti = newTipper
-                    # end if
-
-                    if self.plot_tipper == 'yri' or self.plot_tipper == 'yr':
-                        if ti.mag_real[jj] <= self.arrow_threshold:
-                            txr = ti.mag_real[jj] * ascale * \
-                                np.sin((ti.angle_real[jj]) * np.pi / 180 + adir)
-                            tyr = ti.mag_real[jj] * ascale * \
-                                np.cos((ti.angle_real[jj]) * np.pi / 180 + adir)
-
-                            lpax.arrow(plotx,
-                                       ploty,
-                                       txr,
-                                       tyr,
-                                       width=self.arrow_lw,
-                                       facecolor=self.arrow_color_real,
-                                       edgecolor=self.arrow_color_real,
-                                       length_includes_head=False,
-                                       head_width=self.arrow_head_width,
-                                       head_length=self.arrow_head_length)
-                        else:
-                            pass
-
-                    # plot imaginary tipper
-                    if self.plot_tipper == 'yri' or self.plot_tipper == 'yi':
-                        if ti.mag_imag[jj] <= self.arrow_threshold:
-                            txi = ti.mag_imag[jj] * ascale * \
-                                np.sin((ti.angle_imag[jj]) * np.pi / 180 + adir)
-                            tyi = ti.mag_imag[jj] * ascale * \
-                                np.cos((ti.angle_imag[jj]) * np.pi / 180 + adir)
-
-                            lpax.arrow(plotx,
-                                       ploty,
-                                       txi,
-                                       tyi,
-                                       width=self.arrow_lw,
-                                       facecolor=self.arrow_color_imag,
-                                       edgecolor=self.arrow_color_imag,
-                                       length_includes_head=False,
-                                       head_width=self.arrow_head_width,
-                                       head_length=self.arrow_head_length)
+                        new_tipper_obj.compute_mag_direction()
+                    
+                    if mt_obj.Tipper.tipper is None:
+                        continue
+                    self.add_tipper(new_tipper_obj, plot_x, plot_y)
 
                 # ------------Plot station name------------------------------
-                try:
-                    lpax.text(plotx,
-                              ploty + self.station_pad,
-                              mt.station[self.station_id[0]:self.station_id[1]],
-                              horizontalalignment='center',
-                              verticalalignment='baseline',
-                              fontdict=self.station_font_dict)
-                except AttributeError:
-                    pass
+                if self.plot_station:
+                    try:
+                        name = mt_obj.station[self.station_id[0]:self.station_id[1]]
+                        self.ax_ptm.text(plot_x,
+                                         plot_y + self.station_pad,
+                                         name,
+                                         horizontalalignment='center',
+                                         verticalalignment='baseline',
+                                         fontdict=self.station_font_dict)
+                    except AttributeError:
+                        pass
 
             # ==> print a message if couldn't find the freq
             else:
-                _logger.warn('Did not find {0:.5g} Hz for station {1}'.format(self.plot_freq, mt.station))
+                _logger.warn('Did not find {0:.5g} Hz for station {1}'.format(self.plot_freq, mt_obj.station))
 
         # --> set axes properties depending on map scale------------------------
         if self.mapscale == 'deg':
-            lpax.set_xlabel('Longitude',
+            self.ax_ptm.set_xlabel('Longitude',
                             fontsize=self.font_size,  # +2,
                             fontweight='bold')
-            lpax.set_ylabel('Latitude',
+            self.ax_ptm.set_ylabel('Latitude',
                             fontsize=self.font_size,  # +2,
                             fontweight='bold')
 
         elif self.mapscale == 'm':
-            lpax.set_xlabel('Easting (m)',
+            self.ax_ptm.set_xlabel('Easting (m)',
                             fontsize=self.font_size,  # +2,
                             fontweight='bold')
-            lpax.set_ylabel('Northing (m)',
+            self.ax_ptm.set_ylabel('Northing (m)',
                             fontsize=self.font_size,  # +2,
                             fontweight='bold')
 
         elif self.mapscale == 'km':
-            lpax.set_xlabel('Easting (km)',
+            self.ax_ptm.set_xlabel('Easting (km)',
                             fontsize=self.font_size,  # +2,
                             fontweight='bold')
-            lpax.set_ylabel('Northing (km)',
+            self.ax_ptm.set_ylabel('Northing (km)',
                             fontsize=self.font_size,  # +2,
                             fontweight='bold')
 
         # --> set plot limits
         #    need to exclude zero values from the calculation of min/max!!!!
-        lpax.set_xlim(self.plot_xarr[self.plot_xarr != 0.].min() - self.xpad,
+        self.ax_ptm.set_xlim(self.plot_xarr[self.plot_xarr != 0.].min() - self.xpad,
                       self.plot_xarr[self.plot_xarr != 0.].max() + self.xpad)
-        lpax.set_ylim(self.plot_yarr[self.plot_yarr != 0.].min() - self.xpad,
+        self.ax_ptm.set_ylim(self.plot_yarr[self.plot_yarr != 0.].min() - self.xpad,
                       self.plot_yarr[self.plot_xarr != 0.].max() + self.xpad)
 
         # BM: Now that we have the bounds of the axis, we can plot a
         # background image on the map.
         if self.background_image:
-            plot_geotiff_on_axes(self.background_image, lpax,
+            plot_geotiff_on_axes(self.background_image, self.ax_ptm,
                                  epsg_code=4326, band_number=self.bimg_band,
                                  cmap=self.bimg_cmap)
 
         # --> set tick label format
-        lpax.xaxis.set_major_formatter(FormatStrFormatter(self.tickstrfmt))
-        lpax.yaxis.set_major_formatter(FormatStrFormatter(self.tickstrfmt))
-#       lpax.set_xticklabels(np.round(self.plot_xarr, decimals=2),
+        self.ax_ptm.xaxis.set_major_formatter(FormatStrFormatter(self.tickstrfmt))
+        self.ax_ptm.yaxis.set_major_formatter(FormatStrFormatter(self.tickstrfmt))
+#       self.ax_ptm.set_xticklabels(np.round(self.plot_xarr, decimals=2),
 #                                rotation=45)
-        plt.setp(lpax.get_xticklabels(), rotation=45)
+        plt.setp(self.ax_ptm.get_xticklabels(), rotation=45)
 
         # --> set title in period or freq
         if self.tscale == 'period':
@@ -984,17 +942,13 @@ def plot(self, fig=None, save_path=None, show=True,
         else:
             titlefreq = '{0:.5g} (Hz)'.format(self.plot_freq)
 
-        if not self.plot_title:
-            lpax.set_title('Phase Tensor Map for ' + titlefreq,
-                           fontsize=self.font_size + 2, fontweight='bold')
-        else:
-            lpax.set_title(self.plot_title + titlefreq,
-                           fontsize=self.font_size + 2, fontweight='bold')
+        self.ax_ptm.set_title('{0} {1}'.format(self.plot_title, titlefreq),
+                              fontsize=self.font_size + 2, fontweight='bold')
 
         # --> plot induction arrow scale bar -----------------------------------
         if self.plot_tipper.find('y') == 0:
-            parrx = lpax.get_xlim()
-            parry = lpax.get_ylim()
+            parrx = self.ax_ptm.get_xlim()
+            parry = self.ax_ptm.get_ylim()
             try:
                 axpad = self.arrow_legend_xborderpad
             except AttributeError:
@@ -1008,7 +962,7 @@ def plot(self, fig=None, save_path=None, show=True,
             try:
                 txtpad = self.arrow_legend_fontpad
             except AttributeError:
-                txtpad = .25 * es
+                txtpad = .25 * self.ellipse_size
 
             # make arrow legend postion and arrows coordinates
             if self.arrow_legend_position == 'lower right':
@@ -1050,19 +1004,19 @@ def plot(self, fig=None, save_path=None, show=True,
             pty = 0
             # txy = pay + txtpad
 
-            lpax.arrow(pax,
-                       pay,
-                       ptx,
-                       pty,
-                       width=self.arrow_lw,
-                       facecolor=self.arrow_color_real,
-                       edgecolor=self.arrow_color_real,
-                       length_includes_head=False,
-                       head_width=self.arrow_head_width,
-                       head_length=self.arrow_head_length)
+            self.ax_ptm.arrow(pax,
+                              pay,
+                              ptx,
+                              pty,
+                              width=self.arrow_lw,
+                              facecolor=self.arrow_color_real,
+                              edgecolor=self.arrow_color_real,
+                              length_includes_head=False,
+                              head_width=self.arrow_head_width,
+                              head_length=self.arrow_head_length)
 
             # FZ: what is this '|T|=1'? and the horizontal line?
-            # lpax.text(txa,
+            # self.ax_ptm.text(txa,
             #              txy,
             #              '|T|=1',
             #              horizontalalignment='center',
@@ -1071,13 +1025,14 @@ def plot(self, fig=None, save_path=None, show=True,
         # END: if self.plot_tipper.find('yes') == 0 ---------------------------
 
         # make a grid with color lines
-        lpax.grid(True, alpha=.3, which='both', color=(0.5, 0.5, 0.5))
+        self.ax_ptm.grid(True, alpha=.3, which='both', color=(0.5, 0.5, 0.5))
+        self.ax_ptm.set_axisbelow(True)
         if self.minorticks_on:
             plt.minorticks_on()  # turn on minor ticks automatically
 
         # ==> make a colorbar with appropriate colors
         if self.cb_position is None:
-            lpax2, kw = mcb.make_axes(lpax,
+            self.ax_cb, kw = mcb.make_axes(self.ax_ptm,
                                       orientation=self.cb_orientation,
                                       shrink=.35)
             # FZ: try to fix colorbar h-position
@@ -1089,9 +1044,9 @@ def plot(self, fig=None, save_path=None, show=True,
             # self.ax2 = divider.append_axes("right", size="5%", pad=0.05)
 
         else:
-            lpax2 = lpfig.add_axes(self.cb_position)
+            self.ax_cb = self.fig.add_axes(self.cb_position)
 
-        if cmap == 'mt_seg_bl2wh2rd':
+        if self.ellipse_cmap == 'mt_seg_bl2wh2rd':
             # make a color list
             self.clist = [(cc, cc, 1)
                           for cc in np.arange(0, 1 + 1. / (nseg), 1. / (nseg))] + \
@@ -1102,26 +1057,26 @@ def plot(self, fig=None, save_path=None, show=True,
             mt_seg_bl2wh2rd = colors.ListedColormap(self.clist)
 
             # make bounds so that the middle is white
-            bounds = np.arange(ckmin - ckstep, ckmax + 2 * ckstep, ckstep)
+            bounds = np.arange(self.ellipse_range[0] - self.ellipse_range[2], self.ellipse_range[1] + 2 * self.ellipse_range[2], self.ellipse_range[2])
 
             # normalize the colors
             norms = colors.BoundaryNorm(bounds, mt_seg_bl2wh2rd.N)
 
             # make the colorbar
-            self.cb = mcb.ColorbarBase(lpax2,
+            self.cb = mcb.ColorbarBase(self.ax_cb,
                                        cmap=mt_seg_bl2wh2rd,
                                        norm=norms,
                                        orientation=self.cb_orientation,
                                        ticks=bounds[1:-1])
         else:
-            if cmap in list(mtcl.cmapdict.keys()):
-                cmap_input = mtcl.cmapdict[cmap]
+            if self.ellipse_cmap in list(mtcl.cmapdict.keys()):
+                cmap_input = mtcl.cmapdict[self.ellipse_cmap]
             else:
-                cmap_input = mtcl.cm.get_cmap(cmap)
-            self.cb = mcb.ColorbarBase(lpax2,
+                cmap_input = mtcl.cm.get_cmap(self.ellipse_cmap)
+            self.cb = mcb.ColorbarBase(self.ax_cb,
                                        cmap=cmap_input,  # mtcl.cmapdict[cmap],
-                                       norm=colors.Normalize(vmin=ckmin,
-                                                             vmax=ckmax),
+                                       norm=colors.Normalize(vmin=self.ellipse_range[0],
+                                                             vmax=self.ellipse_range[1]),
                                        orientation=self.cb_orientation)
 
         # label the color bar accordingly
@@ -1144,19 +1099,19 @@ def plot(self, fig=None, save_path=None, show=True,
         show_phi = False  # JingMingDuan does not want to show the black circle - it's not useful
         if show_phi is True:
             ref_ellip = patches.Ellipse((0, .0),
-                                        width=es,
-                                        height=es,
+                                        width=self.ellipse_size,
+                                        height=self.ellipse_size,
                                         angle=0)
             ref_ellip.set_facecolor((0, 0, 0))
-            ref_ax_loc = list(lpax2.get_position().bounds)
+            ref_ax_loc = list(self.ax_cb.get_position().bounds)
             ref_ax_loc[0] *= .95
             ref_ax_loc[1] -= .17
             ref_ax_loc[2] = .1
             ref_ax_loc[3] = .1
-            self.ref_ax = lpfig.add_axes(ref_ax_loc, aspect='equal')
+            self.ref_ax = self.fig.add_axes(ref_ax_loc, aspect='equal')
             self.ref_ax.add_artist(ref_ellip)
-            self.ref_ax.set_xlim(-es / 2. * 1.05, es / 2. * 1.05)
-            self.ref_ax.set_ylim(-es / 2. * 1.05, es / 2. * 1.05)
+            self.ref_ax.set_xlim(-self.ellipse_size / 2. * 1.05, self.ellipse_size / 2. * 1.05)
+            self.ref_ax.set_ylim(-self.ellipse_size / 2. * 1.05, self.ellipse_size / 2. * 1.05)
             plt.setp(self.ref_ax.xaxis.get_ticklabels(), visible=False)
             plt.setp(self.ref_ax.yaxis.get_ticklabels(), visible=False)
             self.ref_ax.set_title(r'$\Phi$ = 1')
@@ -1268,7 +1223,7 @@ def redraw_plot(self):
         use this function if you updated some attributes and want to re-plot.
         """
 
-        plt.close(self.fig)
+        self.fig.clf()
         self.plot()
 
     def export_params_to_file(self, save_path=None):
@@ -1339,7 +1294,7 @@ def export_params_to_file(self, save_path=None):
         tiazmap = np.zeros((xlist.shape[0], ylist.shape[0]))
         stationmap = np.zeros((xlist.shape[0], ylist.shape[0]), dtype='|S8')
 
-        station_location = {}  # a dict to store all MT stations (lan lat)
+        station_location = {}  # a dict to store all mt_obj stations (lan lat)
 
         # put the information into the zeroed arrays
         for ii in range(nx):
diff --git a/mtpy/imaging/phase_tensor_pseudosection.py b/mtpy/imaging/phase_tensor_pseudosection.py
index 7e4ccae7a..9221aad9d 100644
--- a/mtpy/imaging/phase_tensor_pseudosection.py
+++ b/mtpy/imaging/phase_tensor_pseudosection.py
@@ -112,16 +112,9 @@ class PlotPhaseTensorPseudoSection(mtpl.PlotSettings):
                         start and stop of station name indicies.
                         ex: for MT01dr station_id=(0,4) will be MT01
 
-        **rotz** : float or np.ndarray
+        **rotation_angle** : float or np.ndarray
                    angle in degrees to rotate the data clockwise positive.
-                   Can be an array of angle to individually rotate stations or
-                   periods or both.
-                       - If rotating each station by a constant
-                         angle the array needs to have a shape of
-                         (# of stations)
-                        - If rotating by period needs to have shape
-                           # of periods
-                        - If rotating both individually shape=(ns, nf)
+                   a single number for all input files.
                   *Default* is 0
 
         **title** : string
@@ -131,7 +124,7 @@ class PlotPhaseTensorPseudoSection(mtpl.PlotSettings):
                   dots per inch of the resolution. *default* is 300
 
 
-        **fignum** : int
+        **fig_num** : int
                      figure number.  *Default* is 1
 
         **plot_tipper** : [ 'yri' | 'yr' | 'yi' | 'n' ]
@@ -394,19 +387,6 @@ def __init__(self, **kwargs):
         if 'text_offset_y' not in list(self.scale_arrow_dict.keys()):
             self.scale_arrow_dict['text_offset_y'] = 0.
 
-        self._rot_z = kwargs.pop('rot_z', 0)
-        if isinstance(self._rot_z, float) or isinstance(self._rot_z, int):
-            self._rot_z = np.array([self._rot_z] * len(self.mt_list))
-
-        # if the rotation angle is an array for rotation of different
-        # freq than repeat that rotation array to the len(mt_list)
-        elif isinstance(self._rot_z, np.ndarray):
-            if self._rot_z.shape[0] != len(self.mt_list):
-                self._rot_z = np.repeat(self._rot_z, len(self.mt_list))
-
-        else:
-            pass
-
         # --> set induction arrow properties -------------------------------
         self.plot_tipper = kwargs.pop('plot_tipper', 'n')
 
@@ -434,9 +414,6 @@ def rotation_angle(self, value):
         only a single value is allowed
         """
         for ii, mt in enumerate(self.mt_list):
-            # JP: need to set the rotation angle negative for plotting
-            # I think its because the way polar plots work by measuring 
-            # counter clockwise
             mt.rotation_angle = value
             
         self._rotation_angle = value
@@ -458,6 +435,7 @@ def plot(self, show=True):
 
         # create a plot instance
         self.fig = plt.figure(self.fig_num, self.fig_size, dpi=self.fig_dpi)
+        self.fig.clf()
         self.ax = self.fig.add_subplot(1, 1, 1, aspect='equal')
 
         # FZ: control tick rotation=30 not that good
@@ -532,7 +510,7 @@ def plot(self, show=True):
             azimuth = pt.azimuth[::-1]
 
             # if there are induction arrows, flip them as pt
-            if self.plot_tipper.find('y') == 0:
+            if 'y' in self.plot_tipper:
                 tip = mt.Tipper
                 if tip.mag_real is not None:
                     tmr = tip.mag_real[::-1]
@@ -551,26 +529,27 @@ def plot(self, show=True):
 
             # get the properties to color the ellipses by
             if self.ellipse_colorby == 'phimin':
-                colorarray = pt.phimin[::-1]
+                color_array = pt.phimin[::-1]
 
             elif self.ellipse_colorby == 'phimax':
-                colorarray = pt.phimin[::-1]
+                color_array = pt.phimin[::-1]
 
             elif self.ellipse_colorby == 'phidet':
-                colorarray = np.sqrt(abs(pt.det[::-1])) * (180 / np.pi)
+                color_array = np.sqrt(abs(pt.det[::-1])) * (180 / np.pi)
 
             elif self.ellipse_colorby == 'skew' or \
                     self.ellipse_colorby == 'skew_seg':
-                colorarray = pt.beta[::-1]
+                color_array = pt.beta[::-1]
 
             elif self.ellipse_colorby == 'normalized_skew' or \
                     self.ellipse_colorby == 'normalized_skew_seg':
-                colorarray = 2 * pt.beta[::-1]
+                color_array = 2 * pt.beta[::-1]
 
             elif self.ellipse_colorby == 'ellipticity':
-                colorarray = pt.ellipticity[::-1]
+                color_array = pt.ellipticity[::-1]
             elif self.ellipse_colorby in ['strike', 'azimuth']:
-                colorarray = pt.azimuth[::-1] % 180
+                color_array = pt.azimuth[::-1] % 180
+                color_array[color_array > 90] -= 180                
             else:
                 raise NameError(self.ellipse_colorby + ' is not supported')
 
@@ -585,8 +564,8 @@ def plot(self, show=True):
                     plot_periodlist = periodlist
 
             # get min and max of the color array for scaling later
-            minlist.append(min(colorarray))
-            maxlist.append(max(colorarray))
+            minlist.append(min(color_array))
+            maxlist.append(max(color_array))
 
             for jj, ff in enumerate(periodlist):
 
@@ -596,25 +575,28 @@ def plot(self, show=True):
 
                 # create an ellipse scaled by phimin and phimax and orient
                 # the ellipse so that north is up and east is right
-                # need to add 90 to do so instead of subtracting
+                # Ellipse patch assumes the angle measure counterclockwise
+                # therefore need to multiply the azimuth by -1 because
+                # it is measured clockwise positive and subtract 90 to put
+                # it in a coordinate system of N=0, E=90
                 ellipd = patches.Ellipse((offset * self.xstretch,
                                           np.log10(ff) * self.ystretch),
                                          width=ewidth,
                                          height=eheight,
                                          edgecolor='k',
                                          lw=0.5,
-                                         angle=azimuth[jj] + 90)
+                                         angle=90 - azimuth[jj])
 
                 # get ellipse color
                 if cmap.find('seg') > 0:
-                    ellipd.set_facecolor(mtcl.get_plot_color(colorarray[jj],
+                    ellipd.set_facecolor(mtcl.get_plot_color(color_array[jj],
                                                              self.ellipse_colorby,
                                                              cmap,
                                                              ckmin,
                                                              ckmax,
                                                              bounds=bounds))
                 else:
-                    ellipd.set_facecolor(mtcl.get_plot_color(colorarray[jj],
+                    ellipd.set_facecolor(mtcl.get_plot_color(color_array[jj],
                                                              self.ellipse_colorby,
                                                              cmap,
                                                              ckmin,
@@ -683,10 +665,6 @@ def plot(self, show=True):
         self._plot_periodlist = plot_periodlist
         n = len(plot_periodlist)
 
-        # calculate minimum period and maximum period with a stretch factor
-        #        pmin = np.log10(plot_periodlist.min())*self.ystretch
-        #        pmax = np.log10(plot_periodlist.max())*self.ystretch
-
         pmin = int(np.floor(np.log10(plot_periodlist.min())))
         pmax = int(np.ceil(np.log10(plot_periodlist.max())))
 
@@ -698,28 +676,17 @@ def plot(self, show=True):
 
         self.offsetlist = offset_sort['offset']
         self.stationlist = offset_sort['station']
-        #        if self.offsetlist[0] > 0:
-        #            print 'rotating'
-        #            print self.stationlist
-        #            self.stationlist = self.stationlist[::-1]
 
         # set y-ticklabels
         if self.tscale == 'period':
             yticklabels = [mtpl.labeldict[ii]
                            for ii in range(pmin, pmax + 1, 1)]
-            #            yticklabels = ['{0:>4}'.format('{0: .1e}'.format(plot_period_list[ll]))
-            #                            for ll in np.arange(0, n, self.ystep)]+\
-            #                        ['{0:>4}'.format('{0: .1e}'.format(plot_period_list[-1]))]
 
             self.ax.set_ylabel('Period (s)',
                                fontsize=self.font_size + 2,
                                fontweight='bold')
 
         elif self.tscale == 'frequency':
-            #            yticklabels = ['{0:>4}'.format('{0: .1e}'.format(1./plot_period_list[ll]))
-            #                            for ll in np.arange(0, n, self.ystep)]+\
-            #                            ['{0:>4}'.format('{0: .1e}'.format(1./plot_period_list[-1]))]
-            #
             yticklabels = [mtpl.labeldict[-ii]
                            for ii in range(pmin, pmax + 1, 1)]
             self.ax.set_ylabel('Frequency (Hz)',
@@ -732,15 +699,10 @@ def plot(self, show=True):
 
         # --> set tick locations and labels
         # set y-axis major ticks
-        #        self.ax.yaxis.set_ticks([np.log10(plot_periodlist[ll])*self.ystretch
-        #                             for ll in np.arange(0, n, self.ystep)]):
         self.ax.yaxis.set_ticks(np.arange(pmin * self.ystretch,
                                           (pmax + 1) * self.ystretch,
                                           self.ystretch))
 
-        # set y-axis minor ticks
-        #        self.ax.yaxis.set_ticks([np.log10(plot_periodlist[ll])*self.ystretch
-        #                             for ll in np.arange(0, n, 1)],minor=True)
         # set y-axis tick labels
         self.ax.set_yticklabels(yticklabels)
 
@@ -1199,45 +1161,45 @@ def writeTextFiles(self, save_path=None, ptol=0.10):
             tpiazlines.append(tpiazline + '\n')
 
         # write files
-        skfid = file(os.path.join(svpath, 'PseudoSection.skew'), 'w')
+        skfid = open(os.path.join(svpath, 'PseudoSection.skew'), 'w')
         skfid.writelines(sklines)
         skfid.close()
 
-        phiminfid = file(os.path.join(svpath, 'PseudoSection.phimin'), 'w')
+        phiminfid = open(os.path.join(svpath, 'PseudoSection.phimin'), 'w')
         phiminfid.writelines(phiminlines)
         phiminfid.close()
 
-        phimaxfid = file(os.path.join(svpath, 'PseudoSection.phimax'),
+        phimaxfid = open(os.path.join(svpath, 'PseudoSection.phimax'),
                          'w')
         phimaxfid.writelines(phimaxlines)
         phimaxfid.close()
 
-        ellipfid = file(os.path.join(svpath, 'PseudoSection.ellipticity'),
+        ellipfid = open(os.path.join(svpath, 'PseudoSection.ellipticity'),
                         'w')
         ellipfid.writelines(elliplines)
         ellipfid.close()
 
-        azfid = file(os.path.join(svpath, 'PseudoSection.azimuth'),
+        azfid = open(os.path.join(svpath, 'PseudoSection.azimuth'),
                      'w')
         azfid.writelines(azimlines)
         azfid.close()
 
-        tprfid = file(os.path.join(svpath, 'PseudoSection.tipper_mag_real'),
+        tprfid = open(os.path.join(svpath, 'PseudoSection.tipper_mag_real'),
                       'w')
         tprfid.writelines(tprlines)
         tprfid.close()
 
-        tprazfid = file(os.path.join(svpath, 'PseudoSection.tipper_ang_real'),
+        tprazfid = open(os.path.join(svpath, 'PseudoSection.tipper_ang_real'),
                         'w')
         tprazfid.writelines(tprazlines)
         tprazfid.close()
 
-        tpifid = file(os.path.join(svpath, 'PseudoSection.tipper_mag_imag'),
+        tpifid = open(os.path.join(svpath, 'PseudoSection.tipper_mag_imag'),
                       'w')
         tpifid.writelines(tpilines)
         tpifid.close()
 
-        tpiazfid = file(os.path.join(svpath, 'PseudoSection.tipper_ang_imag'),
+        tpiazfid = open(os.path.join(svpath, 'PseudoSection.tipper_ang_imag'),
                         'w')
         tpiazfid.writelines(tpiazlines)
         tpiazfid.close()
diff --git a/mtpy/imaging/plot_mt_response.py b/mtpy/imaging/plot_mt_response.py
index 53dc99c3a..f28e0df7c 100644
--- a/mtpy/imaging/plot_mt_response.py
+++ b/mtpy/imaging/plot_mt_response.py
@@ -476,7 +476,9 @@ def plot(self, show=True, overlay_mt_obj=None):
                   'ellipticity': r'Ellipticity',
                   'skew_seg': r'Skew (deg)',
                   'normalized_skew_seg': r'Normalized Skew (deg)',
-                  'geometric_mean': r'$\sqrt{\Phi_{min} \cdot \Phi_{max}}$'}
+                  'geometric_mean': r'$\sqrt{\Phi_{min} \cdot \Phi_{max}}$',
+                  'strike': r"Strike (deg)",
+                  'azimuth': r"Strike (deg)"}
 
         if self.plot_tipper.find('y') == 0:
             if self.Tipper is None or np.all(self.Tipper.tipper == 0 + 0j):
@@ -542,6 +544,7 @@ def plot(self, show=True, overlay_mt_obj=None):
 
         # make figure instance
         self.fig = plt.figure(self.fig_num, self.fig_size, dpi=self.fig_dpi)
+        self.fig.clf()
 
         # --> make figure for xy,yx components
         if self.plot_num == 1 or self.plot_num == 3:
@@ -762,15 +765,15 @@ def plot(self, show=True, overlay_mt_obj=None):
         # -----plot tipper----------------------------------------------------
         if self.plot_tipper.find('y') == 0:
 
-            txr = self.Tipper.mag_real * np.sin(self.Tipper.angle_real * np.pi / 180 + \
-                                                np.pi * self.arrow_direction)
-            tyr = self.Tipper.mag_real * np.cos(self.Tipper.angle_real * np.pi / 180 + \
-                                                np.pi * self.arrow_direction)
+            txr = self.Tipper.mag_real * \
+                    np.cos(np.deg2rad(self.Tipper.angle_real))
+            tyr = self.Tipper.mag_real * \
+                    np.sin(np.deg2rad(self.Tipper.angle_real))
 
-            txi = self.Tipper.mag_imag * np.sin(self.Tipper.angle_imag * np.pi / 180 + \
-                                                np.pi * self.arrow_direction)
-            tyi = self.Tipper.mag_imag * np.cos(self.Tipper.angle_imag * np.pi / 180 + \
-                                                np.pi * self.arrow_direction)
+            txi = self.Tipper.mag_imag * \
+                    np.cos(np.deg2rad(self.Tipper.angle_imag))
+            tyi = self.Tipper.mag_imag * \
+                    np.sin(np.deg2rad(self.Tipper.angle_imag))
 
             nt = len(txr)
 
@@ -895,23 +898,27 @@ def plot(self, show=True, overlay_mt_obj=None):
             # get the properties to color the ellipses by
             if self.ellipse_colorby == 'phiminang' or \
                     self.ellipse_colorby == 'phimin':
-                colorarray = self.pt.phimin
+                color_array = self.pt.phimin
 
             elif self.ellipse_colorby == 'phimaxang' or \
                     self.ellipse_colorby == 'phimax':
-                colorarray = self.pt.phimax
+                color_array = self.pt.phimax
 
 
             elif self.ellipse_colorby == 'phidet':
-                colorarray = np.sqrt(abs(self.pt.det)) * (180 / np.pi)
+                color_array = np.sqrt(abs(self.pt.det)) * (180 / np.pi)
 
 
             elif self.ellipse_colorby == 'skew' or \
                     self.ellipse_colorby == 'skew_seg':
-                colorarray = self.pt.beta
+                color_array = self.pt.beta
 
             elif self.ellipse_colorby == 'ellipticity':
-                colorarray = self.pt.ellipticity
+                color_array = self.pt.ellipticity
+                
+            elif self.ellipse_colorby in ['strike', 'azimuth']:
+                color_array = self.pt.azimuth % 180
+                color_array[np.where(color_array > 90)] -= 180
 
             else:
                 raise NameError(self.ellipse_colorby + ' is not supported')
@@ -937,14 +944,14 @@ def plot(self, show=True, overlay_mt_obj=None):
 
                 # get ellipse color
                 if cmap.find('seg') > 0:
-                    ellipd.set_facecolor(mtcl.get_plot_color(colorarray[ii],
+                    ellipd.set_facecolor(mtcl.get_plot_color(color_array[ii],
                                                              self.ellipse_colorby,
                                                              cmap,
                                                              ckmin,
                                                              ckmax,
                                                              bounds=bounds))
                 else:
-                    ellipd.set_facecolor(mtcl.get_plot_color(colorarray[ii],
+                    ellipd.set_facecolor(mtcl.get_plot_color(color_array[ii],
                                                              self.ellipse_colorby,
                                                              cmap,
                                                              ckmin,
@@ -1364,7 +1371,7 @@ def redraw_plot(self):
             >>> p1.redraw_plot()
         """
 
-        plt.close(self.fig)
+        self.fig.clf()
         self.plot()
 
     def __str__(self):
diff --git a/mtpy/imaging/plotnresponses.py b/mtpy/imaging/plotnresponses.py
index 98319bdc0..d456da4db 100644
--- a/mtpy/imaging/plotnresponses.py
+++ b/mtpy/imaging/plotnresponses.py
@@ -309,23 +309,6 @@ def __init__(self, **kwargs):
         self.plot_style = kwargs.pop('plot_style', '1')
         self.plot_title = kwargs.pop('plot_title', None)
 
-        # if rotation angle is an int or float make an array the length of
-        # mt_list for plotting purposes
-        self._rot_z = kwargs.pop('rot_z', 0)
-        if isinstance(self._rot_z, float) or isinstance(self._rot_z, int):
-            self._rot_z = np.array([self._rot_z] * len(self.mt_list))
-
-        # if the rotation angle is an array for rotation of different
-        # freq than repeat that rotation array to the len(mt_list)
-        elif isinstance(self._rot_z, np.ndarray):
-            if self._rot_z.shape[0] != len(self.mt_list):
-                self._rot_z = np.repeat(self._rot_z, len(self.mt_list))
-
-        else:
-            pass
-
-        self._set_rot_z(self._rot_z)
-
         # set plot limits
         self.xlimits = kwargs.pop('xlimits', None)
         self.res_limits = kwargs.pop('res_limits', None)
@@ -385,34 +368,20 @@ def __init__(self, **kwargs):
         if self.plot_yn == 'y':
             self.plot()
 
-    # ---rotate data on setting rot_z
-    def _set_rot_z(self, rot_z):
+    #---need to rotate data on setting rotz
+    @property
+    def rotation_angle(self):
+        return self._rotation_angle
+    
+    @rotation_angle.setter
+    def rotation_angle(self, value):
         """
-        need to rotate data when setting z
+        only a single value is allowed
         """
-
-        # if rotation angle is an int or float make an array the length of
-        # mt_list for plotting purposes
-        if isinstance(rot_z, float) or isinstance(rot_z, int):
-            self._rot_z += np.array([rot_z] * len(self.mt_list))
-
-        # if the rotation angle is an array for rotation of different
-        # freq than repeat that rotation array to the len(mt_list)
-        elif isinstance(rot_z, np.ndarray):
-            if rot_z.shape[0] != len(self.mt_list):
-                self._rot_z += np.repeat(rot_z, len(self.mt_list))
-
-        else:
-            pass
-
         for ii, mt in enumerate(self.mt_list):
-            mt.rot_z = self._rot_z[ii]
-
-    def _get_rot_z(self):
-        return self._rot_z
-
-    rot_z = property(fget=_get_rot_z, fset=_set_rot_z,
-                     doc="""rotation angle(s)""")
+            mt.rotation_angle = value
+            
+        self._rotation_angle = value
 
     # --> on setting plot_ make sure to update the order and list
     def _set_plot_tipper(self, plot_tipper):
@@ -501,86 +470,6 @@ def plot(self, show=True):
         # set height ratios of the subplots
         hr = [2, 1.5] + [1] * (len(list(pdict.keys())) - 2)
 
-        #        if self.plot_style == '1':
-        #            self.plotlist = []
-        #
-        #            #--> plot from edi's if given, don't need to rotate because
-        #            #    data has already been rotated by the funcion _set_rot_z
-        ##            if self.fig_size is None:
-        ##                self.fig_size = [6, 6]
-        #            for ii, mt in enumerate(self.mt_list, 1):
-        #                p1 = plotresponse(mt_object=mt,
-        #                                  fig_num=ii,
-        #                                  fig_size=self.fig_size,
-        #                                  plot_num=self.plot_num,
-        #                                  fig_dpi=self.fig_dpi,
-        #                                  plot_yn='n',
-        #                                  plot_tipper=self._plot_tipper,
-        #                                  plot_strike=self._plot_strike,
-        #                                  plot_skew=self._plot_skew,
-        #                                  plot_pt=self._plot_pt)
-        #
-        #                #make sure all the properties are set to match the users
-        #                #line style between points
-        #                p1.xy_ls = self.xy_ls
-        #                p1.yx_ls = self.yx_ls
-        #                p1.det_ls = self.det_ls
-        #
-        #                #outline color
-        #                p1.xy_color = self.xy_color
-        #                p1.yx_color = self.yx_color
-        #                p1.det_color = self.det_color
-        #
-        #                #face color
-        #                p1.xy_mfc = self.xy_mfc
-        #                p1.yx_mfc = self.yx_mfc
-        #                p1.det_mfc = self.det_mfc
-        #
-        #                #maker
-        #                p1.xy_marker = self.xy_marker
-        #                p1.yx_marker = self.yx_marker
-        #                p1.det_marker = self.det_marker
-        #
-        #                #size
-        #                p1.marker_size = 2
-        #
-        #                #set plot limits
-        #                p1.xlimits = self.xlimits
-        #                p1.res_limits = self.res_limits
-        #                p1.phase_limits = self.phase_limits
-        #
-        #                #set font parameters
-        #                p1.font_size = self.font_size
-        #
-        #                #set arrow properties
-        #                p1.arrow_lw = self.arrow_lw
-        #                p1.arrow_head_width = self.arrow_head_width
-        #                p1.arrow_head_length = self.arrow_head_length
-        #                p1.arrow_color_real = self.arrow_color_real
-        #                p1.arrow_color_imag = self.arrow_color_imag
-        #                p1.arrow_direction = self.arrow_direction
-        #                p1.tipper_limits = self.tipper_limits
-        #
-        #                #skew properties
-        #                p1.skew_color = self.skew_color
-        #                p1.skew_marker = self.skew_marker
-        #
-        #                #strike properties
-        #                p1.strike_inv_marker = self.strike_inv_marker
-        #                p1.strike_inv_color = self.strike_inv_color
-        #
-        #                p1.strike_pt_marker = self.strike_pt_marker
-        #                p1.strike_pt_color = self.strike_pt_color
-        #
-        #                p1.strike_tip_marker = self.strike_tip_marker
-        #                p1.strike_tip_color = self.strike_tip_color
-        #
-        #                #--> plot the apparent resistivity and phase
-        #                self.plotlist.append(p1)
-        #
-        #                p1.plot()
-        #
-
         # -----Plot All in one figure with each plot as a subfigure------------
         if self.plot_style == 'all':
 
@@ -802,28 +691,6 @@ def plot(self, show=True):
                     pymin = min(0, min([min(rp.phase_xy), min(rp.phase_yx)]))
                     pymax = max(89.9, max([max(rp.phase_xy), max(rp.phase_yx)]))
                     self.phase_limits = (pymin, pymax)
-                #     self.phase_limits = (pymin, pymax)
-                # else:
-                #     self.phase_limits = (min(self.phase_limits[0], pymin),
-                #                          max(self.phase_limits[1], pymax))
-                # if self.phase_limits is None:
-                #     if min(rp.phasexy) < 0 or min(rp.phase_yx) < 0:
-                #         pymin = min([min(rp.phase_xy),
-                #                      min(rp.phase_yx)])
-                #         if pymin > 0:
-                #             pymin = 0
-                #     else:
-                #         pymin = 0
-                #
-                #     if max(rp.phasexy) > 90 or max(rp.phase_yx) > 90:
-                #         pymax = min([max(rp.phase_xy),  # YG: should use max instead ??
-                #                      max(rp.phase_yx)])
-                #         if pymax < 91:
-                #             pymax = 89.9  # YG: why??
-                #     else:
-                #         pymax = 89.9
-                #
-                #     self.phase_limits = (pymin, pymax)
 
                 # --> set axes properties
                 if ii == 0:
@@ -1182,6 +1049,11 @@ def plot(self, show=True):
 
                     elif self.ellipse_colorby == 'ellipticity':
                         colorarray = pt.ellipticity
+                    elif self.ellipse_colorby in ['strike', 'azimuth']:
+                        colorarray = self.fold_strike(pt.azimuth)
+                        self.ellipse_range = (-90, 90)
+                        ckmin = self.ellipse_range[0]
+                        ckmax = self.ellipse_range[1]
 
                     else:
                         raise NameError(self.ellipse_colorby + ' is not supported')
@@ -1888,31 +1760,6 @@ def plot(self, show=True):
                 # ------plot strike angles----------------------------------------------
                 if self._plot_strike.find('y') == 0:
 
-                    #                    if self._plot_strike.find('i') > 0:
-                    #                        #strike from invariants
-                    #                        zinv = mt.Z.invariants
-                    #                        s1 = zinv.strike
-                    #
-                    #                        #fold angles so go from -90 to 90
-                    #                        s1[np.where(s1>90)] -= -180
-                    #                        s1[np.where(s1<-90)] += 180
-                    #
-                    #                        #plot strike with error bars
-                    #                        ps1 = self.axst.errorbar(mt.period,
-                    #                                                s1,
-                    #                                                marker=mxy[ii % len(mxy)],
-                    #                                                ms=self.marker_size,
-                    #                                                mfc=cst[ii],
-                    #                                                mec=cst[ii],
-                    #                                                mew=self.marker_lw,
-                    #                                                ls='none',
-                    #                                                yerr=zinv.strike_err,
-                    #                                                ecolor=cst[ii],
-                    #                                                capsize=self.marker_size,
-                    #                                                elinewidth=self.marker_lw)
-                    #
-                    #                        stlist.append(ps1[0])
-
                     if self._plot_strike.find('p') > 0:
                         # strike from phase tensor
                         s2 = mt.pt.azimuth
diff --git a/mtpy/imaging/plotpt.py b/mtpy/imaging/plotpt.py
index 7e7bee9ca..04fe4ade5 100755
--- a/mtpy/imaging/plotpt.py
+++ b/mtpy/imaging/plotpt.py
@@ -167,19 +167,25 @@ class PlotPhaseTensor(mtpl.MTEllipse):
     """
 
     def __init__(self, **kwargs):
+        super().__init__()
 
         fn = kwargs.pop('fn', None)
         z_object = kwargs.pop('z_object', None)
         mt_object = kwargs.pop('mt_object', None)
         pt_object = kwargs.pop('pt_object', None)
+        
+        self._rotation_angle = 0
 
         #--> get mt object
         if fn is not None:
             self._mt = mtpl.MTplot(fn=fn)
+            self.pt = self._mt.pt
         elif z_object is not None:
             self._mt = mtpl.MTplot(z_object=z_object)
+            self.pt = self._mt.pt
         elif mt_object is not None:
             self._mt = mt_object
+            self.pt = self._mt.pt
         elif pt_object is not None:
             self.pt = pt_object
             self._mt = mtpl.MTplot()
@@ -222,11 +228,9 @@ def __init__(self, **kwargs):
         # read ellipse dict
         ellipse_dict = kwargs.pop('ellipse_dict', None)
         if ellipse_dict is None:
-            self._ellipse_dict = {'size': .25}
-        else:
-            self._ellipse_dict = ellipse_dict
+            ellipse_dict = {'size': .25}
 
-        self._read_ellipse_dict()
+        self._read_ellipse_dict(ellipse_dict)
 
         self.ellipse_spacing = kwargs.pop('ellipse_spacing', 1)
 
@@ -234,6 +238,31 @@ def __init__(self, **kwargs):
 
         if self.plot_yn == 'y':
             self.plot()
+            
+        #---need to rotate data on setting rotz
+    @property
+    def rotation_angle(self):
+        return self._rotation_angle
+    
+    @rotation_angle.setter
+    def rotation_angle(self, value):
+        """
+        only a single value is allowed
+        """
+        
+        self.pt.rotate(value)
+        self._mt.rotation_angle = value
+            
+        self._rotation_angle = value
+
+    def fold_strike(self, strike):
+        """
+        
+        """
+        strike = strike % 180
+        strike[np.where(strike > 90)] -= 180
+        
+        return strike
 
     def plot(self):
         """
@@ -256,16 +285,6 @@ def plot(self):
         self.fig = plt.figure(self.fig_num, self.fig_size, dpi=self.fig_dpi)
         plt.clf()
 
-        # get phase tensor instance
-        try:
-            self.pt
-            self.pt.rotate(self.rot_z)
-        except AttributeError:
-            self.pt = self._mt.pt
-            self.pt.rotate(self.rot_z)
-            self.zinv = self._mt.Z.invariants
-            self.zinv.rotate(self.rot_z)
-
         cmap = self.ellipse_cmap
         ckmin = self.ellipse_range[0]
         ckmax = self.ellipse_range[1]
@@ -281,29 +300,30 @@ def plot(self):
         # get the properties to color the ellipses by
         if self.ellipse_colorby == 'phiminang' or \
            self.ellipse_colorby == 'phimin':
-            colorarray = self.pt.phimin[0]
+            colorarray = self.pt.phimin
 
         elif self.ellipse_colorby == 'phidet':
-            colorarray = np.sqrt(abs(self.pt.det[0])) * (180 / np.pi)
+            colorarray = np.sqrt(abs(self.pt.det)) * (180 / np.pi)
 
         elif self.ellipse_colorby == 'skew' or\
                 self.ellipse_colorby == 'skew_seg':
-            colorarray = self.pt.beta[0]
+            colorarray = self.pt.beta
 
         elif self.ellipse_colorby == 'ellipticity':
-            colorarray = self.pt.ellipticity[0]
+            colorarray = self.pt.ellipticity
+        elif self.ellipse_colorby in ['strike', 'azimuth']:
+                colorarray = self.fold_strike(self.pt.azimuth)
 
         else:
             raise NameError(self.ellipse_colorby + ' is not supported')
 
         #-------------plotPhaseTensor-----------------------------------
         self.ax1 = self.fig.add_subplot(3, 1, 1, aspect='equal')
-        self._mt._period = 1. / self._mt.freq
 
         for ii, ff in enumerate(self._mt.period):
             # make sure the ellipses will be visable
-            if self.pt.phimax[0][ii] != 0:
-                eheight = self.pt.phimin[0][ii] / self.pt.phimax[0][ii] *\
+            if self.pt.phimax[ii] != 0:
+                eheight = self.pt.phimin[ii] / self.pt.phimax[ii] *\
                     self.ellipse_size
 
                 ewidth = self.ellipse_size
@@ -313,22 +333,13 @@ def plot(self):
                 eheight = 0.01 * self.ellipse_size
                 ewidth = 0.01 * self.ellipse_size
 
-            # ewidth = self.pt.phimax[0][ii]/self.pt.phimax[0][ii]*\
-            #                                                  self.ellipse_size
-
-            # alternative scaling
-            # eheight = self.pt.phimin[0][ii]/max(np.abs(self.pt.phimax[0]))*\
-            #                                                   self.ellipse_size
-            # ewidth = self.pt.phimax[0][ii]/max(np.abs(self.pt.phimax[0]))*\
-            #                                                   self.ellipse_size
-
             # create an ellipse scaled by phimin and phimax and oriented along
             # the azimuth which is calculated as clockwise but needs to
             # be plotted counter-clockwise hence the negative sign.
             ellipd = patches.Ellipse((np.log10(ff) * self.ellipse_spacing, 0),
                                      width=ewidth,
                                      height=eheight,
-                                     angle=90 - self.pt.azimuth[0][ii])
+                                     angle=90 - self.pt.azimuth[ii])
 
             self.ax1.add_patch(ellipd)
 
@@ -352,12 +363,13 @@ def plot(self):
         xlimits = (np.floor(np.log10(self._mt.period[0])),
                    np.ceil(np.log10(self._mt.period[-1])))
 
-        self.ax1.set_xlim(xlimits)
+        self.ax1.set_xlim((xlimits[0] * self.ellipse_spacing,
+                           xlimits[1] * self.ellipse_spacing))
         tklabels = []
         xticks = []
         for tk in self.ax1.get_xticks():
             try:
-                tklabels.append(mtpl.labeldict[tk])
+                tklabels.append(mtpl.labeldict[tk/self.ellipse_spacing])
                 xticks.append(tk)
             except KeyError:
                 pass
@@ -372,7 +384,7 @@ def plot(self):
                       which='major',
                       color=(.25, .25, .25),
                       lw=.25)
-
+        self.ax1.set_axisbelow(True)
         plt.setp(self.ax1.get_yticklabels(), visible=False)
         # add colorbar for PT
         self.cbax = self.fig.add_axes(self.cb_position)
@@ -415,12 +427,8 @@ def plot(self):
 
         #---------------plotStrikeAngle-----------------------------------
         self.ax2 = self.fig.add_subplot(3, 2, 3)
-        az = self.pt.azimuth[0]
-        az_err = self.pt.azimuth[1]
-
-        # put the strike into a coordinate system that goes from -90 to 90
-        az[np.where(az > 90)] -= 180
-        az[np.where(az < -90)] += 180
+        az = self.fold_strike(self.pt.azimuth)
+        az_err = self.pt.azimuth_err
 
         stlist = []
         stlabel = []
@@ -442,12 +450,8 @@ def plot(self):
         stlist.append(ps2[0])
         stlabel.append('PT')
         try:
-            strike = self.zinv.strike
-            strikeerr = np.nan_to_num(self.zinv.strike_err)
-            # put the strike into a coordinate system that goes from -90 to 90
-            strike[np.where(strike > 90)] = strike[np.where(strike > 90)] - 180
-            strike[np.where(strike < -90)
-                   ] = strike[np.where(strike < -90)] + 180
+            strike = self.fold_strike(self.zinv.strike)
+            strike_err = np.nan_to_num(self.zinv.strike_err)
 
             # plot invariant strike
             erxy = self.ax2.errorbar(self._mt.period,
@@ -458,7 +462,7 @@ def plot(self):
                                      mec=self.strike_inv_color,
                                      mew=self.marker_lw,
                                      ls='none',
-                                     yerr=strikeerr,
+                                     yerr=strike_err,
                                      ecolor=self.strike_inv_color,
                                      capsize=self.marker_size,
                                      elinewidth=self.marker_lw)
@@ -468,14 +472,10 @@ def plot(self):
         except AttributeError:
             print('Could not get z_invariants from pt, input z if desired.')
 
-        if self._mt.tipper is not None:
+        if self._mt.Tipper.tipper is not None:
             # strike from tipper
             tp = self._mt.Tipper
-            s3 = tp.angle_real + 90
-
-            # fold to go from -90 to 90
-            s3[np.where(s3 > 90)] = s3[np.where(s3 > 90)] - 180
-            s3[np.where(s3 < -90)] = s3[np.where(s3 < -90)] + 180
+            s3 = self.fold_strike(tp.angle_real)
 
             # plot strike with error bars
             ps3 = self.ax2.errorbar(self._mt.period,
@@ -510,7 +510,7 @@ def plot(self):
         plt.setp(ltext, fontsize=6)    # the legend text fontsize
 
         if self.strike_limits is None:
-            self.strike_limits = (-89.99, 89.99)
+            self.strike_limits = (-94, 94)
 
         self.ax2.set_yscale('linear')
         self.ax2.set_xscale('log', nonposx='clip')
@@ -524,10 +524,10 @@ def plot(self):
         self.ax2.set_title('Strike', fontdict=font_dictt)
 
         #---------plot Min & Max Phase-----------------------------------------
-        minphi = self.pt.phimin[0]
-        minphierr = self.pt.phimin[1]
-        maxphi = self.pt.phimax[0]
-        maxphierr = self.pt.phimax[1]
+        minphi = self.pt.phimin
+        minphierr = self.pt.phimin_err
+        maxphi = self.pt.phimax
+        maxphierr = self.pt.phimax_err
 
         self.ax3 = self.fig.add_subplot(3, 2, 4, sharex=self.ax2)
 
@@ -558,10 +558,10 @@ def plot(self):
                                   elinewidth=self.marker_lw)
 
         if self.pt_limits is None:
-            self.pt_limits = [min([self.pt.phimax[0].min(),
-                                   self.pt.phimin[0].min()]) - 3,
-                              max([self.pt.phimax[0].max(),
-                                   self.pt.phimin[0].max()]) + 3]
+            self.pt_limits = [min([self.pt.phimax.min(),
+                                   self.pt.phimin.min()]) - 3,
+                              max([self.pt.phimax.max(),
+                                   self.pt.phimin.max()]) + 3]
             if self.pt_limits[0] < -10:
                 self.pt_limits[0] = -9.9
             if self.pt_limits[1] > 100:
@@ -595,8 +595,8 @@ def plot(self):
 
         #-----------------------plotSkew---------------------------------------
 
-        skew = self.pt.beta[0]
-        skewerr = self.pt.beta[1]
+        skew = self.pt.beta
+        skewerr = self.pt.beta_err
 
         self.ax4 = self.fig.add_subplot(3, 2, 5, sharex=self.ax2)
         erskew = self.ax4.errorbar(self._mt.period,
@@ -639,8 +639,8 @@ def plot(self):
         self.ax4.set_title('Skew Angle', fontdict=font_dictt)
 
         #----------------------plotEllipticity--------------------------------
-        ellipticity = self.pt.ellipticity[0]
-        ellipticityerr = self.pt.ellipticity[1]
+        ellipticity = self.pt.ellipticity
+        ellipticityerr = self.pt.ellipticity_err
 
         self.ax5 = self.fig.add_subplot(3, 2, 6, sharex=self.ax2)
         erskew = self.ax5.errorbar(self._mt.period,
diff --git a/mtpy/imaging/plotresponse.py b/mtpy/imaging/plotresponse.py
index 318101496..415c859dc 100644
--- a/mtpy/imaging/plotresponse.py
+++ b/mtpy/imaging/plotresponse.py
@@ -615,7 +615,7 @@ def plot(self):
         #phase_yx: Note add 180 to place it in same quadrant as phase_xy
         self.ebyxp = mtpl.plot_errorbar(self.axp, 
                                        self.mt.period, 
-                                       self.mt.Z.phase_yx+180*self.phase_quadrant, 
+                                       self.mt.Z.phase_yx + 180*self.phase_quadrant, 
                                        marker=self.yx_marker, 
                                        ms=self.marker_size,  
                                        color=self.yx_color, 
@@ -626,14 +626,14 @@ def plot(self):
 
         #check the phase to see if any point are outside of [0:90]
         if self.phase_limits == None:
-            if min(self.mt.Z.phase_xy)<0 or min(self.mt.Z.phase_yx)<0:
+            if min(self.mt.Z.phase_xy) < 0 or min(self.mt.Z.phase_yx) < 0:
                 pymin = min([min(self.mt.Z.phase_xy), min(self.mt.Z.phase_yx)])
                 if pymin > 0:
                    pymin = 0
             else:
                 pymin = 0
             
-            if max(self.mt.Z.phase_xy)>90 or max(self.mt.Z.phase_yx)>90:
+            if max(self.mt.Z.phase_xy) > 90 or max(self.mt.Z.phase_yx) > 90:
                 pymax = min([max(self.mt.Z.phase_xy), max(self.mt.Z.phase_yx)])
                 if pymax < 91:
                     pymax = 89.9
@@ -817,8 +817,8 @@ def plot(self):
                 s2 = self.mt.pt.azimuth
                 s2_err = self.mt.pt.azimuth_err
                 #fold angles to go from -90 to 90
-                s2[np.where(s2>90)] -= 180
-                s2[np.where(s2<-90)] += 180
+                s2[np.where(s2 > 90)] -= 180
+                s2[np.where(s2 < -90)] += 180
                 
                 #plot strike with error bars
                 ps2 = mtpl.plot_errorbar(self.axst, 
@@ -839,7 +839,7 @@ def plot(self):
             
             if self._plot_strike.find('t') > 0:
                 #strike from tipper
-                s3 = self.mt.Tipper.angle_real+90
+                s3 = self.mt.Tipper.angle_real + 90
                 
                 #fold to go from -90 to 90
                 s3[np.where(s3 > 90)] -= 180
diff --git a/mtpy/imaging/plotstrike.py b/mtpy/imaging/plotstrike.py
index cd8f5b242..28018b968 100644
--- a/mtpy/imaging/plotstrike.py
+++ b/mtpy/imaging/plotstrike.py
@@ -16,8 +16,6 @@
 import mtpy.imaging.mtplottools as mtpl
 
 #==============================================================================
-
-
 class PlotStrike(object):
     """
     PlotStrike will plot the strike estimated from the invariants, phase tensor
@@ -179,11 +177,12 @@ def __init__(self, **kwargs):
         self.color_tip = (.2, .65, .2)
         self.ring_spacing = 10
         self.ring_limits = None
-        self.plot_orthogonal = True
+        self.plot_orthogonal = False
         
         self.font_size = 7
         self.text_pad = 0.6
         self.text_size = self.font_size
+        self.polar_limits = (np.deg2rad(-180), np.deg2rad(180))
             
         for key, value in kwargs.items():
             setattr(self, key, value)
@@ -227,6 +226,12 @@ def rotation_angle(self, value):
     def make_strike_array(self):
         """
         make strike array
+        
+        ..note:: Polar plots assume the azimuth is an angle measured 
+                counterclockwise positive from x = 0.  Therefore all angles
+                are calculated as 90 - angle to make them conform to the
+                polar plot convention.
+                
         """
         inv_list = []
         pt_list = []
@@ -240,15 +245,12 @@ def make_strike_array(self):
             
             if mt.period.size > nt:
                 nt = mt.period.size
-            #-----------get strike angle from invariants-----------------------
+            #-----------get strike angle from invariants----------------------
             zinv = Zinvariants(mt.Z)
 
-            # add 90 degrees because invariants assume 0 is north, but plotting
-            # assumes that 90 is north and measures clockwise, thus the negative
-            # because the strike angle from invariants is measured
-            # counter-clockwise
-
-            zs = 90 - zinv.strike
+            # subtract 90 because polar plot assumes 0 is on the x an 90 is 
+            # on the y
+            zs = 90 - zinv.strike 
 
             # fold so the angle goes from 0 to 180
             if self.fold == True:
@@ -264,7 +266,9 @@ def make_strike_array(self):
             mdictinv = dict([(ff, jj) for ff, jj in zip(mt.period, zs)])
             inv_list.append(mdictinv)
 
-            #------------get strike from phase tensor strike angle-------------
+            #------------get strike from phase tensor strike angle------------
+            # subtract 90 because polar plot assumes 0 is on the x an 90 is 
+            # on the y
             pt = mt.pt
             az = 90 - pt.azimuth
             az_err = pt.azimuth_err
@@ -296,10 +300,12 @@ def make_strike_array(self):
                                               dtype='complex')
                 tip.compute_components()
 
-            # needs to be negative because measures clockwise
-            tipr = -tip.angle_real
+            # # subtract 90 because polar plot assumes 0 is on the x an 90 is 
+            # on the y
+            tipr = 90 - tip.angle_real
 
-            tipr[np.where(tipr == 180.)] = 0.0
+            tipr[np.where(abs(tipr) == 180.)] = np.nan
+            tipr[np.where(abs(tipr) == 0)] = np.nan
 
             # fold so the angle goes from 0 to 180
             if self.fold == True:
@@ -316,8 +322,10 @@ def make_strike_array(self):
             tip_list.append(tiprdict)
 
         #--> get min and max period
-        self.max_per = np.amax([np.max(list(mm.keys())) for mm in inv_list], axis=0)
-        self.min_per = np.amin([np.min(list(mm.keys())) for mm in pt_list], axis=0)
+        self.max_per = np.amax([np.max(list(mm.keys())) for mm in inv_list],
+                               axis=0)
+        self.min_per = np.amin([np.min(list(mm.keys())) for mm in pt_list],
+                               axis=0)
 
         # make empty arrays to put data into for easy manipulation
         medinv = np.zeros((nt, nc))
@@ -350,8 +358,7 @@ def make_strike_array(self):
         self.med_inv = medinv
         self.med_pt = medpt
         self.med_tip = medtipr
-        
-        
+             
     def get_mean(self, st_array):
         """
         get mean value
@@ -476,9 +483,11 @@ def plot(self, show=True):
                     self.ax_pt = self.fig.add_subplot(n_subplots, nb, jj + nb,
                                                        polar=True)
                 elif 'v' in self.plot_orientation:
-                    self.ax_inv = self.fig.add_subplot(nb, n_subplots, 2*jj-1,
+                    self.ax_inv = self.fig.add_subplot(nb, n_subplots,
+                                                       jj * n_subplots - 2,
                                                        polar=True)
-                    self.ax_pt = self.fig.add_subplot(nb, n_subplots, 2*jj,
+                    self.ax_pt = self.fig.add_subplot(nb, n_subplots, 
+                                                      jj * n_subplots - 1,
                                                       polar=True)
                 ax_list = [self.ax_inv, self.ax_pt]
                 # vertical orientation
@@ -488,8 +497,8 @@ def plot(self, show=True):
                                                            jj + 2 * nb,
                                                            polar=True)
                     elif 'v' in self.plot_orientation:
-                        self.ax_tip = self.fig.add_subplot(n_subplots, nb,
-                                                           2*jj + 2,
+                        self.ax_tip = self.fig.add_subplot(nb, n_subplots, 
+                                                           jj * n_subplots,
                                                            polar=True)
                     ax_list.append(self.ax_tip)
 
@@ -587,13 +596,14 @@ def plot(self, show=True):
 
                     # make a light grid
                     axh.grid(alpha=.25, zorder=0)
+                    axh.set_xlim(self.polar_limits)
 
                     # properties for the invariants
                     if aa == 0:
                         # limits need to be rotate 90 counter clockwise because
                         # we already rotated by 90 degrees so the range is
                         # from -90 to 270 with -90 being east
-                        axh.set_xlim(-90 * np.pi / 180, 270 * np.pi / 180)
+                        
 
                         # label the plot with the mode value of strike
                         # need to subtract 90 again because the histogram is
@@ -615,13 +625,15 @@ def plot(self, show=True):
                         #--> set title of subplot
                         if 'h' in self.plot_orientation:
                             axh.set_title(self.title_dict[bb], fontdict=fd,
-                                          bbox={'facecolor': 'white', 'alpha': .25})
+                                          bbox={'facecolor': 'white',
+                                                'alpha': .25})
 
                             #--> set the title offset
                             axh.titleOffsetTrans._t = (0, .1)
                         elif 'v' in self.plot_orientation:
                             axh.set_ylabel(self.title_dict[bb], fontdict=fd,
-                                           bbox={'facecolor': 'white', 'alpha': .25}, 
+                                           bbox={'facecolor': 'white',
+                                                 'alpha': .25}, 
                                            rotation=0,
                                            labelpad=50)
                             axh.yaxis.set_label_position("right")
@@ -630,7 +642,7 @@ def plot(self, show=True):
                     elif aa == 1:
                         # limits go from -180 to 180 as that is how the angle
                         # is calculated
-                        axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180)
+                        #axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180)
                         
                         pt_median, pt_mode, pt_mean = self.get_stats(plot_pt,
                                                                      pt_hist,
@@ -648,7 +660,7 @@ def plot(self, show=True):
                     # set tipper axes properties
                     elif aa == 2:
                         # limits go from -180 to 180
-                        axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180)
+                        #axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180)
 
                         tr_median, tr_mode, tr_mean = self.get_stats(tr,
                                                                      tr_hist,
@@ -1279,9 +1291,9 @@ def writeTextFiles(self, save_path=None):
             slisttip[kk + 2].append(tpmed)
             slisttip[kk + 3].append(tpmode)
 
-        invfid = file(os.path.join(svpath, 'Strike.invariants'), 'w')
-        ptfid = file(os.path.join(svpath, 'Strike.pt'), 'w')
-        tpfid = file(os.path.join(svpath, 'Strike.tipper'), 'w')
+        invfid = open(os.path.join(svpath, 'Strike.invariants'), 'w')
+        ptfid = open(os.path.join(svpath, 'Strike.pt'), 'w')
+        tpfid = open(os.path.join(svpath, 'Strike.tipper'), 'w')
 
         #---> write strike from the invariants
         # == > mean
diff --git a/mtpy/utils/calculator.py b/mtpy/utils/calculator.py
index 585bdaddb..b59bf7c27 100644
--- a/mtpy/utils/calculator.py
+++ b/mtpy/utils/calculator.py
@@ -42,13 +42,14 @@ def roundsf(number, sf):
     round a number to a specified number of significant figures (sf)
     """
     # can't have < 1 s.f.
-    sf = max(sf,1.)
+    sf = max(sf, 1.)
     rounding = int(np.ceil(-np.log10(number) + sf - 1.))
     
     return np.round(number, rounding)
 
 
-def get_period_list(period_min,period_max,periods_per_decade,include_outside_range=True):
+def get_period_list(period_min, period_max, periods_per_decade,
+                    include_outside_range=True):
     """
     get a list of values (e.g. periods), evenly spaced in log space and 
     including values on multiples of 10
@@ -118,7 +119,8 @@ def nearest_index(val,array):
     return np.where(diff==min(diff))[0][0]
 
 
-def make_log_increasing_array(z1_layer, target_depth, n_layers, increment_factor=0.999):
+def make_log_increasing_array(z1_layer, target_depth, n_layers,
+                              increment_factor=0.999):
     """
     create depth array with log increasing cells, down to target depth,
     inputs are z1_layer thickness, target depth, number of layers (n_layers)
@@ -189,17 +191,17 @@ def invertmatrix_incl_errors(inmatrix, inmatrix_err=None):
                     for l in range(2):
                         #each entry has 4 summands
 
-                        err += np.abs (- inv_matrix[i,k]  * inv_matrix[l,j]  * inmatrix_err[k,l])
-
+                        err += np.abs(- inv_matrix[i,k] * inv_matrix[l,j] *\
+                                      inmatrix_err[k,l])
 
                 inv_matrix_err[i,j] = err
 
- 
     return inv_matrix, inv_matrix_err
 
 def rhophi2z(rho, phi, freq):
     """
-        Convert impedance-style information given in Rho/Phi format into complex valued Z.
+        Convert impedance-style information given in Rho/Phi format into 
+        complex valued Z.
 
         Input:
         rho - 2x2 array (real) - in Ohm m
@@ -228,7 +230,8 @@ def rhophi2z(rho, phi, freq):
     return z 
 
 
-def compute_determinant_error(z_array, z_err_array, method = 'theoretical', repeats=1000):
+def compute_determinant_error(z_array, z_err_array, method='theoretical',
+                              repeats=1000):
     """
     compute the error of the determinant of z using a stochastic method
     seed random z arrays with a normal distribution around the input array
@@ -257,33 +260,44 @@ def compute_determinant_error(z_array, z_err_array, method = 'theoretical', repe
         error = np.std(detlist,axis=0)
     
     else:
-        error = np.abs(z_err_array[:,0,0]*np.abs(z_array[:,1,1]) + z_err_array[:,1,1]*np.abs(z_array[:,0,0]) \
-             - z_err_array[:,0,1]*np.abs(z_array[:,1,0]) - z_err_array[:,1,0]*np.abs(z_array[:,0,1]))
+        error = np.abs(z_err_array[:, 0, 0] * np.abs(z_array[:, 1, 1]) +\
+                       z_err_array[:, 1, 1] * np.abs(z_array[:, 0, 0]) \
+                       - z_err_array[:, 0, 1]*np.abs(z_array[:, 1, 0]) - \
+                       z_err_array[:, 1, 0]*np.abs(z_array[:, 0, 1]))
     
     return error
 
-
-
-
 def propagate_error_polar2rect(r,r_error,phi, phi_error):
     """
-        Find error estimations for the transformation from polar to cartesian coordinates.
-
-        Uncertainties in polar representation define a section of an annulus. Find the 4 corners of this section and additionally the outer boundary point, which is defined by phi = phi0, rho = rho0 + sigma rho.
-        The cartesian "box" defining the uncertainties in x,y is the outer bound around the annulus section, defined by the four outermost points. So check the four corners as well as the outer boundary edge of the section to find the extrema in x znd y. These give you the sigma_x/y. 
+        Find error estimations for the transformation from polar to cartesian
+        coordinates.
+
+        Uncertainties in polar representation define a section of an annulus.
+        Find the 4 corners of this section and additionally the outer boundary
+        point, which is defined by phi = phi0, rho = rho0 + sigma rho.
+        The cartesian "box" defining the uncertainties in x,y is the outer
+        bound around the annulus section, defined by the four outermost
+        points. So check the four corners as well as the outer boundary edge
+        of the section to find the extrema in x znd y. These give you the
+        sigma_x/y. 
 
     """ 
 
-    corners = [ ( np.real(cmath.rect(r-r_error, phi-phi_error)), np.imag(cmath.rect(r-r_error, phi-phi_error))),\
-                 ( np.real(cmath.rect(r+r_error, phi-phi_error)), np.imag(cmath.rect(r+r_error, phi-phi_error))),\
-                 ( np.real(cmath.rect(r+r_error, phi+phi_error)), np.imag(cmath.rect(r+r_error, phi+phi_error))),\
-                 ( np.real(cmath.rect(r-r_error, phi+phi_error)), np.imag(cmath.rect(r-r_error, phi+phi_error))),\
-                 ( np.real(cmath.rect(r+r_error, phi)), np.imag(cmath.rect(r+r_error, phi))) ]
+    corners = [( np.real(cmath.rect(r-r_error, phi-phi_error)),
+                np.imag(cmath.rect(r-r_error, phi-phi_error))),
+               (np.real(cmath.rect(r+r_error, phi-phi_error)),
+                np.imag(cmath.rect(r+r_error, phi-phi_error))),
+               (np.real(cmath.rect(r+r_error, phi+phi_error)),
+                np.imag(cmath.rect(r+r_error, phi+phi_error))),
+               (np.real(cmath.rect(r-r_error, phi+phi_error)),
+                np.imag(cmath.rect(r-r_error, phi+phi_error))),
+               (np.real(cmath.rect(r+r_error, phi)),
+                np.imag(cmath.rect(r+r_error, phi))) ]
 
     lo_x = [i[0] for i in corners]
     lo_y = [i[1] for i in corners]
 
-    point =  (np.real(cmath.rect(r, phi)), np.imag(cmath.rect(r, phi)) )
+    point =  (np.real(cmath.rect(r, phi)), np.imag(cmath.rect(r, phi)))
     lo_xdiffs = [ abs(point[0] - i) for i in lo_x]
     lo_ydiffs = [ abs(point[1] - i) for i in lo_y]
     
@@ -292,25 +306,24 @@ def propagate_error_polar2rect(r,r_error,phi, phi_error):
 
     return xerr, yerr
 
-
-
-
 def propagate_error_rect2polar(x,x_error,y, y_error):
     
     # x_error, y_error define a  rectangular uncertainty box  
     
-    # rho error is the difference between the closest and furthest point of the box (w.r.t. the origin)
+    # rho error is the difference between the closest and furthest point 
+    # of the box (w.r.t. the origin)
     # approximation: just take corners and midpoint of edges 
-    lo_points = [ (x + x_error, y), (x - x_error, y), (x, y - y_error ), (x, y + y_error ),\
-                  (x - x_error, y - y_error) ,(x + x_error, y - y_error) ,(x + x_error, y + y_error) ,(x - x_error, y + y_error) ]
-
+    lo_points = [(x + x_error, y), (x - x_error, y),
+                 (x, y - y_error ), (x, y + y_error ),
+                 (x - x_error, y - y_error), (x + x_error, y - y_error),
+                 (x + x_error, y + y_error), (x - x_error, y + y_error)]
 
     #check, if origin is within the box:
     origin_in_box = False
     if x_error >= np.abs(x) and y_error >= np.abs(y):
         origin_in_box = True
 
-    lo_polar_points = [ cmath.polar(np.complex(*i)) for i in lo_points ]
+    lo_polar_points = [cmath.polar(np.complex(*i)) for i in lo_points]
 
     lo_rho = [i[0] for i in lo_polar_points ]
     lo_phi = [math.degrees(i[1])%360 for i in lo_polar_points ]
@@ -319,14 +332,13 @@ def propagate_error_rect2polar(x,x_error,y, y_error):
     phi_err = 0.5*(max(lo_phi) - min(lo_phi))
 
     if (270 < max(lo_phi) < 360 ) and (0 < min(lo_phi) < 90 ):
-        tmp1 = [ i for i in lo_phi if (0 < i < 90) ]
-        tmp4 = [ i for i in lo_phi if (270 < i < 360) ]
-        phi_err = 0.5*((max(tmp1) - min(tmp4))%360)
+        tmp1 = [i for i in lo_phi if (0 < i < 90)]
+        tmp4 = [i for i in lo_phi if (270 < i < 360)]
+        phi_err = 0.5 * ((max(tmp1) - min(tmp4)) % 360)
 
 
     if phi_err > 180:
-        #print phi_err,' -> ',(-phi_err)%360
-        phi_err = (-phi_err)%360
+        phi_err = (-phi_err) % 360
 
     if origin_in_box is True:
         #largest rho:
@@ -336,8 +348,6 @@ def propagate_error_rect2polar(x,x_error,y, y_error):
 
     return rho_err, phi_err
 
-
-
 def z_error2r_phi_error(z_real, z_imag, error):
     """
     Error estimation from rectangular to polar coordinates.
@@ -366,9 +376,9 @@ def z_error2r_phi_error(z_real, z_imag, error):
     
     res_rel_err = 2.*z_rel_err
     
-    #if the relative error of the amplitude is >=100% that means that the relative 
-    #error of the resistivity is 200% - that is then equivalent to an uncertainty 
-    #in the phase angle of 90 degrees:
+    # if the relative error of the amplitude is >=100% that means that the
+    # relative error of the resistivity is 200% - that is then equivalent
+    # to an uncertainty in the phase angle of 90 degrees
     if np.iterable(z_real):
         phi_err = np.degrees(np.arctan(z_rel_err))   
         phi_err[res_rel_err > 1.] = 90.
@@ -378,58 +388,50 @@ def z_error2r_phi_error(z_real, z_imag, error):
             phi_err = 90
         else:
             phi_err = np.degrees(np.arctan(z_rel_err))    
-    
-    
+
     return res_rel_err, phi_err
     
-    
 
 def old_z_error2r_phi_error(x,x_error,y, y_error):
     """
-        Error estimation from rect to polar, but with small variation needed for 
-        MT: the so called 'relative phase error' is NOT the relative phase error,
-        but the ABSOLUTE uncertainty in the angle that corresponds to the relative
-        error in the amplitude. 
-
-        So, here we calculate the transformation from rect to polar coordinates, 
-        esp. the absolute/length of the value. Then we find the uncertainty in 
-        this length and calculate the relative error of this. The relative error of
-        the resistivity will be double this value, because it's calculated by taking 
-        the square of this length.
-        
-        The relative uncertainty in length defines a circle around (x,y) 
-        (APPROXIMATION!). The uncertainty in phi is now the absolute of the 
-        angle beween the vector to (x,y) and the origin-vector tangential to the
-        circle.
-        BUT....since the phase angle uncertainty is interpreted with regard to 
-        the resistivity and not the Z-amplitude, we have to look at the square of
-        the length, i.e. the relative error in question has to be halfed to get
-        the correct relationship between resistivity and phase errors!!.
+    Error estimation from rect to polar, but with small variation needed for 
+    MT: the so called 'relative phase error' is NOT the relative phase error,
+    but the ABSOLUTE uncertainty in the angle that corresponds to the relative
+    error in the amplitude. 
+
+    So, here we calculate the transformation from rect to polar coordinates, 
+    esp. the absolute/length of the value. Then we find the uncertainty in 
+    this length and calculate the relative error of this. The relative error of
+    the resistivity will be double this value, because it's calculated by taking 
+    the square of this length.
+    
+    The relative uncertainty in length defines a circle around (x,y) 
+    (APPROXIMATION!). The uncertainty in phi is now the absolute of the 
+    angle beween the vector to (x,y) and the origin-vector tangential to the
+    circle.
+    BUT....since the phase angle uncertainty is interpreted with regard to 
+    the resistivity and not the Z-amplitude, we have to look at the square of
+    the length, i.e. the relative error in question has to be halfed to get
+    the correct relationship between resistivity and phase errors!!.
 
     """
 
     # x_error, y_error define a  rectangular uncertainty box  
     
-    # rho error is the difference between the closest and furthest point of the box (w.r.t. the origin)
-    # approximation: just take corners and midpoint of edges 
-    lo_points = [ (x + x_error, y), (x - x_error, y), (x, y - y_error ), 
-                    (x, y + y_error ), (x - x_error, y - y_error) ,
-                    (x + x_error, y - y_error) ,(x + x_error, y + y_error) ,
-                    (x - x_error, y + y_error) ]
+    # rho error is the difference between the closest and furthest point
+    # of the box (w.r.t. the origin) approximation: just take corners and
+    # midpoint of edges 
+    lo_points = [(x + x_error, y), (x - x_error, y),
+                 (x, y - y_error ), (x, y + y_error ), 
+                 (x - x_error, y - y_error), (x + x_error, y - y_error),
+                 (x + x_error, y + y_error), (x - x_error, y + y_error)]
 
+    lo_polar_points = [cmath.polar(np.complex(*i)) for i in lo_points]
 
-    #check, if origin is within the box:
-    origin_in_box = False
-    if x_error >= np.abs(x) and y_error >= np.abs(y):
-        origin_in_box = True
-
-    lo_polar_points = [ cmath.polar(np.complex(*i)) for i in lo_points ]
-
-    lo_rho = [i[0] for i in lo_polar_points ]
-    lo_phi = [math.degrees(i[1])%360 for i in lo_polar_points ]
+    lo_rho = [i[0] for i in lo_polar_points]
 
     #uncertainty in amplitude is defined by half the diameter of the box around x,y
-    rho_err = 0.5*(max(lo_rho) - min(lo_rho) )
+    rho_err = 0.5 * (max(lo_rho) - min(lo_rho) )
 
     rho = cmath.polar(np.complex(x,y))[0] 
     try:
@@ -443,14 +445,11 @@ def old_z_error2r_phi_error(x,x_error,y, y_error):
     if rel_error_rho > 1.:
         phi_err = 90
     else:
-        phi_err = np.degrees(np.arcsin( rel_error_rho))
+        phi_err = np.degrees(np.arcsin(rel_error_rho))
 
 
     return rho_err, phi_err
 
-
-
-
 #rotation:
 #1. rotation positive in clockwise direction
 #2. orientation of new X-axis X' given by rotation angle
@@ -470,106 +469,142 @@ def old_z_error2r_phi_error(x,x_error,y, y_error):
 # That is NOT the same as the rotated error matrix Zerr (although the result is similar)
 
 
-def rotatematrix_incl_errors(inmatrix, angle, inmatrix_err = None) :
-   
-    if inmatrix is None :
-        raise MTex.MTpyError_inputarguments('Matrix AND eror matrix must be defined')
+def rotate_matrix_with_errors(in_matrix, angle, error=None):
+    """
+    
+    Rotate a matrix including errors clockwise given an angle in degrees.
+    
+    :param in_matrix: A n x 2 x 2  matrix to rotate 
+    :type inmatrix: np.ndarray
+    
+    :param angle: Angle to rotate by assuming clockwise positive from
+                 0 = north
+    :type angle: float
+    
+    :param error: A n x 2 x 2 matrix of associated errors,
+                        defaults to None
+    :type error: np.ndarray, optional
+    
+    :raises MTex: If input array is incorrect
+    
+    :return: rotated matrix
+    :rtype: np.ndarray
+    
+    :return: rotated matrix errors
+    :rtype: np.ndarray
 
-    if (inmatrix_err is not None) and (inmatrix.shape != inmatrix_err.shape):
-        raise MTex.MTpyError_inputarguments('Matrix and err-matrix shapes do not match: %s - %s'%(str(inmatrix.shape), str(inmatrix_err.shape)))
+    """
+   
+    if in_matrix is None:
+        raise MTex.MTpyError_inputarguments('Matrix must be defined')
 
+    if (error is not None) and (in_matrix.shape != error.shape):
+        msg = 'matricies are not the same shape in_matrix={0}, err={1}'.format(
+               in_matrix.shape, error.shape)
+        raise MTex.MTpyError_inputarguments(msg)
 
     try:
-        degreeangle = angle % 360
-    except:
-        raise MTex.MTpyError_inputarguments('"Angle" must be a valid number (in degrees)')
-
-    phi = math.radians(degreeangle)
+        phi = np.deg2rad(float(angle) % 360)
+    except TypeError:
+        raise MTex.MTpyError_inputarguments('"Angle" must be a float')
 
     cphi = np.cos(phi)
     sphi = np.sin(phi)
 
-    # JP: Changed the rotation matrix to be formulated to rotate
-    # counter clockwise, I cannot find a good reason for this except that 
-    # when you plot the strike and phase tensors the look correct with this
-    # formulation.
-    rotmat = np.array([[ cphi, sphi], [-sphi, cphi]])
-    # rotmat = np.array([[ cphi, -sphi], [sphi, cphi]])
-    rotated_matrix = np.dot(np.dot(rotmat, inmatrix), np.linalg.inv(rotmat))
+    # clockwise rotation matrix is given by [[cos, -sin], [sin, cos]]
+    rot_mat = np.array([[ cphi, -sphi], [sphi, cphi]])
+    rotated_matrix = np.dot(np.dot(rot_mat, in_matrix), np.linalg.inv(rot_mat))
 
-    errmat  = None
-    if (inmatrix_err is not None) :
-        err_orig = np.real(inmatrix_err) 
-        errmat = np.zeros_like(inmatrix_err)
+    err_mat  = None
+    if (error is not None) :
+        err_orig = np.real(error) 
+        err_mat = np.zeros_like(error)
 
         # standard propagation of errors:
-        errmat[0,0] = np.sqrt( (cphi**2 * err_orig[0,0])**2 + \
-                              (cphi * sphi * err_orig[0,1])**2 + \
-                                  (cphi * sphi * err_orig[1,0])**2 + \
-                                      (sphi**2 * err_orig[1,1])**2 )
-        errmat[0,1] = np.sqrt( (cphi**2 * err_orig[0,1])**2 + \
-                              (cphi * sphi * err_orig[1,1])**2 + \
-                                  (cphi * sphi * err_orig[0,0])**2 + \
-                                      (sphi**2 * err_orig[1,0])**2 )
-        errmat[1,0] = np.sqrt( (cphi**2 * err_orig[1,0])**2 + \
-                              (cphi * sphi * err_orig[1,1])**2 +\
-                                  (cphi * sphi * err_orig[0,0])**2 + \
-                                      (sphi**2 * err_orig[0,1])**2 )
-        errmat[1,1] = np.sqrt( (cphi**2 * err_orig[1,1])**2 + \
-                              (cphi * sphi * err_orig[0,1])**2 + \
-                                  (cphi * sphi * err_orig[1,0])**2 + \
-                                      (sphi**2 * err_orig[0,0])**2 )
-
-    return rotated_matrix, errmat
-
-
-def rotatevector_incl_errors(invector, angle, invector_err = None):
-    #check for row or column vector 
-    
-    if invector is None :
-        raise MTex.MTpyError_inputarguments('Vector AND error-vector must be defined')
-
-    if (invector_err is not None) and (invector.shape != invector_err.shape):
-        raise MTex.MTpyError_inputarguments('Vector and errror-vector shapes do not match: %s - %s'%(str(invector.shape), str(invector_err.shape)))
+        err_mat[0,0] = np.sqrt((cphi**2 * err_orig[0, 0])**2 + \
+                               (cphi * sphi * err_orig[0, 1])**2 + \
+                               (cphi * sphi * err_orig[1, 0])**2 + \
+                               (sphi**2 * err_orig[1, 1])**2)
+        err_mat[0,1] = np.sqrt((cphi**2 * err_orig[0, 1])**2 + \
+                               (cphi * sphi * err_orig[1, 1])**2 + \
+                               (cphi * sphi * err_orig[0, 0])**2 + \
+                               (sphi**2 * err_orig[1, 0])**2)
+        err_mat[1,0] = np.sqrt((cphi**2 * err_orig[1, 0])**2 + \
+                               (cphi * sphi * err_orig[1, 1])**2 +\
+                               (cphi * sphi * err_orig[0, 0])**2 + \
+                               (sphi**2 * err_orig[0, 1])**2)
+        err_mat[1,1] = np.sqrt((cphi**2 * err_orig[1, 1])**2 + \
+                               (cphi * sphi * err_orig[0, 1])**2 + \
+                               (cphi * sphi * err_orig[1, 0])**2 + \
+                               (sphi**2 * err_orig[0, 0])**2)
+
+    return rotated_matrix, err_mat
+
+
+def rotate_vector_with_errors(in_vector, angle, error=None):
+    """
+    
+    Rotate a vector including errors clockwise given an angle in degrees.
+    
+    :param in_matrix: A n x 1 x 2  vector to rotate 
+    :type invector: np.ndarray
+    
+    :param angle: Angle to rotate by assuming clockwise positive from
+                 0 = north
+    :type angle: float
+    
+    :param error: A n x 1 x 2 vector of associated errors,
+                        defaults to None
+    :type error: np.ndarray, optional
+    
+    :raises MTex: If input array is incorrect
+    
+    :return: rotated vector
+    :rtype: np.ndarray
+    
+    :return: rotated vector errors
+    :rtype: np.ndarray
 
-    try:
-        degreeangle = angle%360
-    except:
-        raise MTex.MTpyError_inputarguments('"Angle" must be a valid number (in degrees)')
+    """
+    
+    if in_vector is None:
+        raise MTex.MTpyError_inputarguments('Vector AND error-vector must'+
+                                            ' be defined')
 
-    phi = math.radians(degreeangle)
+    if (error is not None) and (in_vector.shape != error.shape):
+        msg = 'matricies are not the same shape in_vector={0}, err={1}'.format(
+               in_vector.shape, error.shape)
+        raise MTex.MTpyError_inputarguments(msg)
     
+    try:
+        phi = np.deg2rad(float(angle) % 360)
+    except TypeError:
+        raise MTex.MTpyError_inputarguments('"Angle" must be a float')
+
     cphi = np.cos(phi)
     sphi = np.sin(phi)
-
-    # JP: Changed the rotation matrix to be formulated to rotate
-    # counter clockwise, I cannot find a good reason for this except that 
-    # when you plot the strike and phase tensors the look correct with this
-    # formulation.
-    rotmat = np.array([[ cphi, sphi],[-sphi, cphi]])
-    # rotmat = np.array([[ cphi, -sphi],[sphi, cphi]])
-
-    if invector.shape == (1, 2):
-        rotated_vector = np.dot( invector, np.linalg.inv(rotmat) )
-    else:
-        rotated_vector = np.dot( rotmat, invector )
     
-    
-    errvec = None
-    if (invector_err is not None) :   
-        errvec = np.zeros_like(invector_err)
+    rot_mat = np.array([[ cphi, -sphi],[sphi, cphi]])
 
-        if invector_err.shape == (1, 2):
-            errvec = np.dot(invector_err, np.abs(np.linalg.inv(rotmat)))
-        else:
-            errvec = np.dot(np.abs(rotmat), invector_err )
+    if in_vector.shape == (1, 2):
+        rotated_vector = np.dot(in_vector, np.linalg.inv(rot_mat) )
+    else:
+        rotated_vector = np.dot(rot_mat, in_vector)
 
+    err_vector = None
+    if (error is not None):   
+        err_vector = np.zeros_like(error)
 
-    return rotated_vector, errvec
+        if error.shape == (1, 2):
+            err_vector = np.dot(error, np.abs(np.linalg.inv(rot_mat)))
+        else:
+            err_vector = np.dot(np.abs(rot_mat), error)
 
+    return rotated_vector, err_vector
 
 
-def multiplymatrices_incl_errors(inmatrix1, inmatrix2, inmatrix1_err = None,inmatrix2_err = None ):
+def multiplymatrices_incl_errors(inmatrix1, inmatrix2, 
+                                 inmatrix1_err=None, inmatrix2_err=None):
 
     if inmatrix1 is None or inmatrix2 is None:
         raise MTex.MTpyError_inputarguments('ERROR - two 2x2 arrays needed as input')
@@ -601,20 +636,20 @@ def multiplymatrices_incl_errors(inmatrix1, inmatrix2, inmatrix1_err = None,inma
 
 def reorient_data2D(x_values, y_values, x_sensor_angle = 0 , y_sensor_angle = 90):
     """
-        Re-orient time series data of a sensor pair, which has not been in default (x=0, y=90) orientation.
+    Re-orient time series data of a sensor pair, which has not been in default (x=0, y=90) orientation.
 
-        Input:
-        - x-values - Numpy array
-        - y-values - Numpy array
-        Note: same length for both! - If not, the shorter length is taken 
+    Input:
+    - x-values - Numpy array
+    - y-values - Numpy array
+    Note: same length for both! - If not, the shorter length is taken 
 
-        Optional:
-        - Angle of the x-sensor - measured in degrees, clockwise from North (0) 
-        - Angle of the y-sensor - measured in degrees, clockwise from North (0) 
+    Optional:
+    - Angle of the x-sensor - measured in degrees, clockwise from North (0) 
+    - Angle of the y-sensor - measured in degrees, clockwise from North (0) 
 
-        Output:
-        - corrected x-values (North)
-        - corrected y-values (East)
+    Output:
+    - corrected x-values (North)
+    - corrected y-values (East)
     """
 
     x_values = np.array(x_values)
diff --git a/tests/analysis/test_geometry.py b/tests/analysis/test_geometry.py
index 62ae3279b..d8dd69600 100644
--- a/tests/analysis/test_geometry.py
+++ b/tests/analysis/test_geometry.py
@@ -33,24 +33,26 @@ def test_get_dimensionality_from_edi_file(self):
 
     def test_get_eccentricity_from_edi_file(self):
         mt_obj = MT(os.path.normpath(os.path.join(TEST_MTPY_ROOT, "examples/data/edi_files/pb42c.edi")))
-        eccentricity_pb42c = (np.array([0.01675639, 0.01038589, 0.00527011, 0.00638819, 0.01483804,
-                                        0.00385233, 0.00513294, 0.00403781, 0.02862114, 0.02689821,
-                                        0.01425044, 0.05686524, 0.05742524, 0.02696736, 0.0275285,
-                                        0.03647819, 0.04721932, 0.06336521, 0.12789841, 0.16409303,
-                                        0.20630821, 0.34261225, 0.3967886, 0.51629705, 0.56645987,
-                                        0.52558696, 0.46954261, 0.48028767, 0.47490701, 0.44927612,
-                                        0.45185046, 0.44143159, 0.43570377, 0.41537978, 0.40546014,
-                                        0.38785478, 0.37174031, 0.34534557, 0.35510941, 0.32282644,
-                                        0.28501461, 0.22463964, 0.20683855]),
-                              np.array([0.01648335, 0.01850986, 0.02356282, 0.02351531, 0.0245181 ,
-                                        0.0295284 , 0.03050838, 0.03262934, 0.04868243, 0.04905988,
-                                        0.04742146, 0.04365671, 0.06595206, 0.09502736, 0.11485881,
-                                        0.1330897 , 0.16488362, 0.18365985, 0.2233313 , 0.32718984,
-                                        0.35370171, 0.38087385, 0.55555244, 0.59755086, 0.60360704,
-                                        0.82315127, 0.75538697, 0.73747297, 0.46115785, 0.85300128,
-                                        0.95439195, 0.4485694 , 0.42871072, 0.43533434, 0.43440266,
-                                        0.4950619 , 0.56743593, 0.64693597, 0.74768197, 0.90040195,
-                                        1.01527094, 1.15623132, 1.34441602]))
+        eccentricity_pb42c = (np.array(
+            [0.01675639, 0.01038589, 0.00527011, 0.00638819, 0.01483804,
+             0.00385233, 0.00513294, 0.00403781, 0.02862114, 0.02689821,
+             0.01425044, 0.05686524, 0.05742524, 0.02696736, 0.0275285 ,
+             0.03647819, 0.04721932, 0.06336521, 0.12789841, 0.16409303,
+             0.20630821, 0.34261225, 0.3967886 , 0.51629705, 0.56645987,
+             0.52558696, 0.46954261, 0.48028767, 0.47490701, 0.44927612,
+             0.45185046, 0.44143159, 0.43570377, 0.41537978, 0.40546014,
+             0.38785478, 0.37174031, 0.34534557, 0.35510941, 0.32282644,
+             0.28501461, 0.22463964, 0.20683855]),
+            np.array(
+                [1.09157801, 1.88513021, 2.52461839, 2.82831998, 2.48686079,
+                 2.78766939, 2.37250674, 1.16767867, 2.82659526, 1.41909434,
+                 1.57463881, 2.80007766, 2.73439937, 2.71321983, 2.82682667,
+                 2.25051091, 1.435991  , 0.43231073, 0.81874286, 1.62593021,
+                 2.4864899 , 2.91852175, 3.11049976, 3.48724331, 3.6630486 ,
+                 3.54297703, 3.40317309, 3.4175019 , 3.3910377 , 3.28570849,
+                 3.33766706, 3.31489955, 3.32710981, 3.29740796, 3.28179559,
+                 3.24496124, 3.2052965 , 3.13248919, 3.11445279, 2.98672687,
+                 2.85242595, 2.79044042, 2.73406595]))
 
         eccentricity = mtg.eccentricity(z_object=mt_obj.Z)
         for i in range(2):
diff --git a/tests/core/test_z.py b/tests/core/test_z.py
index 1119a9b24..ba640fea8 100644
--- a/tests/core/test_z.py
+++ b/tests/core/test_z.py
@@ -18,14 +18,15 @@ def setUp(self):
     def test_rotate(self):
 
         alpharad = np.deg2rad(self.rotation_angle)
-        rotation_matrix = np.matrix([[np.cos(alpharad),np.sin(alpharad)],
-                                     [-np.sin(alpharad),np.cos(alpharad)]])
+        rotation_matrix = np.matrix([[np.cos(alpharad), -np.sin(alpharad)],
+                                     [np.sin(alpharad), np.cos(alpharad)]])
         
         z_array = self.MT.Z.z.copy()
         
         ztest = z_array[0]
         
-        ztest_rot = np.ma.dot(np.ma.dot(rotation_matrix,ztest),rotation_matrix.T)
+        ztest_rot = np.ma.dot(np.ma.dot(rotation_matrix,ztest), 
+                              rotation_matrix.T)
         
         self.MT.Z.rotate(self.rotation_angle)