diff --git a/camb/reionization.py b/camb/reionization.py
index b1549c5d..24f934b6 100644
--- a/camb/reionization.py
+++ b/camb/reionization.py
@@ -10,6 +10,12 @@ class ReionizationModel(F2003Class):
 
     _fields_ = (
         ("Reionization", c_bool, "Is there reionization? (can be off for matter power, which is independent of it)"),
+        (
+            "include_heating",
+            c_bool,
+            "Whether to include smooth heating up to ~1e4K following x_e shape (default False)",
+        ),
+        ("heating_temperature", c_double, "Heating temperature in K applied during reionization (default 1e4 K)"),
     )
 
 
diff --git a/fortran/classes.f90 b/fortran/classes.f90
index 5ce85872..a33b6822 100644
--- a/fortran/classes.f90
+++ b/fortran/classes.f90
@@ -109,6 +109,8 @@ module classes
 
     Type, extends(TCambComponent) :: TReionizationModel
         logical  :: Reionization = .true.
+        logical  :: include_heating = .false.
+        real(dl) :: heating_temperature = 1.0d4
     contains
     procedure :: Init => TReionizationModel_Init
     procedure :: x_e => TReionizationModel_xe
diff --git a/fortran/results.f90 b/fortran/results.f90
index cbbf6382..c5792928 100644
--- a/fortran/results.f90
+++ b/fortran/results.f90
@@ -1709,6 +1709,8 @@ subroutine Thermo_Init(this, State,taumin)
     real(dl), allocatable ::  taus(:)
     real(dl), allocatable :: xe_a(:), sdotmu(:), opts(:)
     real(dl), allocatable :: scale_factors(:), times(:), dt(:)
+    ! Variables for simple reionization heating model and cs^2
+    real(dl) :: T_reion, denom, yheat, Tg, muinv, cs2_orig, cs2_heat, Tb_orig
     Type(TCubicSpline) :: dotmuSp
     integer ninverse, nlin
     real(dl) dlna, zstar_min, zstar_max
@@ -1985,17 +1987,41 @@ subroutine Thermo_Init(this, State,taumin)
             this%xe(i)=xe_a(i)
         end if
 
-        !  approximate Baryon sound speed squared (over c**2).
-        !  Use pre-reionization ionization fraction for cs2-related terms for consistency
-        !  (not correct, but avoids odd behaviour at very high k)
-        !  https://github.com/cmbant/CAMB/issues/171
+        ! Reionization heating model: smoothly raise Tb to T_reion following x_e shape
+        ! and map cs^2 smoothly from the original formula to ideal-gas form.
+        ! Use mu^{-1} = (1 - 0.75 Y_He) + (1 - Y_He) x_e.
+        T_reion = CP%Reion%heating_temperature
+
+        ! Original cs^2 (pre-reionization / default)
         fe=(1._dl-CP%yhe)*xe_a(i)/(1._dl-0.75d0*CP%yhe+(1._dl-CP%yhe)*xe_a(i))
         dtbdla=-2._dl*this%tb(i)
         if (a*this%tb(i)-CP%tcmb < -1e-8) then
             dtbdla= dtbdla -thomc0*fe/adot*(a*this%tb(i)-CP%tcmb)/a**3
         end if
         barssc=barssc0*(1._dl-0.75d0*CP%yhe+(1._dl-CP%yhe)*xe_a(i))
-        this%cs2(i)=barssc*this%tb(i)*(1-dtbdla/this%tb(i)/3._dl)
+        cs2_orig = barssc*this%tb(i)*(1-dtbdla/this%tb(i)/3._dl)
+
+        Tb_orig = this%tb(i)
+        if (CP%Reion%Reionization .and. CP%Reion%include_heating .and. tau > State%reion_tau_start) then
+            ! yheat follows the reionization x_e shape, 0 before reionization, ->1 when x_e~1
+            denom = max(1.0d-12, 1._dl - xe_a(i))
+            yheat = (this%xe(i) - xe_a(i))/denom
+            if (yheat < 0) yheat = 0
+            if (yheat > 1) yheat = 1
+            if (yheat > 0) then
+                Tg = Tb_orig + yheat*(T_reion - Tb_orig)
+                muinv = (1._dl-0.75d0*CP%yhe+(1._dl-CP%yhe)*this%xe(i))
+                cs2_heat = (5._dl/3._dl) * barssc0 * muinv * Tg
+                ! Smooth mapping
+                this%cs2(i) = (1._dl - yheat)*cs2_orig + yheat*cs2_heat
+                this%tb(i) = Tg
+            else
+                this%cs2(i) = cs2_orig
+            end if
+        else
+            ! Before reionization (or if heating disabled), keep original values
+            this%cs2(i) = cs2_orig
+        end if
 
         ! Calculation of the visibility function
         this%dotmu(i)=this%xe(i)*State%akthom/a2