The scalar radiative transfer equation in the presence of thermal radiation source is solved in detail, using the adding-doubling method; Planck functions within any given layer are assumed to possess constant, linear, or exponential parameterizations with optical thickness. The radiance profile in any zenith direction is calculated directly in terms of matrix inversions. The inputs to the model are the inherent optical properties (layer total single-scattering albedos, scattering phase functions, and optical thickness) along with temperature and altitude profiles, and the top of the atmosphere and ground surface boundary conditions. The algorithm is implemented in a state-of-the-art MATLAB program, with the cosmic microwave background as the source at the upper boundary and a Lambertian surface reflection at the lower boundary. The simulations are validated against the VLIDORT discrete ordinate radiative transfer model. Results are compared in detail for cases with linear and exponential Planck function parameterizations.