%% Uppgift 8
clear all
clf
clc

global CO2Emissions  CO2ConcRCP45 A tau_init k M_0

% Ladda upp data

utslappRCP45;
koncentrationerRCP45;
beta = 0.35;
omvandlingsfaktor = 0.469;
A = [0.113, 0.213, 0.258, 0.273, 0.1430];
tau_init = [2.0, 12.2, 50.4, 243.3, inf];
M_0 = 278.05158;
k = 2*3.06e-3;

%Förindustriella begynnelsevärden

B_10 = 600;
B_20 = 600;
B_30 = 1500;
B_40 = 38000;
F_120 = 60;
F_210 = 15;
F_310 = 45;
F_230 = 45;
alpha_12 = F_120/B_10;
alpha_13 = 0;
alpha_21 = F_210/B_20;
alpha_23 = F_230/B_20;
alpha_31 = F_310/B_30;
alpha_32 = 0;
NPP0 = F_120;

B1 = zeros(1, length(CO2Emissions));
B2 = zeros(1, length(CO2Emissions));
B3 = zeros(1, length(CO2Emissions));
B4 = zeros(1, length(CO2Emissions));
B1(1)=B_10; % Begynnelsevärden
B2(1)=B_20; % Begynnelsevärden
B3(1)=B_30; % Begynnelsevärden
B4(1)=B_40; % Begynnelsevärden
B1_dot = zeros(1, length(CO2Emissions)+1);
B1_dot_tilde = zeros(1, length(CO2Emissions));
B2_dot = 0;
B3_dot = 0;
inre_summa = zeros(1, 5);
yttre_summa = zeros(1,length(CO2Emissions));
impulssvar = zeros(1,length(CO2Emissions));
dot = 0;

for t = 1:length(CO2Emissions)

    NPP = NPP0*(1 + beta*log(B1(t)/B_10));
    B1_dot(t) = alpha_31*B3(t) + alpha_21*B2(t) - NPP + CO2Emissions(t);
    B2_dot = NPP -(alpha_23+alpha_21)*B2(t);
    B2(t+1) = B2(t) + B2_dot;
    B3_dot = alpha_23*B2(t)-alpha_31*B3(t);
    B3(t+1) = B3(t)+ B3_dot;

    for tilde = 1:t
        for i = 1:5
            inre_summa(i) = sum(A(i)*exp(-(t-tilde)/(tau_init(i)*k*sum(B1_dot(1:t)))));
        end
        impulssvar(tilde) = sum(inre_summa)*B1_dot(tilde);
    end

    if t ~= length(CO2Emissions)
        B1(t+1) = B_10 + sum(impulssvar) + sum(A) * (alpha_31*B3(t+1) + alpha_21*B2(t+1) - NPP + CO2Emissions(t+1));
    else
        B1(t+1) = B_10 + sum(impulssvar);
    end

end


% Ta bort sista elementet för att få längd 736

B1 = B1(1:end-1);
B2 = B2(1:end-1);
B3 = B3(1:end-1);

T = linspace(1765,2500,736);

subplot(1,2,1)
plot(T,B1)
hold on
plot(T,B2)
xlim([min(T)  1765+length(T)])
ylim([600 1600])
title('CO2 och kol i box 1 & 2 sedan 1765.')
xlabel('Tid (År)')
ylabel('Kolstock (GtC)')
legend({'B_1 = Atmosfär','B_2 = Ovan jord'},'Location','best')
set(gca,'FontSize',10,'TickLabelInterpreter','latex')
grid,box on

subplot(1,2,2)
plot(T,B3)
xlim([min(T)  1765+length(T)])
ylim([min(B3) max(B3)])
title('Kol i box 3 sedan 1765.')
xlabel('Tid (År)')
ylabel('Kolstock (GtC)')
legend({'B_3 = Under jord'},'Location','best')
set(gca,'FontSize',10,'TickLabelInterpreter','latex')
grid,box on
hold off

%% Figur för koncentrationen
clf
konc = B1*omvandlingsfaktor;

plot(T,CO2ConcRCP45)
hold on
plot(T,konc)
xlim([min(T)  1765+length(T)])
%ylim([min(konc) max(CO2ConcRCP45)+10])
title('CO2 koncentration i atmosfären sedan 1765.')
xlabel('Tid (År)')
ylabel('CO_2-koncentration(PPM)')
set(gca,'FontSize',10,'TickLabelInterpreter','latex')
legend({'Data','Modell'},'Location','northwest')
grid,box on
hold off

%% havets kolstock vid år 2100
H = 38000 + sum(CO2Emissions(1:335)) - (B1(336)-B1(1)) - (B2(336)-B2(1)) - (B3(336)-B3(1))