% Neil Edelman (110121860)
function b = sbolt(max, steps);
% computes the Stephen-Boltzmann constant numerically directly
% from the formula using approximation over a range up to max
% with steps steps e.g.
% >> sbolt(10000,100000)
% Estimated a remainder of:     0
% ans = 5.6709e-08

dn = max / steps;
int = 0;
n = 0;
while n < max
    % move n to the middle of the step
    n = n + dn / 2;
    % calculate the value of the integrand at n
    val = n.^3 / (exp(n) - 1);
    % add the value to the integral
    int = int + val * dn;
    % calculate the derivative there
    % diff = -n.^2*((n-3)*exp(n)+3)/(exp(n)-1).^2;
    % move n to the end of the step
    n = n + dn / 2;
end
% assume that the function is bounded below f(x) <= exp(-x)
% then we have int(c*exp(-x), x=n..infinity) = c*exp(-n);
int = int + val * exp(-n) / 2;

disp('Estimated a remainder of:');
disp(val * exp(-n) / 2);

k = 1.3807 * 10.^(-23);
c = 2.9979 * 10.^8;
hbar = 1.0546 * 10.^(-34);

b = k.^4 / (4 * pi.^2 * c.^2 * hbar.^3) * int;
