diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index dd431ad9b..14e2affb9 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -282,9 +282,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
     ce = ce + 1
     call thermal_conduction_getSource(Tdot,1,ce)
     scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
-                            + thermal_conduction_getMassDensity (ce)* &
-                              thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k)  - &
-                                                                          T_current(i,j,k))&
+                            + homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
                             + mu_ref*T_current(i,j,k)
   enddo; enddo; enddo
 
@@ -314,7 +312,7 @@ subroutine updateReference
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
     K_ref  = K_ref  + thermal_conduction_getConductivity(ce)
-    mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce)
+    mu_ref = mu_ref + homogenization_thermal_mu_T(ce)
   enddo; enddo; enddo
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 7a931c793..591bfc04e 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -143,15 +143,10 @@ module homogenization
       real(pReal), dimension(3,3) :: K
     end function thermal_conduction_getConductivity
 
-    module function thermal_conduction_getSpecificHeat(ce) result(c_P)
+    module function homogenization_thermal_mu_T(ce) result(mu_T)
       integer, intent(in) :: ce
-      real(pReal) :: c_P
-    end function thermal_conduction_getSpecificHeat
-
-    module function thermal_conduction_getMassDensity(ce) result(rho)
-      integer, intent(in) :: ce
-      real(pReal) :: rho
-    end function thermal_conduction_getMassDensity
+      real(pReal) :: mu_T
+    end function homogenization_thermal_mu_T
 
     module subroutine homogenization_thermal_setField(T,dot_T, ce)
       integer, intent(in) :: ce
@@ -194,9 +189,8 @@ module homogenization
   public ::  &
     homogenization_init, &
     materialpoint_stressAndItsTangent, &
-    thermal_conduction_getSpecificHeat, &
+    homogenization_thermal_mu_T, &
     thermal_conduction_getConductivity, &
-    thermal_conduction_getMassDensity, &
     thermal_conduction_getSource, &
     damage_nonlocal_getMobility, &
     damage_nonlocal_getSourceAndItsTangent, &
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 66acf40c2..95b58ec52 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -45,7 +45,7 @@ module subroutine thermal_init()
 
 
   print'(/,a)', ' <<<+-  homogenization:thermal init  -+>>>'
-  print'(/,a)', ' <<<+-  homogenization:thermal:isotemperature   init  -+>>>'
+  print'(/,a)', ' <<<+-  homogenization:thermal:isotemperature init  -+>>>'
 
 
 
@@ -128,10 +128,20 @@ module function thermal_conduction_getConductivity(ce) result(K)
 end function thermal_conduction_getConductivity
 
 
+module function homogenization_thermal_mu_T(ce) result(mu_T)
+  
+  integer, intent(in) :: ce
+  real(pReal) :: mu_T
+
+  mu_T = c_P(ce) * rho(ce)
+
+end function homogenization_thermal_mu_T
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns homogenized specific heat capacity
 !--------------------------------------------------------------------------------------------------
-module function thermal_conduction_getSpecificHeat(ce) result(c_P)
+function c_P(ce)
 
   integer, intent(in) :: ce
   real(pReal) :: c_P
@@ -139,21 +149,20 @@ module function thermal_conduction_getSpecificHeat(ce) result(c_P)
   integer :: co
 
 
-  c_P = 0.0_pReal
-
-  do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
+  c_P = lattice_c_p(material_phaseID(1,ce))
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     c_P = c_P + lattice_c_p(material_phaseID(co,ce))
   enddo
 
   c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end function thermal_conduction_getSpecificHeat
+end function c_P
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief returns homogenized mass density
 !--------------------------------------------------------------------------------------------------
-module function thermal_conduction_getMassDensity(ce) result(rho)
+function rho(ce)
 
   integer, intent(in) :: ce
   real(pReal) :: rho
@@ -161,15 +170,14 @@ module function thermal_conduction_getMassDensity(ce) result(rho)
   integer :: co
 
 
-  rho = 0.0_pReal
-
-  do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
+  rho = lattice_rho(material_phaseID(1,ce))
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     rho = rho + lattice_rho(material_phaseID(co,ce))
   enddo
 
   rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end function thermal_conduction_getMassDensity
+end function rho