diff --git a/examples/ConfigFiles/Phase_Isotropic_FreeSurface.config b/examples/ConfigFiles/Phase_Isotropic_FreeSurface.config
new file mode 100644
index 000000000..80f65ff46
--- /dev/null
+++ b/examples/ConfigFiles/Phase_Isotropic_FreeSurface.config
@@ -0,0 +1,25 @@
+[Air]
+
+## Isotropic Material model to simulate free surfaces ##
+## For more information see paper Maiti+Eisenlohr2018, Scripta Materialia, 
+## "Fourier-based spectral method solution to finite strain crystal plasticity with free surfaces"
+
+elasticity              hooke
+plasticity              isotropic
+
+/dilatation/
+
+(output)                flowstress
+(output)                strainrate
+
+lattice_structure       isotropic
+c11                     110.9e9
+c12                     0.0
+taylorfactor            3
+tau0                    31e6
+gdot0                   0.001
+n                       20
+h0                      75e6
+tausat                  63e6
+w0                      2.25
+atol_resistance         1
\ No newline at end of file
diff --git a/processing/post/addSpectralDecomposition.py b/processing/post/addSpectralDecomposition.py
index 3c44a16bc..f02133f1c 100755
--- a/processing/post/addSpectralDecomposition.py
+++ b/processing/post/addSpectralDecomposition.py
@@ -22,7 +22,13 @@ parser.add_option('-t','--tensor',
                   dest = 'tensor',
                   action = 'extend', metavar = '<string LIST>',
                   help = 'heading of columns containing tensor field values')
+parser.add_option('--no-check',
+                  dest = 'rh',
+                  action = 'store_false',
+                  help = 'skip check for right-handed eigenvector basis')
 
+parser.set_defaults(rh = True,
+                    )
 (options,filenames) = parser.parse_args()
 
 if options.tensor is None:
@@ -57,8 +63,10 @@ for name in filenames:
       if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type))
       else:
         items[type]['column'].append(table.label_index(what))
-        table.labels_append(['{}_eigval({})'.format(i+1,what) for i in range(3)])                   # extend ASCII header with new labels
-        table.labels_append(['{}_eigvec({})'.format(i+1,what) for i in range(9)])                   # extend ASCII header with new labels
+        for order in ['Min','Mid','Max']:
+          table.labels_append(['eigval{}({})'.format(order,what)])                                         # extend ASCII header with new labels
+        for order in ['Min','Mid','Max']:
+          table.labels_append(['{}_eigvec{}({})'.format(i+1,order,what) for i in range(3)])                # extend ASCII header with new labels
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:
@@ -78,10 +86,10 @@ for name in filenames:
     for type, data in items.iteritems():
       for column in data['column']:
         (u,v) = np.linalg.eigh(np.array(map(float,table.data[column:column+data['dim']])).reshape(data['shape']))
-        if np.dot(np.cross(v[:,0], v[:,1]), v[:,2]) < 0.0 : v[:, 2] *= -1.0                         # ensure right-handed coordinate system
-        table.data_append(list(u))
-        table.data_append(list(v.transpose().reshape(data['dim'])))
-    outputAlive = table.data_write()                                                                # output processed line
+        if options.rh and np.dot(np.cross(v[:,0], v[:,1]), v[:,2]) < 0.0 : v[:, 2] *= -1.0          # ensure right-handed eigenvector basis
+        table.data_append(list(u))                                                                  # vector of max,mid,min eigval
+        table.data_append(list(v.transpose().reshape(data['dim'])))                                 # 3x3=9 combo vector of max,mid,min eigvec coordinates
+    outputAlive = table.data_write()                                                                # output processed line in accordance with column labeling
 
 # ------------------------------------------ output finalization -----------------------------------
 
diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90
index 36529630c..0b10f5e3f 100644
--- a/src/plastic_isotropic.f90
+++ b/src/plastic_isotropic.f90
@@ -60,15 +60,18 @@ module plastic_isotropic
      flowstress, &
      accumulatedShear
  end type
+
  type, private :: tIsotropicAbsTol                                                                  !< internal alias for abs tolerance in state
    real(pReal), pointer :: &                                                                        ! scalars
      flowstress, &
      accumulatedShear
  end type
+
  type(tIsotropicState), allocatable, dimension(:), private :: &                                       !< state aliases per instance
    state, &
    state0, &
    dotState
+
  type(tIsotropicAbsTol), allocatable, dimension(:), private :: &                                       !< state aliases per instance
    stateAbsTol
 
@@ -140,7 +143,7 @@ subroutine plastic_isotropic_init(fileUnit)
    extmsg    = ''
  character(len=64) :: &
    outputtag = ''
-  integer(pInt) :: NipcMyPhase
+ integer(pInt) :: NipcMyPhase
 
  write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_ISOTROPIC_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@@ -293,9 +296,9 @@ subroutine plastic_isotropic_init(fileUnit)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-     sizeState    = 2_pInt                                                                           ! flowstress, accumulated_shear
-     sizeDotState = sizeState                                                                        ! both evolve
+     sizeDotState   = 2_pInt                                                                         ! flowstress, accumulated_shear
      sizeDeltaState = 0_pInt                                                                         ! no sudden jumps in state
+     sizeState      = sizeDotState + sizeDeltaState
      plasticState(phase)%sizeState = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
      plasticState(phase)%sizeDeltaState = sizeDeltaState
@@ -418,7 +421,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
              * ( sqrt(1.5_pReal) * norm_Tstar_dev / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
              **param(instance)%n
 
-   Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/param(instance)%fTaylor
+   Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/param(instance)%fTaylor 
 
    if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
        .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &