clc
clear all
close all

load('results_3P_Side_C.mat')
load('Match_Parameters_3P_Side_C.mat')

tic
Data_p = length(T); 

mwell_i = [1];

xg = xg';

count = 1;
TC = length(mwell_i)*length(mwell_p);


inj = 35000*ones(Data_p,length(mwell_i));
qp = -q_p(:,1:Data_p)';
cwf = cs_hist(:,1:Data_p)';
cinj = 350;
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)])

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)])


figure(4)
plot([T(1);T(1000:1000:end)],[cwf(1,4);cwf(1000:1000:end,4)],'o','LineWidth',1.5,'Color','#7E2F8E')
hold on
plot(T(1:5000),cwf_prod(1:5000,4),'--r','LineWidth',1.5)
hold on
plot(T(5001:end),cwf_prod(5001:end,4),'--','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(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);
Err4 = cwf_prod(:,4)-cwf(:,4);
Err1 = (sqrt(Err1'*Err1)/length(Err1));
Err2 = (sqrt(Err2'*Err2)/length(Err2));
Err3 = (sqrt(Err3'*Err3)/length(Err3));
Err4 = (sqrt(Err4'*Err4)/length(Err4));