%% 2a - Tidigare kod som krävs.
clf
clc
clear all
global CO2Emissions  CO2ConcRCP45 A tau_init k M_0 beta

% Ladda upp data

utslappRCP45;
koncentrationerRCP45;
radiativeForcingRCP45;
NASATAnomali;

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 = 3.06e-3;

% lös diffekv för olika scenarion
[B1, ~, ~] = uppg2a(0);
[B_i, ~, ~] = uppg2a(1);
[B_ii, ~, ~] = uppg2a(2);
[B_iii, ~, ~] = uppg2a(3);

% beräkna RF

CO_2Conc = B1*omvandlingsfaktor; % Koncentrationsprofil fr?n modell i labb 1
CO_2Conc_i = B_i*omvandlingsfaktor;
CO_2Conc_ii = B_ii*omvandlingsfaktor;
CO_2Conc_iii = B_iii*omvandlingsfaktor;

RF_CO2 = zeros(1, 2100-1765+1); % W/m^2 = J/(s m^2)
RF_CO2_i = zeros(1, 2100-1765+1);
RF_CO2_ii = zeros(1, 2100-1765+1);
RF_CO2_iii = zeros(1, 2100-1765+1);

init_CO2_conc = M_0;

T = 1765:2100;

for i = 1:length(2100-1765+1)
    RF_CO2(i) = 5.35*log(CO_2Conc(i)/init_CO2_conc);
    RF_CO2_i(i) = 5.35*log(CO_2Conc_i(i)/init_CO2_conc);
    RF_CO2_ii(i) = 5.35*log(CO_2Conc_ii(i)/init_CO2_conc);
    RF_CO2_iii(i) = 5.35*log(CO_2Conc_iii(i)/init_CO2_conc);
end

s = 0.9;
RF_andra_amnen = totRadForcExclCO2AndAerosols(1:(2100-1765+1));

RF_aerosoler = totRadForcAerosols(1:(2100-1765+1));

RF_total = RF_andra_amnen + s*RF_aerosoler + RF_CO2;
RF_total_i = RF_andra_amnen + s*RF_aerosoler + RF_CO2_i;
RF_total_ii = RF_andra_amnen + s*RF_aerosoler + RF_CO2_ii;
RF_total_iii = RF_andra_amnen + s*RF_aerosoler + RF_CO2_iii;

% beräkna temperaturer
lambda = 0.8;
kappa = 9;
c = 4186;
rho = 1020;
h = 50;
d = 200;
global C1 C2
C1 = c *  rho * h / (60 * 60 * 24 * 365);
C2 = c *  rho * d / (60 * 60 * 24 * 365);

% Lös differentialekvation
T = 0:(2100-1765+1);

[Psi1, ~] = solve2(T, RF_total, lambda, kappa);
[Psi_i, ~] = solve2(T, RF_total_i, lambda, kappa);
[Psi_ii, ~] = solve2(T, RF_total_ii, lambda, kappa);
[Psi_iii, ~] = solve2(T, RF_total_iii, lambda, kappa);

meanTemp = 0; % beräknad medeltemperatur mellan 1951-1980 från modellen

for t = (1951-1765+1):(1980-1765+1)
    meanTemp = meanTemp + Psi1(t);
end
meanTemp = meanTemp / (1980 - 1951 + 1);

T = (1765-1765+1):(2100-1765+1);
plot(1765:2100, Psi1(T) - meanTemp)
hold on;
plot(1765:2100, Psi_i(T) - meanTemp)
plot(1765:2100, Psi_ii(T) - meanTemp)
plot(1765:2100, Psi_iii(T) - meanTemp)
title('Temperaturförändringen över tid.')
xlim([1765 2100])
xlabel('Tid (År)')
ylabel('\Delta T')
legend({'Tidigare modell','Modell i)','Modell ii)', 'Modell iii)'},'Location','Best')

function [B1, B2, B3] = uppg2a(scenario)
    global CO2Emissions A tau_init k beta
    %Förindustriella begynnelsevärden
    B_10 = 600;
    B_20 = 600;
    B_30 = 1500;
    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));
    B1(1)=B_10; % Begynnelsev?rden
    B2(1)=B_20; % Begynnelsev?rden
    B3(1)=B_30; % Begynnelsev?rden
    B1_dot = zeros(1, length(CO2Emissions)+1);
    B2_dot = 0;
    B3_dot = 0;
    inre_summa = zeros(1, 5);
    impulssvar = zeros(1,length(CO2Emissions));

    U_2020 = CO2Emissions(2020-1765+1);
    U = zeros(1, length(2100 - 1765 + 1));
    for t = 1:(2020-1765+1)
        U(t) = CO2Emissions(t);
    end

    if scenario == 0 % utsläpp som vanligt
        for t = (2021-1765+1):(2100-1765+1)
            U(t) = CO2Emissions(t);
        end

    elseif scenario == 1 % i)
        decrement = -U_2020 / (2070-2020);
        for t = (2021-1765+1):(2070-1765+1)
            U(t) = U(t-1) - decrement;
        end
        for t = (2070-1765+1):(2100-1765+1)
            U(t) = 0;
        end

    elseif scenario == 2 % ii)
        for t = (2021-1765+1):(2100-1765+1)
            U(t) = U_2020;
        end

    elseif scenario == 3 % iii)
        increment = 2 * U_2020 / (2100-1765);
        for t = (2021-1765+1):(2100-1765+1)
            U(t) = U(t-1) + increment;
        end

    end

    for t = 1:(2100-1765+1)

        NPP = NPP0*(1 + beta*log(B1(t)/B_10));
        B1_dot(t) = alpha_31*B3(t) + alpha_21*B2(t) - NPP + U(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);
end

% solve eq. (2, 3) using the forward Euler method
function [Psi1, Psi2] = solve2(T, RF_total, lambda, kappa)
    global C1 C2
    Psi1 = zeros(1, length(T));
    Psi2 = zeros(1, length(T));
    Psi1(1) = 0;
    Psi2(1) = 0;
    for t = 1:max(T)
        Psi1(t+1) = Psi1(t) + (RF_total(t) - Psi1(t)/lambda - kappa*(Psi1(t) - Psi2(t))) / C1;
        Psi2(t+1) = Psi2(t) + kappa*(Psi1(t) - Psi2(t)) / C2;
    end
end