Radiant cooling approach to calculate the age of the Earth

AgeOfEarth.m

function [age,ri,sio,sii,vi,dr,A,T,checkIndex,egen] = ageOfEarth(n,dt)
clc
format compact
fprintf('Running ageOfEarth with parameters %i cells and a time step of %i\n',n,dt); 


% This function will use an implicit time differencing scheme to solve for
% the age of the earth given some stated initial and boundary conditions

% Passed Variables
% n = number of cells
% dt = time step [x years]

% CONSTANTS (assuming no variation with T, pressure)
r = 6.378*10^6; % Radius of Earth [m]
% a = 1.2*10^(-6); % Alpha, diffusivity of Earth [(m^2)/s] - not useful as seconds
a = 37.8432; % Alpha, diffusivity of Earth [(m^2)/yr] - per year is better
p = 3.0*10^3; % Rho, density of Earth [kg/(m^3)]
c = 1.0*10^3; % Heat capacity of Earth [J/((kg^2)K)]
k = 3.7; % Thermal conductivity of Earth [W/(m*K)] 
nu = 2.2;%/(10^6); % An approximation of the concentration of uranium in Earth
nth = 8;%/(10^6); % An approximation of the concentration of thorium in Earth
radStart = -4540000000; % Approximate the start of radioactive decay as roughly 5 billion years ago [(Billion years)]

% CONDITIONS 
To = 2000; % Initial Temp of Earth is a uniform 2000 °C, in Kelvin
Ts = 0; % Surface Temp of Earth is a constant 0 °C, in Kelvin
% The crust thickness we'll use is 25km with a temp gradient of 0.037 °C/m
% So at ri = r - 25000; T ~= 1°C ~= 274 K (since we set Ts = 0 °C)
% To find i where ri = r - 25000, we set 
checkIndex = 0; % and then save the first index where ri >= r - 25000 in the loop below
depth = 15000; % Approximate epth of the solid crust [meters]


% Getting some useful things
dr = r/(n+0.5); % Get radius step = r(Earth) / number of cells + half cell width
ri = zeros(n,1); % Make empty radius step vector of length n
sio = zeros(n,1); % Make empty cell inner surface area vector of length n
sii = zeros(n,1); % Make empty cell outer surface area vector of length n 
vi = zeros(n,1); % Make empty cell volume vector of length n
fprintf('Populating ri,sio,sii,vi vectors...\n');
for i = 1:n; % Loop through every cell and populate the vectors
    ri(i) = (i-.5)*dr;
    if checkIndex == 0; % check and make sure we haven't saved an index yet
        if ri(i) >= r - depth; % check and see if the radius is what we're interested in
            checkIndex = i; % save the index if the above two conditions are met
        end
    end
    sio(i) = 4*pi*(dr^2)*(i^2);
    sii(i) = 4*pi*(dr^2)*((i-1)^2);
    vi(i) = 4*pi*(dr^2)*(i*(i-1)+(1/3));
end
fprintf('Done.\n');

saveEgen = 0;


% Build matrix [A] by leading, upper, and lower diagonals
A = zeros(n); % Make empty matrix of size n x n
fprintf('Populating matrix [A]\n');
for i = 1:n; % Loop through every cell and populate the leading diagonal
    A(i,i) = 1+((a*dt)/(dr^2))*(((i^2)+(i-1)^2)/(i*(i-1)+(1/3)));
end
for i = 1:n-1; % Loop through every cell to populate the upper and lower diagonals
    A(i,i+1) = -((a*dt)/(dr^2))*((i^2)/(i*(i-1)+(1/3))); % Calc the upper diag
    A(i+1,i) = -((a*dt)/(dr^2))*(((i-1)^2)/(i*(i-1)+(1/3))); % Calc the lower diag
end
A = sparse(A);
fprintf('Done.\n');

% Build initial temperature vector [T^n]
T = zeros(n,1); % Make an empty n x 1 temperature array
T(:,1) = To; % Fill the temperature array with the initial temperature
T(n,1) = Ts; % Surface temperature is maintained at 0 °C.


% Now we do the actual math to determine every [T^(n+1)] vector
% We want to loop this calculation, taking every result for T^(n+1) and
% recalculating by plugging it back in as T^n until we reach the point
% where T(checkIndex) <= 1 °C. This will mean that we've looped through
% enough time steps that the total time elapsed has allowed the crust to
% cool to approximately the same state as it exist today (by our estimation
% above).

% Here, we'll save every result as a new column in the T arry so that the
% progression of temperatures can be seen. This will be memory intensive. 
% Then we can just return the time elapsed as dt*number of colums.

fprintf('Starting the calculation...\n');
year = 1; % t = 1 for the first calculation, this is our loop counter
while T(checkIndex,year) > Ts+(.0037*depth); % Run the calculation loop until the crust is cool enough 
    %egen = nu*(0.01275*exp(-0.9789*(radStart+((year*dt)))/10^9)+0.3066*exp(-0.1551*(radStart+((year*dt)))/10^9))+nth*0.08297*exp(-0.0493*(radStart+((year*dt)))/10^9); % Get the raw heat generation per billion years assuming uniform radioactive unit concentrations
    %saveEgen = saveEgen + (dt)*(egen/(p*c)); % Track the total heat generated
    %T(:,year) = T(:,year)+(dt)*(egen/(p*c)); % Add the radioactive generated heat to the Temperature vector
    T(:,year+1) = zeros(n,1); % Append a zero column to the matrix to hold the T^(n+1) solution
    T(:,year+1) = A\T(:,year); % Solve the matrix system, store the result as a new column in the temperature matrix
    T(n,year+1) = Ts; % Keep the surface temperature at 0 °C
    T(1,year+1) = To; % Keep core temperature at 2000 °C
    fprintf('Step %i | Crust depth temperature is %.1f °C | Current age: %i \n',year,T(checkIndex,year),dt*year);
    year = year + 1; % Increment the count
    for i = 1:n-1
        if T(i,year) < T(i+1,year);
            T(i,year) = To;
        end
    end
end
age = dt*year;
fprintf('Finished.\nThe age of the Earth is %i years\n',age);
% surf(T);colormap('hot');colorbar

end % Close the function