clc
clear
xdel(winsid())

function thetaIMenosMeio = calculaThetaRecuado(Q, pos)
    if((pos == 2) | ((Q(pos)-Q(pos-1)) == 0)) then
        thetaIMenosMeio = return 0;
    end
    thetaIMenosMeio = return (Q(pos-1)-Q(pos-2))/(Q(pos)-Q(pos-1));
end

function thetaIMaisMeio = calculaThetaAvancado(Q, pos)
    if((Q(pos+1)-Q(pos)) == 0) then
        thetaIMaisMeio = return 0;
    end
    thetaIMaisMeio = return (Q(pos)-Q(pos-1))/(Q(pos+1)-Q(pos));
end

function calcPsi = calculaPsiSuperbee(theta)
    calcPsi = max(0, min(1, 2*theta), min(2, theta));
end

function calcPsi = calculaPsiMc(theta)
    calcPsi = max(0, min(((1+theta)/2), 2, 2*theta));
end

function calcPsi = calculaPsiVanleer(theta)
    calcPsi = (theta + abs(theta)) / (1 + abs(theta));
end

function executaMetodo(calculaPsi)
    I = 2000;
    L = 2;
    uBarra = 0.2;
    deltaT = 0.005;
    tempoTotal = 4;
    c0 = 0.4;
    C = 0.7;
    deltaX = L/I;
    courant = (uBarra * deltaT) / (deltaX);
    disp(courant);

    qAntigo = zeros(1, I + 2);
    qNovo = zeros(1, I + 2);

    i = 2;
    while(i <= I+1) then
        qAntigo(i) = c0;
        i = i + 1;
    end

    qAntigo(1) = C;
    qAntigo(I+2) = qAntigo(I-1);

    tempoAtual = 0;
    while(tempoAtual < tempoTotal) then
        tempoAtual = tempoAtual+deltaT;
        i = 2;
        while(i <= I + 1) then
            thetaIMenosMeio = calculaThetaRecuado(qAntigo, i);
            thetaIMaisMeio = calculaThetaAvancado(qAntigo, i);
            qNovo(i) = (qAntigo(i) - courant*(qAntigo(i)-qAntigo(i-1)))-0.5*courant*(1-courant)*(calculaPsi(thetaIMaisMeio)*(qAntigo(i+1)-qAntigo(i))-calculaPsi(thetaIMenosMeio)*(qAntigo(i)-                      qAntigo(i-1)));
            i = i + 1;
        end
        qAntigo(I+2) = qAntigo(I+1);
        i = 2;
        while(i <= I+1) then
            qAntigo(i) = qNovo(i);
            i = i + 1;
        end
    end

    x = linspace(0, I, I);
    qNovo = qNovo(1, 2:I+1);
    figure();
    plot2d(x, qNovo, rect=[0, 0, I, 1]);
    xlabel('Número de células');
    ylabel('Concentração do traçador');
endfunction

executaMetodo(calculaPsiSuperbee);

executaMetodo(calculaPsiMc);

executaMetodo(calculaPsiVanleer);