clc
clear all
close all

Results = xlsread('Results_Case_CMG_1.xlsx');
load('Match_Parameters_3I2P_CMG_3.mat')

tic
inj = Results(:,2:4);
qp = Results(:,5:6);
cwf = Results(:,7:8);
cinj = Results(:,9);
T = Results(:,1);
ni = 3; 
np = 2; 
TC = ni*np;

Data_p = ceil((1/3)*size(Results,1));

xg = xg';

count = 1;
TC = ni*np;
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

NDiv = round(xg(4*TC+1:end));

hi = [1153565.888683;1167041.425060;1165981.176032];
ii = [-24665171.424087;-24657424.130630;-24664917.158948];
hp = [1160762.307872;1164372.385021];
ip = [-24663256.643936;-24660351.459976];

Lwell = CalculateLength(hi,ii,hp,ip);

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:Data_p),cwf_prod(1:Data_p,1),'--r','LineWidth',1.5)
hold on
plot(T(Data_p+1:end),cwf_prod(Data_p+1: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)])

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:Data_p),cwf_prod(1:Data_p,2),'--r','LineWidth',1.5)
hold on
plot(T(Data_p+1:end),cwf_prod(Data_p+1: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(4)
plot(T,inj(:,1),'--b','LineWidth',1.5)
hold on
plot(T,inj(:,2),'-.k','LineWidth',1.5)
hold on
plot(T,inj(:,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)
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');
set(l,'location','southeast')

figure(6)
plot(T,cinj,'--b','LineWidth',1.5)
xlabel('Time (Days)','FontSize',12,'FontWeight','bold')
ylabel('Injected Contaminant Concentration (g/L)','FontSize',12,'FontWeight','bold')
set(gca,'FontSize',12,'LineWidth',1.5,'FontWeight','bold') 
ylim([100 400])

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