% solve  Phi0 = F(Kt)  to find the age of the earth, where
%     F(z) = (2 theta0/L) sum_{m=0}^infty e^{-(2m+1)^2 pi^2 z/(2L)^2}

MM = 501;  % highest integer in sum; can show a few hundred terms needed
           % to get reasonably close to infinite sum

%%% reasonable constants
L = 6357.0e3;  % m;      source: "earth radius" wikipedia page
K = 1.0e-6;    % m^2/s; "typical value for rock"
               %         source: "thermal diffusivity" wikipedia page
theta0 = 1811; % K; melting point of iron 
               %         source: "melting points of the elements" wikipedia page
Phi0 = 0.025;  % K/m; "geothermal gradient in crust"    source:
               %    http://gpc.edu/~pgore/Earth&Space/GPS/earthinterior.html
Kt = 10.^(7:0.1:12);

% %%% fake constants
% L = 1;  theta0 = 1;  K = 1;  Phi0 = 1;
% Kt = 0:0.01:1;
% %%% solution for t about 0.282 when MM = 5 and 0.281 when MM = 1

ss = zeros(size(Kt));
for m=1:2:MM
   ss = ss + exp(- m * m * pi * pi * Kt / (4 * L * L));
end
FF = (2 * theta0 / L) * ss;

loglog(Kt, Phi0*ones(size(Kt)), '--', Kt, FF)
axis([min(Kt) max(Kt) 1e-3 1])
xlabel('Kt'), ylabel('temperature gradient')
legend('\Phi_0', 'F(Kt)')

Kt_cross = 10^9.22;  % result from blowing up the log-log graph
t_cross = Kt_cross / K;
t_cross_year = t_cross / 31556926  % divide by number of seconds in a year