%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Constrained Least Squares Formulation %%%%
%%%%%%   For Thermo Properties %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% By Jack Ziegler Jan. 2008
% jackalak@caltech.edu
%
% Ax = b --> least squares minumization
% Aeq*x = beq --> constraint equations
% x = [a1; a2]
% a1 = [a10; a11; ...; a16]
% a2 = [a20; a21; ...; a26];
% Solved using Matlab's lsqlin 
% iterative constrained least squares function .
% (internally uses Lagrange multipliers)
%
% without the constraints the form of the linear least squares
% minumization is A*x = b where A is mxn, x is nx1, b is mx1
% would be solved using (A' A) * x = (A' b)

function interp_thermo(substance, composition, Cantera_output_text, DATA, Tmin, Tmid, Tmax, so_over_R_298,ho_over_RT_298)

T = [Tmin:100:Tmax]';
index_mid = (Tmid-Tmin)/100+1;
index_max = (Tmax-Tmin)/100+1;

% **** THERMOCHEMICAL PROPERTIES ****     (PRESSURE = 1.0000 atm)
cp_over_R = DATA(:,3)/1.987; % R = 1.987 cal K-1 mol-1
h_over_RT =  DATA(:,5)/1.987*1000 ./ T;
s0_over_R = DATA(:,4)/1.987;


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Entropy and enthapy of formation matching
% %%% Now fix the so^f and ho^f values so that they match those at 298K = ~300K
h_over_R_correction =(ho_over_RT_298 - h_over_RT((300-Tmin)/100+1))*300;
h_over_RT = (DATA(:,5)/1.987*1000 +h_over_R_correction )./ T;
s0_over_R =(so_over_R_298 - s0_over_R((300-Tmin)/100+1))+s0_over_R;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% RHS
b = [ cp_over_R(1:index_mid); cp_over_R(index_mid:index_max)
    h_over_RT(1:index_mid); h_over_RT(index_mid:index_max)
    s0_over_R(1:index_mid);  s0_over_R(index_mid:index_max) ];

%Zero out A
A = zeros(((index_max+1)*3),2*7);  % index_max+1 data points for each thermo property (one middle point is shared) 
% A --> Organization %
% cp_over_R first range
% cp_over_R second range
% h_over_RT first range
% h_over_RT second range
% s0_over_R first range
% s0_over_R second range

% cp_over_R
for i = 1:index_mid
    for j = 0:4
        A(i,j+1) = T(i).^(j);
    end
end
for i = index_mid:index_max
    for j = 0:4
        A(i+1,j+1+7) = T(i).^(j);
    end
end

% h_over_RT
for i = 1:index_mid
    for j = 0:4
        A(i+index_max+1,j+1) = 1/(j+1)*T(i).^(j);
    end
        A(i+index_max+1,6) = 1./T(i);
end
for i = index_mid:index_max
    for j = 0:4
        A(i+index_max+1+1,j+1+7) = 1/(j+1)*T(i).^(j);
    end
        A(i+index_max+1+1,6+7) = 1./T(i);
end

% s0_over_R
for i = 1:index_mid  
    A(i+(index_max+1)*2,1) = log(T(i));
    for j = 1:4
        A(i+(index_max+1)*2,j+1) = 1/(j)*T(i).^(j);
    end
        A(i+(index_max+1)*2,7) = 1;
end
for i = index_mid:index_max
     A(i+(index_max+1)*2+1,1+7) = log(T(i));
    for j = 1:4
        A(i+(index_max+1)*2+1,j+1+7) = 1/(j)*T(i).^(j);
    end
        A(i+(index_max+1)*2+1,7+7) = 1;
end

% This is what one could use if constraints were not considered
%x = (A'*A)\(A'*b)

%%%% constraint equation Aeq*x = beq %%%%
Aeq = zeros(15,14);
beq = zeros(15,1);

% C^0 at end points
beq(1) = cp_over_R(1);
beq(2) = cp_over_R(index_max);
beq(3) =h_over_RT(1);
beq(4) = h_over_RT(index_max);
%%% Below not used because of dependences
% beq(5) = s0_over_R(1);
% beq(6) = s0_over_R(index_max);


%%% cp_over_R %%%
% left endpoint
    for j = 0:4
        Aeq(1,j+1) = T(1).^(j);
    end
% right endpoint
    for j = 0:4
        Aeq(2,j+1+7) = T(index_max).^(j);
    end    
%continuous C^0 at knot
    for j = 0:4
        Aeq(7,j+1) = T(index_mid).^(j);
    end
    for j = 0:4
        Aeq(7,j+1+7) = -T(index_mid).^(j);
    end    
    %1st derivative C^1 at knot
    for j = 1:4
        Aeq(8,j+1) = j*T(index_mid).^(j-1);
    end
    for j = 1:4
        Aeq(8,j+1+7) = -j*T(index_mid).^(j-1);
    end
    % Below causes dependences in constraints
        %2nd derivative C^2 at knot
%     for j = 2:4
%         Aeq(9,j+1) = j*(j-1)*T(index_mid).^(j-2);
%     end
%     for j = 2:4
%         Aeq(9,j+1+7) = -j*(j-1)*T(index_mid).^(j-2);
%     end
%     
%%% h_over_RT %%%
% left endpoint
    for j = 0:4
        Aeq(3,j+1) = 1/(j+1)*T(1).^(j);
    end
        Aeq(3,6) = 1./T(1);
% right endpoint        
    for j = 0:4
        Aeq(4,j+1+7) = 1/(j+1)*T(index_max).^(j);
    end
        Aeq(4,6+7) = 1./T(index_max);
 % C^0 at Knot
% continuous
    for j = 0:4
        Aeq(10,j+1) = 1/(j+1)*T(index_mid).^(j);
    end
        Aeq(10,6) = 1./T(index_mid);
    for j = 0:4
        Aeq(10,j+1+7) = -1/(j+1)*T(index_mid).^(j);
    end
        Aeq(10,6+7) = -1./T(index_mid);        
% Below causes dependences in constraints
%     %C^1 at knot
%     for j = 1:4
%         Aeq(11,j+1) = j/(j+1)*T(index_mid).^(j-1);
%     end
%         Aeq(11,6) = -1./T(index_mid)^2;
%     for j = 1:4
%         Aeq(11,j+1+7) = -j/(j+1)*T(index_mid).^(j-1);
%     end
%         Aeq(11,6+7) = 1./T(index_mid)^2;
%  %C^2 at knot
%     for j = 1:4
%         Aeq(12,j+1) = (j-1)*j/(j+1)*T(index_mid).^(j-2);
%     end
%         Aeq(12,6) = 2./T(index_mid)^3;
%     for j = 1:4
%         Aeq(12,j+1+7) =-(j-1)*j/(j+1)*T(index_mid).^(j-2);
%     end
%         Aeq(12,6+7) = -2./T(index_mid)^3;
    
%%% s0_over_R %%%
    A(5,1) = log(T(1));
    for j = 1:4
        A(5,j+1) = 1/(j)*T(1).^(j);
    end
        A(5,7) = 1;
     Aeq(6,1+7) = log(T(index_max));
    for j = 1:4
        Aeq(6,j+1+7) = 1/(j)*T(index_max).^(j);
    end
        Aeq(6,7+7) = 1;        
   %C^0 at knot
    Aeq(13,1) = log(T(index_mid));
    for j = 1:4
        Aeq(13,j+1) = 1/(j)*T(index_mid).^(j);
    end
        Aeq(13,7) = 1;
     Aeq(13,1+7) = -log(T(index_mid));
    for j = 1:4
        Aeq(13,j+1+7) = -1/(j)*T(index_mid).^(j);
    end
        Aeq(13,7+7) = -1;
% Below causes dependences in constraints for s0/R
% 1st derivative
% 2nd derivative 

%%% Set up constraint Equations %%%
% Note: not all possible constraints are used

% C^1 Continuity for cp/R, C^0 for h/R/T and s0/R
% Also beginning and end boundary point constraints for cp/R and h/R/T
Aeq = [Aeq(1:4,:);   Aeq(7:8,:);  Aeq(10,:); Aeq(13,:)];
beq = [beq(1:4,:);    beq(7:8,:); beq(10,:); beq(13,:)];

%%% Alternate Simplier constraint choice
% % C^1 Continuity for cp/R, C^0 for h/R/T and s0/R
% Aeq = [ Aeq(7:8,:);  Aeq(10,:); Aeq(13,:)];
% beq = [ beq(7:8,:); beq(10,:); beq(13,:)];  



%%% perform the matlab constrained least squares optimization
%%% Note: For some reason matlab is only using one iteration for the
% medium scale algorithm (large scale is only for constraints with greater
% or less than or equals rather than equals (our case))
%%% To deal with this, the lsqlin function is called multiple times using
%%% the previous solution as an inital guess!  Check the error at the
%%% matching points in the end to make sure enough iterations have been
%%% carried out.

options = optimset('LargeScale','off','Display', 'off');
lb = ones(14,1)*-Inf; % lower bound for x
ub = ones(14,1)*Inf; % upper bound for x

% First iteration let matlab guess x0
[x,resnorm,residual,exitflag,output,lambda]  = lsqlin(A,b,[ ],[ ],Aeq,beq,lb,ub,[ ],options) ;
% Now call 10 times more using the previous x solution as the guess "x0"
for i = 1:10
    xo = x;
    [x,resnorm,residual,exitflag,output,lambda]  = lsqlin(A,b,[ ],[ ],Aeq,beq,lb,ub,xo,options) ;
end
% Now Turn the output options on and call once more to display the
% diagnotistics
options = optimset('LargeScale','off','Display', 'on','Diagnostics', 'on');
[x,resnorm,residual,exitflag,output,lambda]  = lsqlin(A,b,[ ],[ ],Aeq,beq,lb,ub,xo,options) ;
    

% first "polynomial"
a1 = x(1:7);
% second "polynomial"
a2 = x(8:14);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%% Evaluate the "polynomials" for plotting %%%%
T1= [Tmin:100:Tmid];
for j = 1:index_mid
    cp_over_R_1(j) = 0;
    for i = 1:5
        cp_over_R_1(j) = a1(i).*T1(j).^(i-1) +cp_over_R_1(j);
    end
end
T2= [Tmid:100:Tmax];
for j = 1:index_max-index_mid+1
    cp_over_R_2(j) = 0;
    for i = 1:5
        cp_over_R_2(j) = a2(i)*T2(j)^(i-1) +cp_over_R_2(j);
    end
end
for j = 1:index_mid
    h_over_RT_1(j) = 0;
    for i = 1:5
        h_over_RT_1(j) = 1/i*a1(i)*T1(j)^(i-1) +h_over_RT_1(j);
    end
    h_over_RT_1(j) = a1(6)/T1(j)+h_over_RT_1(j);
end
for j = 1:index_max-index_mid+1
    h_over_RT_2(j) = 0;
    for i = 1:5
        h_over_RT_2(j) = 1/i*a2(i)*T2(j)^(i-1) +h_over_RT_2(j);
    end
    h_over_RT_2(j) = a2(6)/T2(j)+h_over_RT_2(j);
end

for j = 1:index_mid
    s0_over_R_1(j) = a1(1)*log(T1(j));
    for i = 2:5
        s0_over_R_1(j) = 1/(i-1)*a1(i)*T1(j)^(i-1) +s0_over_R_1(j);
    end
    s0_over_R_1(j) = a1(7)+s0_over_R_1(j);
end
for j = 1:index_max-index_mid+1
    s0_over_R_2(j) = a2(1)*log(T2(j));
    for i = 2:5
        s0_over_R_2(j) = 1/(i-1)*a2(i)*T2(j)^(i-1) +s0_over_R_2(j);
    end
    s0_over_R_2(j) = a2(7)+s0_over_R_2(j);
end


%%%%% Plots %%%%%%
figure
subplot(2,2,1)
hold on
plot(T1,cp_over_R_1,'b');
plot(T2,cp_over_R_2,'b');
plot(T,cp_over_R,'.r');
hold off
title('cp/R')
subplot(2,2,3)
hold on
plot(T1,h_over_RT_1,'b');
plot(T2,h_over_RT_2,'b');
plot(T,h_over_RT,'.r');
hold off
title('h/RT')
subplot(2,2,4)
hold on
plot(T1,s0_over_R_1,'b');
plot(T2,s0_over_R_2,'b');
plot(T,s0_over_R,'.r');
hold off
title('s0/R (pressure independent)')
%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%
%%%% Output data and Diagonostics %%%%
%%%%%%%%%%%%%%%%%%%%%%%
standard_deviation_cp_over_R = std([cp_over_R_1'; cp_over_R_2(2:index_max-index_mid+1)']-cp_over_R);
standard_deviation_h_over_RT = std([h_over_RT_1'; h_over_RT_2(2:index_max-index_mid+1)']-h_over_RT);
standard_deviation_s0_over_R = std([s0_over_R_1'; s0_over_R_2(2:index_max-index_mid+1)']-s0_over_R);
disp(['standard_deviation_cp_over_R : ',num2str(standard_deviation_cp_over_R)]);
disp(['standard_deviation_h_over_RT : ',num2str(standard_deviation_h_over_RT)]);
disp(['standard_deviation_s0_over_R : ',num2str(standard_deviation_s0_over_R)]);
Norm_of_Residual = resnorm;
disp(['Norm_of_Residual ',num2str(Norm_of_Residual)]);

%%% Chi Squared goodness of fit data
[h_cp,p_cp,Chi_stats0_cp_over_R] = chi2gof([cp_over_R_1'; cp_over_R_2(2:index_max-index_mid+1)']-cp_over_R);
[h_h,p_h,Chi_stats0_h_over_RT] = chi2gof([h_over_RT_1'; h_over_RT_2(2:index_max-index_mid+1)']-h_over_RT);
[h_s,p_s,Chi_stats0_s0_over_R] = chi2gof([s0_over_R_1'; s0_over_R_2(2:index_max-index_mid+1)']-s0_over_R);

disp(['Chi_statistics0_cp_over_R :']);
Chi_stats0_cp_over_R
disp(['Chi_statistics0_h_over_RT :']);
Chi_stats0_h_over_RT
disp(['Chi_statistics0_s0_over_R :']);
Chi_stats0_s0_over_R

disp(['Min_Temperature : ',num2str(Tmin)]);
disp(['Matched_Temperature : ',num2str(Tmid)]);
disp(['Max_Temperature : ',num2str(Tmax)]);

disp(['Coefficients from Tmin to Tmid']);
for i = 1:7
    disp(['n = ', num2str(i-1), ',  ', num2str(x(i),16)]);
end
disp(['Coefficients from Tmid to Tmax']);
for i = 8:14
    disp(['n = ', num2str(i-1-7), ',  ', num2str(x(i),16)]);
end


disp(['Error_at_cp_over_R_match : ',num2str(cp_over_R_1((index_mid))-cp_over_R_2(1))]);
disp(['Error_at_h_over_RT_match : ',num2str(h_over_RT_1((index_mid))-h_over_RT_2(1))]);
disp(['Error_at_s0_over_R_match : ',num2str(s0_over_R_1((index_mid))-s0_over_R_2(1))]);

name_field = substance;
name_char_number = size(substance,2);
number_spaces = 18 - name_char_number;
for i = 1:number_spaces
    name_field =[name_field, ' '];
end

%Chemkin format output
outputfilename =[substance, '_output'];
outputfile = fopen(outputfilename, 'w');
fprintf(outputfile, [name_field, 'CIT/08',composition, '0G   %.3f  %.3f %.1f       1\n'], [Tmin, Tmax,Tmid])
for i = 8:12
    fprintf(outputfile, '%+.8E', x(i));
end
fprintf(outputfile, '    2\n');
for i = 13:14
    fprintf(outputfile, '%+.8E', x(i));
end
for i = 1:3
    fprintf(outputfile, '%+.8E', x(i));
end
fprintf(outputfile, '    3\n');
for i = 4:7
    fprintf(outputfile, '%+.8E', x(i));
end
fprintf(outputfile,'                   4\n');
fclose(outputfile);

outputfile = fopen(outputfilename, 'r');
outputstring = fscanf(outputfile,'%c' );
fclose(outputfile);

if ispc, outputstring= strrep(outputstring, 'E+0', 'E+');   end
if ispc, outputstring= strrep(outputstring, 'E-0', 'E-');   end

outputfile = fopen(outputfilename, 'w');
fprintf(outputfile,outputstring);
fclose(outputfile);

%Cantera Format output
outputfilename =[substance, '_Cantera_output'];
outputfile = fopen(outputfilename, 'w');
fprintf(outputfile, ['species( ', Cantera_output_text, '\n\t', 'thermo = (\n\t', 'NASA( [ ', '%.3f ', ', ', '%.3f', ' ], [ '], [Tmin, Tmid] )
for i = 1:6
    fprintf(outputfile, '%+.8E', x(i));
    fprintf(outputfile, ', ');
end
    fprintf(outputfile, '%+.8E', x(7));
        fprintf(outputfile, [ ' ] ),', '\n\t', 'NASA( [ ', '%.3f ', ', ', '%.3f ', ' ], [ ' ], [Tmid, Tmax]);
for i = 8:13
    fprintf(outputfile, '%+.8E', x(i));
    fprintf(outputfile, ', ');
end
    fprintf(outputfile, '%+.8E', x(14));
        fprintf(outputfile, [  ' ] )', '\n\t', '),', '\n\t', 'note = "CIT/08" ) ' ]);

outputfile = fopen(outputfilename, 'r');
outputstring = fscanf(outputfile,'%c' );
fclose(outputfile);

if ispc, outputstring= strrep(outputstring, 'E+0', 'E+');   end
if ispc, outputstring= strrep(outputstring, 'E-0', 'E-');   end

outputfile = fopen(outputfilename, 'w');
fprintf(outputfile,outputstring);
fclose(outputfile);