clc
clear all
close all

load('results_3I3P_Large.mat')
load('Match_Parameters_3I3P_Large_2.mat')

tic
Data_p = length(T); 

mwell_i = sort(mwell_i);

xg = xg';

count = 1;
TC = length(mwell_i)*length(mwell_p);
for L = 1:TC
    xg1(count,1) = xg(L);
    xg1(count+1,1) = xg(L+TC);
    xg1(count+2,1) = xg(L+TC+TC);
    xg1(count+3,1) = xg(L+TC+TC+TC);
    xg1(count+4,1) = xg(L+TC+TC+TC+TC);
    count = count+5;
end


for t = 1:length(T)
    if t<=1000
        cinj_p(t,1) = 250;
    elseif t>1000 && t<=2000
        cinj_p(t,1) = 200;
    elseif t>2000 && t<=2500
        cinj_p(t,1) = 350;
    elseif t>2500 && t<=3500
        cinj_p(t,1) = 170;
    elseif t>3500 && t<=5000
        cinj_p(t,1) = 200;
    elseif t>5000 && t<=7500
        cinj_p(t,1) = 250;
    elseif t>7500 && t<=10000
        cinj_p(t,1) = 350;
    elseif t>10000 && t<=15000
        cinj_p(t,1) = 220;
    elseif t>15000 && t<=18000
        cinj_p(t,1) = 300;
    elseif t>18000
        cinj_p(t,1) = 265;
    end
end


inj = q_i(:,1:Data_p)';
qp = -q_p(:,1:Data_p)';
cwf = cs_hist(:,1:Data_p)';
cinj = 250;
T = T(1:Data_p)';
ni = length(mwell_i);
np = length(mwell_p);
TC = ni*np;

NDiv = round(xg(4*TC+1:end));
Lwell = CalculateLength(mwell_i,mwell_p,H,I,dx,dy);

count = 0;
cwf_prod = zeros(length(T),np);
for ii = 1:ni
    len1 = np*(ii-1)+1;
    len2 = np*ii;
    fij = xg(len1:len2); %fij for corresponding injector with all producer
    for pp = 1:np
        sum1 = 0;
        f = xg(pp:np:TC); %fij for corresponding producer with all injector
        for i1 = 1:ni
            sum1 = sum1 + (f(i1)*inj(:,i1));
        end
        qp1 = (qp(:,pp)*f(ii).*inj(:,ii))./sum1;
        inj1 = inj(:,ii)*fij(pp);
        len3 = TC + ((ii-1)*np) + pp;
        len4 = 2*TC + ((ii-1)*np) + pp;
        len5 = 3*TC + ((ii-1)*np) + pp;
        Vp = xg(len3);
        D = xg(len4);
        qf = xg(len5);
        cwf_prod(:,pp) = cwf_prod(:,pp) + ADE_model(qp1,inj1,T,cinj,NDiv(count+1),Lwell(count+1),Vp,D,qf);
        count = count + 1;
    end
end
toc

%% Plots
figure(1)
plot([T(1);T(1000:1000:end)],[cwf(1,1);cwf(1000:1000:end,1)],'o','LineWidth',1.5,'Color','#7E2F8E')
hold on
plot(T(1:5000),cwf_prod(1:5000,1),'--r','LineWidth',1.5)
hold on
plot(T(5001:end),cwf_prod(5001:end,1),'--','LineWidth',1.5,'Color','#0072BD')
xlabel('Time (Days)','FontSize',12,'FontWeight','bold')
ylabel('Concentration (g/L)','FontSize',12,'FontWeight','bold')
set(gca,'FontSize',12,'LineWidth',1.5,'FontWeight','bold') %'FontWeight','bold',

l=legend('Full Scale Model','ISADE Model - History Matched','ISADE Model - Predicted');
set(l,'location','southeast')
xlim([T(1) T(end)])
% ylim([0 400])

figure(2)
plot([T(1);T(1000:1000:end)],[cwf(1,2);cwf(1000:1000:end,2)],'o','LineWidth',1.5,'Color','#7E2F8E')
hold on
plot(T(1:5000),cwf_prod(1:5000,2),'--r','LineWidth',1.5)
hold on
plot(T(5001:end),cwf_prod(5001:end,2),'--','LineWidth',1.5,'Color','#0072BD')
xlabel('Time (Days)','FontSize',12,'FontWeight','bold')
ylabel('Concentration (g/L)','FontSize',12,'FontWeight','bold')
set(gca,'FontSize',12,'LineWidth',1.5,'FontWeight','bold') %'FontWeight','bold',

l=legend('Full Scale Model','ISADE Model - History Matched','ISADE Model - Predicted');
set(l,'location','southeast')
xlim([T(1) T(end)])

figure(3)
plot([T(1);T(1000:1000:end)],[cwf(1,3);cwf(1000:1000:end,3)],'o','LineWidth',1.5,'Color','#7E2F8E')
hold on
plot(T(1:5000),cwf_prod(1:5000,3),'--r','LineWidth',1.5)
hold on
plot(T(5001:end),cwf_prod(5001:end,3),'--','LineWidth',1.5,'Color','#0072BD')
xlabel('Time (Days)','FontSize',12,'FontWeight','bold')
ylabel('Concentration (g/L)','FontSize',12,'FontWeight','bold')
set(gca,'FontSize',12,'LineWidth',1.5,'FontWeight','bold') %'FontWeight','bold',

l=legend('Full Scale Model','ISADE Model - History Matched','ISADE Model - Predicted');
set(l,'location','southeast')
xlim([T(1) T(end)])


T1 = [T(20:1999);T(2020:4999);T(5020:9999);T(10020:end)];
inj1 = [inj(20:1999,:);inj(2020:4999,:);inj(5020:9999,:);inj(10020:end,:)];

figure(4)
plot(T1,inj1(:,1),'-b','LineWidth',1.5)
hold on
plot(T1,inj1(:,2),'-.k','LineWidth',1.5)
hold on
plot(T1,inj1(:,3),':','color',[0.8500 0.3250 0.0980],'LineWidth',1.5)
xlabel('Time (Days)','FontSize',12,'FontWeight','bold')
ylabel('Injection Rate (bbls/day)','FontSize',12,'FontWeight','bold')
set(gca,'FontSize',12,'LineWidth',1.5,'FontWeight','bold') %'FontWeight','bold',
xlim([T(1) T(end)])
% ylim([1000 8500])
l=legend('I1','I2','I3');
set(l,'location','northeast')


figure(5)
plot(T(3:end),qp(3:end,1),'-','color',[0.4660 0.6740 0.1880],'LineWidth',1.5)
hold on
plot(T(3:end),qp(3:end,2),'-','color',[0.3010 0.7450 0.9330],'LineWidth',1.5)
hold on
plot(T(3:end),qp(3:end,2),':r','LineWidth',1.5)
xlabel('Time (Days)','FontSize',12,'FontWeight','bold')
ylabel('Production Rate (bbls/day)','FontSize',12,'FontWeight','bold')
set(gca,'FontSize',12,'LineWidth',1.5,'FontWeight','bold') %'FontWeight','bold',
xlim([T(1) T(end)])
l=legend('P1','P2','P3');
set(l,'location','southeast')

Err1 = cwf_prod(:,1)-cwf(:,1);
Err2 = cwf_prod(:,2)-cwf(:,2);
Err3 = cwf_prod(:,3)-cwf(:,3);
Err1 = (sqrt(Err1'*Err1)/length(Err1));
Err2 = (sqrt(Err2'*Err2)/length(Err2));
Err3 = (sqrt(Err3'*Err3)/length(Err3));