diff --git a/src/constants.f90 b/src/constants.f90
index 69bfefdb04569334a61b26950ef9d6ef7fe24e12..52a98cfbbd549e0af1653294892e3b8631f9180d 100644
--- a/src/constants.f90
+++ b/src/constants.f90
@@ -13,7 +13,7 @@ module constants
   real(real12), parameter :: hbar = 1.05457148e-34_real12 !Reduced Planck constant
   real(real12), parameter :: atomic_mass=1.67262158e-27_real12 !Atomic mass of hydrogen
   real(real12), parameter :: avogadros=6.022e23_real12 !Avogadro's number
-  real(real12), parameter :: StefBoltz=5.670373e-8_real12 !Stefan-Boltzmann constant
+  
   real(real12), parameter :: TINY=1.0e-12_real12 !Tiny number to avoid division by zero
   integer(int12) , parameter :: fields=6_int12 !Number of fields in the system
 end module constants
diff --git a/src/constructions.f90 b/src/constructions.f90
index ae076aee95b5181952b4308979e728affe8ad291..5550305b96f230bce8f44354664a3ad9bd2a4de8 100644
--- a/src/constructions.f90
+++ b/src/constructions.f90
@@ -27,7 +27,7 @@ module constructions
      real(real12) :: volume !volume of the block
      real(real12), dimension(3) :: Length !length of the block in 3 dimensions
      integer(int12) :: iheater !whether the block is a heater
-     real(real12) :: kappa, rho, heat_capacity, tau, em !physical properties of the material
+     real(real12) :: kappa, rho, heat_capacity, tau !physical properties of the material
   end type heatblock
 !!Defines the material to be used in the simulation
   type material
@@ -39,7 +39,6 @@ module constructions
      real(real12) :: rho !density of the material
      real(real12) :: sound_speed   !speed of sound in the material
      real(real12) :: tau !relaxation time
-     real(real12) :: em !emissivity of the material
      logical :: source !whether the material is a source of heat ??
   end type material
 end module constructions
diff --git a/src/heating.f90 b/src/heating.f90
index 79ef2c0e9d49c02f27cae39853205ba5a1d2872d..b40d9e2e39c7af6844e04334ebba39bd1e094f8f 100644
--- a/src/heating.f90
+++ b/src/heating.f90
@@ -11,9 +11,9 @@
 !!! Author: Harry Mclean, Frank Davies, Steven Hepplestone
 !!!#################################################################################################
 module Heating
-  use constants, only: real12, int12, pi, StefBoltz
+  use constants, only: real12, int12, pi
   use globe_data, only: Temp_p, Temp_pp, Heat, heated_volume
-  use inputs, only: nx,ny,nz, grid, NA, power_in, time_step, T_System, freq, ntime, T_Bath
+  use inputs, only: nx,ny,nz, grid, NA, power_in, time_step, T_System, freq, ntime
   use materials, only: material
   implicit none
 contains
@@ -35,7 +35,7 @@ contains
     Q = 0._real12
     POWER = power_in
     time = time_step * real(itime,real12)
-    time_pulse = 1E-12_real12
+    time_pulse = 0.5_real12
     heated_volume=0.0
     heated_num=0
 
@@ -89,20 +89,20 @@ contains
                 !------------------------------
                 ! AC oscillatory heating raw, with power correction
                 !------------------------------
-                Q(IA) = POWER * (sin(time * 2.0_real12 * PI * freq)**2.0_real12)
+                Q(IA) = POWER * (sin(time * 2.0_real12 * PI * freq)**2.0_real12)&
+                  +POWER*2.0_real12*PI*freq*tau*sin(2.0_real12*time*2.0_real12*PI*freq)
                 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
              case(5)
                 !------------------------------
-                !  AC oscillatory heating raw, with power correction
+                ! Radiative heating 
                 !------------------------------
-                Q(IA) = POWER * (sin(time * 2.0_real12 * PI * freq)**2.0_real12)&
-                  +POWER*2.0_real12*PI*freq*tau*sin(2.0_real12*time*2.0_real12*PI*freq)
+                ! Q(IA)= e*eps* T**4
                 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
              case(6)
                 !------------------------------
-                ! Pulsed heating
+                ! Step one heating
                 !------------------------------
-               if ((itime .le. time_pulse)) then
+                if (itime == 1) then
                    Q(IA) = POWER
                 else
                    Q(IA) = 0.0
@@ -111,15 +111,12 @@ contains
              end select
 
              !------------------------------
-             ! If emissitivity is not zero, then calculate the radiative heating
+             ! convert from Q~density to Q
              !------------------------------
-              
-               Q(IA) = Q(IA) - grid(ix,iy,iz)%em * grid(ix,iy,iz)%length(1)*&
-                       grid(ix,iy,iz)%length(2)*StefBoltz &
-                       * ((Temp_p(IA)**4.0_real12) - (T_Bath**4.0_real12)) 
+             !Qdens(IA)=Q(IA)
+             Q(IA)=Q(IA)
              !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-
              
              
              !------------------------------
diff --git a/src/inputs.f90 b/src/inputs.f90
index fb40836f8e560973258c1894c4eb1930582f3e1a..160c9968254897f9baf3354a110a8d3dcc361c3f 100644
--- a/src/inputs.f90
+++ b/src/inputs.f90
@@ -525,7 +525,7 @@ subroutine read_mat(unit)
     type(material), dimension(100) :: dum_mat
     character(1024) :: buffer
     integer :: reason, j
-    integer, dimension(8) :: readvarmat
+    integer, dimension(7) :: readvarmat
     integer :: i, index
 
     i=0
@@ -571,7 +571,6 @@ subroutine read_mat(unit)
        CALL assignD(buffer,"rho"          ,dum_mat(i)%rho          ,readvarmat(5))! assign rho
        CALL assignD(buffer,"sound_speed"  ,dum_mat(i)%sound_speed  ,readvarmat(6))! assign sound_speed
        CALL assignD(buffer,"tau"          ,dum_mat(i)%tau          ,readvarmat(7))! assign tau
-       CALL assignD(buffer,"em"           ,dum_mat(i)%em           ,readvarmat(8))! assign emissivity
     end do read
     
     ! Check for duplicate indices
diff --git a/src/material.f90 b/src/material.f90
index 5af32ba2ec974ce0c4f42114480b97afc075e7a6..86e1077c7706aa7c5bc6cf7c0ac06ae23025e186 100644
--- a/src/material.f90
+++ b/src/material.f90
@@ -22,11 +22,11 @@ module materials
 
 contains
 
-subroutine material(imaterial_type,kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em)
+subroutine material(imaterial_type,kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau)
   
    integer(int12), intent(in) ::imaterial_type
    integer(int12) :: i, tmp
-   real(real12), intent(inout) :: kappa3D, kappa, h_conv, heat_capacity, sound_speed, rho, tau, em
+   real(real12), intent(inout) :: kappa3D, kappa, h_conv, heat_capacity, sound_speed, rho, tau
    logical :: found
 
 !!!=============================================
@@ -66,7 +66,6 @@ subroutine material(imaterial_type,kappa,kappa3D,h_conv,heat_capacity,rho,sound_
             rho           = input_materials(i)%rho
             sound_speed   = input_materials(i)%sound_speed
             tau           = input_materials(i)%tau
-            em            = input_materials(i)%em
 
             exit mat_loop
          end if
diff --git a/src/setup.f90 b/src/setup.f90
index a1591b79887db3c6b327df2e1f0b1e9474b58e22..59d2cd91c26546338b264b238e2f445c792f3bf2 100644
--- a/src/setup.f90
+++ b/src/setup.f90
@@ -32,7 +32,7 @@ module setup
 !!!#################################################################################################
    subroutine set_global_variables()
       integer(int12) :: ix,iy,iz,index
-      real(real12) :: kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em
+      real(real12) :: kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau
       allocate(Temp_cur(nx, ny, nz))
       allocate(Temp_p(NA))
       allocate(Temp_pp(NA))
@@ -51,12 +51,11 @@ module setup
             do ix = 1, nx
             index = index + 1
             CALL material(grid(ix,iy,iz)%imaterial_type,&
-                 kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em)
+                 kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau)
             grid(ix,iy,iz)%kappa = kappa
             grid(ix,iy,iz)%rho = rho
             grid(ix,iy,iz)%heat_capacity = heat_capacity
             grid(ix,iy,iz)%tau = tau*inverse_time*inverse_time
-            grid(ix,iy,iz)%em = em
             lin_rhoc(index) = rho*heat_capacity
             if (Check_Stability) CALL stability(kappa, rho, heat_capacity, ix, iy, iz)
             end do