% Coded by NFR on 11.7.17, updated on 4.5.2019 
% Solves DOE Problem
%Part1
CCIval = ccdesign(2,'type','inscribed');
CCCval = ccdesign(2,'type','circumscribed');
m = 16;
%Randomize run order and set up real world limits
runorder = randperm(16);
bounds = [4 80; 2 10];
RV1 = zeros(size(CCIval));
RV2 = zeros(size(CCCval));
for i = 1:2
    zmax = max(CCIval(:,i));
    zmin = min(CCIval(:,i));
    RV1(:,i) = interp1([zmin zmax],bounds(i,:),CCIval(:,i));
    %zmax = max(CCCval(:,i));
    %zmin = min(CCCval(:,i));
    % The CCC design should go beyond the design bounds, keep zmax and zmin
    % at 1 and -1 respectively.
    RV2(:,i) = interp1([zmin zmax],bounds(i,:),CCCval(:,i));
    % Looks like the code we used before DOES NOT do a good job showing the
    % axial points.  For the sake of time, I will adjust these manually.
    AxialD1 = 1.4142*(bounds(1,2)-bounds(1,1))/2;
    AxialD2 = 1.4142*(bounds(2,2)-bounds(2,1))/2;
    RV2(5,1) = (bounds(1,2)-bounds(1,1))/2-AxialD1;
    RV2(6,1) = AxialD1+(bounds(1,2)-bounds(1,1))/2;
    RV2(7,2) = (bounds(2,2)-bounds(2,1))/2-AxialD2;
    RV2(8,2) = AxialD2+(bounds(2,2)-bounds(2,1))/2;
end
disp('########## Part 1 - CCI Design ###########')  % NOTE to grader - FOr part 1 please focus only on the CCI design.  THe CCC design went beyond scope of the problem.
disp({'Run Number','Pressure','Ratio'})
disp(sortrows([runorder' RV1]))
disp('########## Part 1 - CCC Design ###########')
disp({'Run Number','Pressure','Ratio'})
disp(sortrows([runorder' RV2]))
%NOTE: At these bounds you cannot use the CCC design (the low axial points for
%pressure and ratio go into negative values.  Stick with the CCI design.
%
% Part 2
D = xlsread('DataFormatted.xlsx');
TR1 = D(:,4);
TR2 = D(:,5);
RealValue = D(:,2:3);
% Fit to both outputs to quadratic models
mdl_uniform = fitlm(RealValue,TR1,'quadratic') % Don't suppress and you get model parameters
mdl_stress = fitlm(RealValue,TR2,'quadratic')
% Look at p values to determine dominant terms.  From the model fit, we
% will keep the first FOUR parameters for the uniformity response model and
% we will keep 1-3,5 parameters for the stressresponse model.  
%
% Remove the term that has the largest p-value
New_mdl_uniform = removeTerms(mdl_uniform,'x2:x2')
New_mdl_stress = removeTerms(mdl_stress,'x1:x2')
% Inspect again and remove any other term that is still above p-pvalue
% threshold
New_mdl_uniform_2 = removeTerms(New_mdl_uniform,'x1:x1')
New_mdl_stress_2 = removeTerms(New_mdl_stress,'x2:x2')
%
% Here are the equations (where x1 is pressure and x2 is ratio of feed gases):
Pu = New_mdl_uniform_2.Coefficients.Estimate;
Ps = New_mdl_stress_2.Coefficients.Estimate;
% The reduced model fits:
f_uniform = @(x1,x2) Pu(1)+Pu(2).*x1+Pu(3).*x2+Pu(4).*x1.*x2;
f_stress = @(x1,x2) Ps(1)+Ps(2).*x1+Ps(3).*x2+Ps(4).*x1.^2;
%
% Use fsolve to find the conditions settings for these two desired output
% parameters
xo = [50,5];
Solution = fsolve(@func1,xo,[],Pu,Ps);
disp('Settings for uniformity = 5.7 and stress = 7.7:')
Pressure_setting = Solution(1)
Ratio_Setting = Solution(2)
% Plot them!
xvec = linspace(4,80,100);
yvec = linspace(2,10,100);
[xmat,ymat] = meshgrid(xvec,yvec);
z_u = f_uniform(xmat,ymat);
z_s = f_stress(xmat,ymat);
[C,h] = contour(xvec,yvec,z_u);
clabel(C,h)
hold on
[C,h] = contour(xvec,yvec,z_s,'--');
clabel(C,h)
ylabel('Ratio H2/WF6')
xlabel('Pressure (Torr)')
title('Solid Line = Uniformity, Dashed Line = Stress, o = set point')
plot(Solution(1),Solution(2),'o')
hold off

function F = func1(X,Pu,Ps)
x1 = X(1);
x2 = X(2);
% Turn the model functions into a root solving problem, by subtracting the
% optimal point:
F(1) = Pu(1)+Pu(2)*x1+Pu(3)*x2+Pu(4)*x1*x2-5.7;
F(2) = Ps(1)+Ps(2)*x1+Ps(3)*x2+Ps(4)*x1^2-7.7;
end