clc;
clear all;
close all;
% Read exle file
filename = 'C:\PhD in MCM\Dofasco project\MeOH\lapse\sensitivityv1.xlsx';
Excel = actxserver('Excel.Application');
Workbooks = Excel.Workbooks;
Excel.Visible =1;

sheet = 9;
x1Range = 'D1:G6';
x2Range = 'N1:R7';
x3Range = 'D1:H1';

[A D both] = xlsread(filename,sheet,x1Range);
[B C both] = xlsread(filename,sheet,x2Range);   % B read numbers; C read text. both means read both
network_station_coordinates = both;
E = xlsread (filename,10,x3Range);
ECOG = 214288224;   % status-quo electricity generated from COG. kwh


% Define the input grid
[Xelec,Xmeoh] =meshgrid (0:50:300,200:80:680);



for j = 1:1:5           %AF factor
     
for i = 1:1:5       % Case study
 figure (i + (j-1)*5)  
% subplot(3,2,i)
zlim([-10 30]);           % Set limit for Z value

TCO21 = -0.8845*A(1,1)*Xelec/((A(5,1)+ECOG*B(6,i)/1E9)*B(5,i)) ...
    +0.8845*A(2,1)*Xmeoh/(A(5,1)+ECOG*B(6,i)/1E9)...
    -A(4,1)*1E3/(B(1,i)*(A(5,1)+ECOG*B(6,i)/1E9))...
    -A(3,1)*1E3/(B(1,i)*E(1,j)*(A(5,1)+ECOG*B(6,i)/1E9));

  
NPV1 = -0.8845*E(1,j)*A(1,1)*B(3,i)/(B(5,i)*1E3)...
    + 0.8845*E(1,j)*A(2,1)*B(4,i)/1E3...
    - E(1,j)*(A(5,1)+ECOG*B(6,i)/1E9)*B(2,i)/1E3...
    - E(1,j)*A(4,1)/B(1,i)...
    - A(3,1)/B(1,i)

TCO22 = -0.8845*A(1,2)*Xelec/((A(5,2)+ECOG*B(6,i)/1E9)*B(5,i)) ...
    +0.8845*A(2,2)*Xmeoh/(A(5,2)+ECOG*B(6,i)/1E9)...
    -A(4,2)*1E3/(B(1,i)*(A(5,2)+ECOG*B(6,i)/1E9))...
    -A(3,2)*1E3/(B(1,i)*E(1,j)*(A(5,2)+ECOG*B(6,i)/1E9));
   
NPV2 = -0.8845*E(1,j)*A(1,2)*B(3,i)/(B(5,i)*1E3)...
    + 0.8845*E(1,j)*A(2,2)*B(4,i)/1E3...
    - E(1,j)*(A(5,2)+ECOG*B(6,i)/1E9)*B(2,i)/1E3...
    - E(1,j)*A(4,2)/B(1,i)...
    - A(3,2)/B(1,i)
  
TCO23 = -0.8845*A(1,3)*Xelec/((A(5,3)+ECOG*B(6,i)/1E9)*B(5,i)) ...
    +0.8845*A(2,3)*Xmeoh/(A(5,3)+ECOG*B(6,i)/1E9)...
    -A(4,3)*1E3/(B(1,i)*(A(5,3)+ECOG*B(6,i)/1E9))...
    -A(3,3)*1E3/(B(1,i)*E(1,j)*(A(5,3)+ECOG*B(6,i)/1E9));

  
NPV3 = -0.8845*E(1,j)*A(1,3)*B(3,i)/(B(5,i)*1E3)...
    + 0.8845*E(1,j)*A(2,3)*B(4,i)/1E3...
    - E(1,j)*(A(5,3)+ECOG*B(6,i)/1E9)*B(2,i)/1E3...
    - E(1,j)*A(4,3)/B(1,i)...
    - A(3,3)/B(1,i)

TCO24 = -0.8845*A(1,4)*Xelec/((A(5,4)+ECOG*B(6,i)/1E9)*B(5,i)) ...
    +0.8845*A(2,4)*Xmeoh/(A(5,4)+ECOG*B(6,i)/1E9)...
    -A(4,4)*1E3/(B(1,i)*(A(5,4)+ECOG*B(6,i)/1E9))...
    -A(3,4)*1E3/(B(1,i)*E(1,j)*(A(5,4)+ECOG*B(6,i)/1E9));
   
NPV4 = -0.8845*E(1,j)*A(1,4)*B(3,i)/(B(5,i)*1E3)...
    + 0.8845*E(1,j)*A(2,4)*B(4,i)/1E3...
    - E(1,j)*(A(5,4)+ECOG*B(6,i)/1E9)*B(2,i)/1E3...
    - E(1,j)*A(4,4)/B(1,i)...
    - A(3,4)/B(1,i)


[NPV,I] = max ([NPV1 NPV2  NPV3  NPV4])     % Find the maximum NPV

z = -0.8845*A(1,I)*Xelec/((A(5,I)+ECOG*B(6,i)/1E9)*B(5,i)) ...
    +0.8845*A(2,I)*Xmeoh/(A(5,I)+ECOG*B(6,i)/1E9)...
    -A(4,I)*1E3/(B(1,i)*(A(5,I)+ECOG*B(6,i)/1E9))...
    -A(3,I)*1E3/(B(1,i)*E(1,j)*(A(5,I)+ECOG*B(6,i)/1E9));

PY =   (A(3,I)/B(1,i))/...
    ( 0.8845*A(2,I)*B(4,i)/1E3...
    - 0.8845*A(1,I)*B(3,i)/(B(5,i)*1E3)...
    - B(2,i)*(A(5,I)+ECOG*B(6,i)/1E9)/1E3...
    - A(4,I)/B(1,i)...
    )    % Payback period yr

CCA = NPV*1E3/((A(5,I)+ECOG*B(6,i)/1E9)*30)    % The best scenario CCA

% quiver3, get arrow/direction of the surface increasing 
[u,v,w]=surfnorm(Xelec, Xmeoh, z);
q = quiver3(Xelec, Xmeoh, z,-u,-v,-w,0.05)
hold on


% Visualize the two surfaces
h1 = surface(Xelec, Xmeoh, TCO21, 'FaceColor', 'y', 'EdgeColor', 'none')
alpha(0.5)

h2 = surface(Xelec, Xmeoh, TCO22, 'FaceColor', 'r', 'EdgeColor', 'none')
alpha(0.5)

h3 = surface(Xelec, Xmeoh, TCO23, 'FaceColor', 'm', 'EdgeColor', 'none')
alpha(0.5)

h4 = surface(Xelec, Xmeoh, TCO24, 'FaceColor', 'c', 'EdgeColor', 'none')
alpha(0.5)

view(35,-45)
hold on
 
% legend(D(1,2));

xlabel('\itX_{Elec.} ($/MWH)');
% xlabel('\itX_{Elec.} (E6 $/MWH)');
ylabel('\itX_{MeOH} ($/Mton)');
zlabel('CO_{2}\rmTax ($/Mton)');
title(C(1,i));
view(3); camlight; axis vis3d



% case study
hold on;
h5 = plot3(B(3,i),B(4,i),B(2,i), 'r*');     % Case Study
hold on;

% TCO24 = -0.8845*A(1,4)*100*Xelec/(A(5,4)*B(5,i)) ...
%     +0.8845*A(2,4)*Xmeoh/A(5,4)...
%     -A(4,4)*1E3/(B(1,i)*A(5,4))...
%     -A(3,4)*1E3/(B(1,i)*E(1,j)*A(5,4));
%    
% NPV4 = -0.8845*E(1,j)*A(1,4)*100*B(3,i)/(B(5,i)*1E3)...
%     + 0.8845*E(1,j)*A(2,4)*B(4,i)/1E3...
%     - E(1,j)*A(5,4)*B(2,i)/1E3...
%     - E(1,j)*A(4,4)/B(1,i)...
%     - A(3,4)/B(1,i)
  
dist = abs( -0.8845*A(1,I)*B(3,i)*E(1,j)/(B(5,i)*1E3)...
    + 0.8845*A(2,I)*B(4,i)*E(1,j)/1E3...
    - (A(5,I)+ECOG*B(6,i)/1E9)*B(2,i)*E(1,j)/1E3...
    -A(4,I)*E(1,j)/B(1,i)...
    -A(3,I)/B(1,i))/...
    sqrt( (-0.8845*A(1,I)*E(1,j)/(B(5,i)*1E3))^(2)...
    + (0.8845*A(2,I)*E(1,j)/1E3)^(2)...
    + (-(A(5,I)+ECOG*B(6,i)/1E9)*E(1,j)/1E3)^(2));            % Point to plane distance




% triggle point of carbon tax
X3 = B(3,i);        %X value of a point in the plane: electricity price
Y3 = B(4,i);        %Y value of a point in the plane: MeOH price
f3 = @(Z3)  -0.8845*A(1,I)*X3/((A(5,I)+ECOG*B(6,i)/1E9)*B(5,i)) ...
    +0.8845*A(2,I)*Y3/(A(5,I)+ECOG*B(6,i)/1E9)...
    -A(4,I)*1E3/(B(1,i)*(A(5,I)+ECOG*B(6,i)/1E9))...
    -A(3,I)*1E3/(B(1,i)*E(1,j)*(A(5,I)+ECOG*B(6,i)/1E9))...
    - Z3;

z3 = fzero(f3,1);          % Z value of a point in the plane: that is the CO2 tax triger point

% triggle point of elect price.
Z1 = B(2,i);        %Z value of a point in the plane
Y1 = B(4,i);        %Y value of a point in the plane
f1 = @(X1) A(2,I)*B(5,i)*Y1/(A(1,I))...
    - (A(5,I)+ECOG*B(6,i)/1E9)*B(5,i)*Z1/(0.8845*A(1,I))...
    - A(4,I)*B(5,i)*1E3/(0.8845*B(1,i)*A(1,I))...
    - A(3,I)*B(5,i)*1E3/(0.8845*B(1,i)*E(1,j)*A(1,I))...
    - X1;

x1 = fzero(f1,1)          

% triggle point of MeOH price.
Z2 = B(2,i);        %Z value of a point in the plane
X2 = B(3,i);        %X value of a point in the plane
f2 = @(Y2) A(1,I)*X2/(B(5,i)*A(2,I))...
    + (A(5,I)+ECOG*B(6,i)/1E9)*Z2/(0.8845*A(2,I))...
    + A(4,I)*1E3/(0.8845*B(1,i)*A(2,I))...
    + A(3,I)*1E3/(0.8845*B(1,i)*A(2,I)*E(1,j))...    
    - Y2;

y2 = fzero(f2,1);    


tp = [x1;y2;z3]                    %triggle point

p1 = [x1; Y1; Z1]                 % electricity price trigger point on the plan
p3 = [X2;y2,;Z2]                    % MeOH price triggle point

% Normal vector of the plane is [a,b,c]. where ax+by+cz=d.

N = [-0.8845*A(1,I)*E(1,j)/(B(5,i)*1E3); 0.8845*A(2,I)*E(1,j)/1E3; -(A(5,I)+ECOG*B(6,i)/1E9)*E(1,j)/1E3];

p2 = [B(3,i); B(4,i); B(2,i)]     % Cast study point
% http://www.nabla.hr/CG-LinesPlanesIn3DB5.htm

lamda = ( E(1,j)*A(4,I)/B(1,i) + A(3,I)/B(1,i) -N'*p2 )/(N'*N)
   
% Projection point on the plane: 
px= lamda*N(1,1)+p2(1,1); 
py= lamda*N(2,1)+p2(2,1); 
pz= lamda*N(3,1)+p2(3,1);
p= [px;py;pz]

NPVa = -0.8845*E(1,j)*A(1,I)*p(1,1)/(B(5,i)*1E3)...
    + 0.8845*E(1,j)*A(2,I)*p(2,1)/1E3...
    - E(1,j)*(A(5,I)+ECOG*B(6,i)/1E9)*p(3,1)/1E3...
    - E(1,j)*A(4,I)/B(1,i)...
    - A(3,I)/B(1,i)            %Check if the plane means NPV = 0.
   
% check if p-p2 is perpendicular to p1-p3
g = (p-p2)'*(p1-p3)

% xlswrite(filename,A,sheet,xlRange) writes to the specified worksheet and range.
  
    if j == 1
        startingColumn = 'D';
        
    elseif j==2
        startingColumn = 'E';
%         nextColumn = 2 + 4*(i-1);
%         range = sprintf('%s%d',startingColumn,nextColumn);
%         xlswrite(filename,NPV,10,range);
    elseif j==3
        startingColumn = 'F';  
%         nextColumn = 2 + 4*(i-1);
%         range = sprintf('%s%d',startingColumn,nextColumn);
%         xlswrite(filename,NPV,10,range);
    elseif j==4
        startingColumn = 'G';  
%         nextColumn = 2 + 4*(i-1);
%         range = sprintf('%s%d',startingColumn,nextColumn);
%         xlswrite(filename,NPV,10,range);
    else 
        startingColumn = 'H';  
%         nextColumn = 2 + 4*(i-1);
%         range = sprintf('%s%d',startingColumn,nextColumn);
%         xlswrite(filename,NPV,10,range);
    end
        nextrownpv = 2 + 6*(i-1);
        rangenpv = sprintf('%s%d',startingColumn,nextrownpv);
        xlswrite(filename,round(NPV),10,rangenpv);

% Write sensitivity analysis result into Excel
        tN = sprintf(' %.4f ',N);      % Put all numbers into a string   
        tN = sprintf('[%s]',tN);    % Enclose in square brackets
        strN = {tN};
        nextrownorm = 3 + 6*(i-1);
        rangenorm = sprintf('%s%d',startingColumn,nextrownorm);
        xlswrite(filename,strN,10,rangenorm);

% Triggle point
% 
        ttp = sprintf(' %.2f ',tp);      % Pua all numbers into a string; keep two decimals   
        ttp = sprintf('[%s]',ttp);    % Enclose in square brackets
        strtp = {ttp};

        nextrowtp = 4 + 6*(i-1);           
        rangetp = sprintf('%s%d',startingColumn,nextrowtp);
        xlswrite(filename,strtp,10,rangetp);

% payback period
        nextrowy = 5 + 6*(i-1);
        rangey = sprintf('%s%d',startingColumn,nextrowy);
        xlswrite(filename,round(PY),10,rangey);
% Best Process
        nextrowI = 6 + 6*(i-1);
        rangeI = sprintf('%s%d',startingColumn,nextrowI);
        xlswrite(filename,I,10,rangeI);
% CCA $/Mton CO2e
        nextrowcca = 7 + 6*(i-1);
        rangecca = sprintf('%s%d',startingColumn,nextrowcca);
        xlswrite(filename,CCA,10,rangecca);


  
  
h6 = plot3 (p(1,1), p(2,1), p(3,1), 'k*')          %Check the position of the assumed point in plane and the projection point on the plane
hold on;
h7 = plot3 ([p2(1,1) p(1,1)],[p2(2,1) p(2,1)],[p2(3,1) p(3,1)],'k');     %Draw the line according to the case study point and its projection on the plane

hold on;
str1 = 'NPV =';
str2 = num2str(round(NPV));
str3 = '\times10^{6}$';
% str3 = 'E6 $';
str4 = 'Y= ';
str5 = num2str(round(PY*100)/100);
str6 = ' year'
stra = [str1 str2 str3];
strb = [str4 str5 str6];
str = {stra,strb};
% annotation. Put the distance to the perpendicular line
k = [p(1,1) B(3,i)];
m = [p(2,1) B(4,i)];
n = [p(3,1) B(2,i)];
text(k(end),m(end),n(end),str);

legend([q,h1,h2,h3,h4,h5,h6,h7], {'NPV Increasing Direction','uCOG','uNG','uBFG','uNG-WGS','Case Study','Projection', 'Distance'});

grid on


end

end







