% Esercizio 0: Estimate the radiogenic heat production of a granite between
% 0 and 4 Gyr and at specific time of 1.5 Gyr

% Radiogenic heat from each radioactive isotope (from literature) 

A_U238= 96.7;   % Uranium 238 (microWatt/kg)
A_Th232= 26.3;   % Thorium 232 (microWatt/kg)
A_K40= 3.5E-03;   % Potassium 40 (microWatt/kg)

% Half time of each radioactive isotope (from literature) 

tau_U238=4.468E09; % half time  Uranium 238
tau_U235=0.704E09;  % half time  Uranium 235
tau_Th232=14E09;  % half time  Thorium 232
tau_K40=1.25E09;    % half time  Potassium 40

% Write the abundance of each isotope
P_U238=0.993; P_U235=0.0071; P_Th232=1; P_k40=0.00012;

% Write the rock's density (kg/m3)
rho=2700; 

t=0:1E06:4E09; % time

c_U238=4.3*1E-06; % Uranium 238 concentration
c_Th232=12.2*1E-06; % Thorium concentration
c_K40=5.1*1E-02; % Potassium concentration

A_El= [A_U238, A_Th232, A_K40]; % Radiogenic heat

tau=[tau_U238, tau_Th232,tau_K40]; % Half time of the three radioactive isotopes

c(1:3)=[c_U238, c_Th232,c_K40]; % Concentration of radioactive isotopes

% A for two specific times
t1=1.5E09;

% Concentration at a specific time t1
ct1= c.*exp(t1*log(2)./tau); 
ct0=c;

A0=rho*ct0*A_El';
A1=rho*ct1*A_El';

disp(['A al tempo attuale: ',num2str(A0),' microW/m3; A al tempo ',num2str(t1/1E09),' Gyr : ',num2str(A1),' microW/m3'])


% Initialize the matrix
MagA_element=zeros(length(t),4);
Magct=zeros(length(t),3);

figure

% Estimate the variation of the radiogenic heat during geological time

for j=1:3

 % Estimate the variation of concentration for each isotope with time
 
   ct= c(j).*exp(t'.*log(2)./tau(j)); 
   Magct(:,j)=ct;
   
% Estimate the variation of radiogenic heat for each isotope with time

  A_element=rho.*ct.*A_El(j);
  MagA_element(:,j)= A_element;
  
% Plot the results
  
  plot(t./10^9,MagA_element(:,j),'LineWidth',2)
  hold on
  
end

% Estimate the sum of the radiogenic heat production of the three isotopes and plot it

Asum=sum(MagA_element');
MagA_element(:,4)=Asum;
  
hold on
plot(t./10^9,MagA_element(:,4),'LineWidth',2)
  
xlabel('time (Gyr)','FontSize',20,'FontWeight','bold')
ylabel('A (microW/m3)','FontSize',20,'FontWeight','bold')
title('Radiogenic Heat production','FontSize',20,'FontWeight','bold')
legend('A U238', 'A Th232', 'A K40', 'A (summed)','FontSize',20,'FontWeight','bold')