diff --git a/code/IO.f90 b/code/IO.f90
index 59632cd0a..86961a785 100644
--- a/code/IO.f90
+++ b/code/IO.f90
@@ -262,6 +262,7 @@ endsubroutine
    IO_hybridIA(1,i) = deltas(1)*(mod(bin/(steps(3)*steps(2)),steps(1))+center)  ! phi1
    IO_hybridIA(2,i) = deltas(2)*(mod(bin/ steps(3)          ,steps(2))+center)  ! Phi
    IO_hybridIA(3,i) = deltas(3)*(mod(bin                    ,steps(3))+center)  ! phi2
+!    write(6,'(i6,x,a3,i5,x,a5,x,3(f8.5,x))') i,'bin',bin,'euler',IO_hybridIA(:,i)
    binSet(j) = binSet(i)
  enddo
  close(999)
diff --git a/code/crystallite.f90 b/code/crystallite.f90
index 3973c3b19..560515477 100644
--- a/code/crystallite.f90
+++ b/code/crystallite.f90
@@ -86,7 +86,8 @@ subroutine crystallite_init(Temperature)
  use debug, only:            debug_info, &
                              debug_reset
  use math, only:             math_I3, &
-                             math_EulerToR
+                             math_EulerToR, &
+                             math_RToEuler
  use FEsolving, only:        FEsolving_execElem, &
                              FEsolving_execIP
  use mesh, only:             mesh_element, &
@@ -1043,7 +1044,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
               k = kl(1,mod(cycleCounter/iJacoStiffness,9)+1,g,i,e)
               l = kl(2,mod(cycleCounter/iJacoStiffness,9)+1,g,i,e)      
               crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg                             ! ... use tangent dP_ij/dFg_kl
-            elseif (.not. storedConvergenceFlag(g,i,e)) then                                                                        ! if crystallite didnÕt converge before...
+            elseif (.not. storedConvergenceFlag(g,i,e)) then                                                                        ! if crystallite didn�t converge before...
               crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e)                                             ! ... use (elastic) fallback
             endif
       enddo; enddo; enddo
@@ -1549,6 +1550,8 @@ use prec, only:                       pInt, &
                                       pReal
 use math, only:                       math_pDecomposition, &
                                       math_RtoQuaternion, &
+                                      math_RtoEuler, &
+                                      math_QuaternionToEuler, &
                                       math_QuaternionDisorientation, &
                                       inDeg, &
                                       math_qConj
@@ -1557,6 +1560,7 @@ use FEsolving, only:                  FEsolving_execElem, &
 use IO, only:                         IO_warning
 use material, only:                   material_phase, &
                                       homogenization_Ngrains, &
+                                      material_EulerAngles, &
                                       phase_constitution
 use mesh, only:                       mesh_element, &
                                       mesh_ipNeighborhood, &
@@ -1588,6 +1592,7 @@ logical error
 
 
 !$OMP PARALLEL DO
+  n = 0_pInt
   do e = FEsolving_execElem(1),FEsolving_execElem(2)
     do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
       do g = 1,homogenization_Ngrains(mesh_element(3,e))
diff --git a/code/material.f90 b/code/material.f90
index d841001f0..a42207c3c 100644
--- a/code/material.f90
+++ b/code/material.f90
@@ -667,6 +667,9 @@ subroutine material_populateGrains()
            constituentGrain = constituentGrain + myNorientations
 
          endif
+!          do j = 1,myNorientations
+!            write(6,'(i6,x,a5,x,3(f8.5,x))') j,'euler',orientationOfGrain(:,grain+j)
+!          enddo
 ! ----------------------------------------------------------------------------
          symExtension = texture_symmetry(textureID) - 1_pInt
          if (symExtension > 0_pInt) then                                      ! sample symmetry
@@ -687,6 +690,7 @@ subroutine material_populateGrains()
        enddo  ! constituent
 
 ! ----------------------------------------------------------------------------
+
        if (.not. microstructure_elemhomo(micro)) then                           ! unless element homogeneous, reshuffle grains
          do i=1,myNgrains-1                                                     ! walk thru grains
            call random_number(rnd)
@@ -720,6 +724,7 @@ subroutine material_populateGrains()
                material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g)
                material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g)
              end forall
+
              grain = grain + FE_Nips(mesh_element(2,e)) * dGrains               ! wind forward by Nips*NgrainsPerIP
            endif
            ! write (6,*) e
diff --git a/code/math.f90 b/code/math.f90
index 3e2ac16f2..49a3eaeb2 100644
--- a/code/math.f90
+++ b/code/math.f90
@@ -431,7 +431,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
  
 !**************************************************************************
 ! matrix multiplication 33x33 = 3x3
-!**************************************************************************´
+!**************************************************************************�
  pure function math_mul33x33(A,B)  
 
  use prec, only: pReal, pInt
@@ -1256,9 +1256,9 @@ pure function math_transpose3x3(A)
  if(val >  1.0_pReal) val =  1.0_pReal
  if(val < -1.0_pReal) val = -1.0_pReal
      
- math_RtoEuler(2) = acos(val)
+ math_RtoEuler(2) = dacos(val)
 
- if(math_RtoEuler(2) < 1.0e-30_pReal) then
+ if(math_RtoEuler(2) < 1.0e-8_pReal) then
 ! calculate phi2
      math_RtoEuler(3) = 0.0_pReal
 ! calculate phi1
@@ -1266,7 +1266,7 @@ pure function math_transpose3x3(A)
      if(val >  1.0_pReal) val =  1.0_pReal
      if(val < -1.0_pReal) val = -1.0_pReal
      
-     math_RtoEuler(1) = acos(val)
+     math_RtoEuler(1) = dacos(val)
      if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1)
  else
 ! calculate phi2
@@ -1274,14 +1274,14 @@ pure function math_transpose3x3(A)
      if(val >  1.0_pReal) val =  1.0_pReal
      if(val < -1.0_pReal) val = -1.0_pReal
      
-     math_RtoEuler(3) = acos(val)
+     math_RtoEuler(3) = dacos(val)
      if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3)
 ! calculate phi1
      val=-R(3,2)/sin(math_RtoEuler(2))
      if(val >  1.0_pReal) val =  1.0_pReal
      if(val < -1.0_pReal) val = -1.0_pReal
      
-     math_RtoEuler(1) = acos(val)
+     math_RtoEuler(1) = dacos(val)
      if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1)
  end if
  return
@@ -1297,43 +1297,48 @@ pure function math_transpose3x3(A)
  use prec, only: pReal, pInt
  implicit none
 
- real(pReal), dimension (3,3), parameter :: evenPerm = reshape((/1,2,3,2,3,1,3,1,2/),(/3,3/))
-
  real(pReal), dimension (3,3), intent(in) :: R
- real(pReal), dimension(4) :: math_RtoQuaternion
- real(pReal) largestDiagElem,root
- integer(pInt) i,largest
+ real(pReal), dimension(4)   :: absQ,math_RtoQuaternion
+ integer(pInt), dimension(1) :: largest
  
- ! math adopted from http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternions_and_other_representations_of_rotations
- ! assuming that above corresponds to active rotation, we negate rotation angle
+ ! math adopted from http://code.google.com/p/mtex/source/browse/trunk/geometry/geometry_tools/mat2quat.m
+
+ math_RtoQuaternion = 0.0_pReal
+
+ absQ(1) = 0.5_pReal * dsqrt(1_pReal+R(1,1)+R(2,2)+R(3,3))
+ absQ(2) = 0.5_pReal * dsqrt(1_pReal+R(1,1)-R(2,2)-R(3,3))
+ absQ(3) = 0.5_pReal * dsqrt(1_pReal-R(1,1)+R(2,2)-R(3,3))
+ absQ(4) = 0.5_pReal * dsqrt(1_pReal-R(1,1)-R(2,2)+R(3,3))
  
- largestDiagElem = 0.0_pReal
- largest = 0_pInt
+ largest = maxloc(absQ)
+ select case(largest(1))
+   case (1)
+
+      math_RtoQuaternion(2) = R(2,3)-R(3,2)
+      math_RtoQuaternion(3) = R(3,1)-R(1,3)
+      math_RtoQuaternion(4) = R(1,2)-R(2,1)
+   
+   case (2)
+      math_RtoQuaternion(1) = R(2,3)-R(3,2)
+
+      math_RtoQuaternion(3) = R(1,2)+R(2,1)
+      math_RtoQuaternion(4) = R(3,1)+R(1,3)
+   
+   case (3)
+      math_RtoQuaternion(1) = R(3,1)-R(1,3)
+      math_RtoQuaternion(2) = R(1,2)+R(2,1)
+
+      math_RtoQuaternion(4) = R(2,3)+R(3,2)
+   
+   case (4)
+      math_RtoQuaternion (1) = R(1,2)-R(2,1)
+      math_RtoQuaternion (2) = R(3,1)+R(1,3)
+      math_RtoQuaternion (3) = R(3,2)+R(2,3)
  
- do i = 1,3
-   if (dabs(R(i,i)) >= largestDiagElem) then
-     largestDiagElem = dabs(R(i,i))
-     largest = i
-   endif
- enddo
- 
- root = dsqrt( max(0.0_pReal, &
-                   1.0_pReal - (R(1,1)+R(2,2)+R(3,3)) + 2.0_pReal*R(largest,largest)) )
- if (root < 1.0e-8_pReal) then
-   math_RtoQuaternion = 0.0_pReal
-   math_RtoQuaternion(1) = 1.0_pReal
- else
-   math_RtoQuaternion(1) = 0.5_pReal / root * &
-                           ( R(evenPerm(2,largest),evenPerm(3,largest)) - &
-                             R(evenPerm(3,largest),evenPerm(2,largest)) )
-   math_RtoQuaternion(1+evenPerm(1,largest)) = 0.5_pReal * root
-   math_RtoQuaternion(1+evenPerm(2,largest)) = 0.5_pReal / root * &
-                                               ( R(evenPerm(1,largest),evenPerm(2,largest)) + &
-                                                 R(evenPerm(2,largest),evenPerm(1,largest)) )
-   math_RtoQuaternion(1+evenPerm(3,largest)) = 0.5_pReal / root * &
-                                               ( R(evenPerm(3,largest),evenPerm(1,largest)) + &
-                                                 R(evenPerm(1,largest),evenPerm(3,largest)) )
- endif
+ end select
+
+ math_RtoQuaternion = math_RtoQuaternion*0.25_pReal/maxval(absQ)
+ math_RtoQuaternion(largest(1)) = maxval(absQ)
 
  return
  
@@ -1515,19 +1520,24 @@ pure function math_transpose3x3(A)
  real(pReal), dimension(4), intent(in) :: Q
  real(pReal), dimension(3) :: math_QuaternionToEuler
 
+ math_QuaternionToEuler(2) = dacos(1.0_pReal-2.0_pReal*(Q(2)*Q(2)+Q(3)*Q(3)))
 
- math_QuaternionToEuler(1) = atan2(Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)-Q(3)*Q(4))
- if (math_QuaternionToEuler(1) < 0.0_pReal) &
-   math_QuaternionToEuler(1) = math_QuaternionToEuler(1) + 2.0_pReal * pi
+ if (dabs(math_QuaternionToEuler(2)) < 1.0e-3_pReal) then
+   math_QuaternionToEuler(1) = 2.0_pReal*dacos(Q(1))
+   math_QuaternionToEuler(3) = 0.0_pReal
+ else
+   math_QuaternionToEuler(1) = datan2(Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)-Q(3)*Q(4))
+   if (math_QuaternionToEuler(1) < 0.0_pReal) &
+     math_QuaternionToEuler(1) = math_QuaternionToEuler(1) + 2.0_pReal * pi
+
+   math_QuaternionToEuler(3) = datan2(-Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)+Q(3)*Q(4))
+   if (math_QuaternionToEuler(3) < 0.0_pReal) &
+     math_QuaternionToEuler(3) = math_QuaternionToEuler(3) + 2.0_pReal * pi
+ endif
 
- math_QuaternionToEuler(2) = acos(1.0_pReal-2.0_pReal*(Q(2)*Q(2)+Q(3)*Q(3)))
  if (math_QuaternionToEuler(2) < 0.0_pReal) &
    math_QuaternionToEuler(2) = math_QuaternionToEuler(2) + pi
 
- math_QuaternionToEuler(3) = atan2(-Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)+Q(3)*Q(4))
- if (math_QuaternionToEuler(3) < 0.0_pReal) &
-   math_QuaternionToEuler(3) = math_QuaternionToEuler(3) + 2.0_pReal * pi
- 
  return
  
  endfunction