% Shock and Detonation Toolbox Utility Program
%
% thermo_refit.m - Plots thermodynamic data from .cti input file and 
% evaluates standard states for specified list of species in file
% "BadSpecies.txt" generated by running thermo_check.py
% Options:
% 1. Refits CP data and recomputes NASA-7 coefficients to eliminate jumps at midpoints
% 2. Output data in Cantera format
% 3. Output data in NASA 7 format or NASA 9 format
% 4. Reads a list of input species from file BadSpecies.txt created by
%    checkThermo.py
% 5. Runs in noninteractive mode
%
% ################################################################################
% Theory, numerical methods and applications are described in the following report:
% 
%     Numerical Solution Methods for Shock and Detonation Jump Conditions, S.
%     Browne, J. Ziegler, and J. E. Shepherd, GALCIT Report FM2006.006 - R3,
%     California Institute of Technology Revised September, 2018.
% 
% Please cite this report and the website if you use these routines. 
% 
% Please refer to LICENCE.txt or the above report for copyright and disclaimers.
% 
% http://shepherd.caltech.edu/EDL/PublicResources/sdt/
% 
% ################################################################################ 
% Updated September 2018
% Tested with: 
%     MATLAB 2017b and 2018a, Cantera 2.3 and 2.4
% Under these operating systems:
%     Windows 8.1, Windows 10, Linux (Debian 9)
%
clc;clear;close all;
silent = true;  %  runs in noninteractive mode from the thermo_check.py output file, no plot output
file = true;  % read list of species from the file created by thermo_check.py
refit = true; % refit data using least squares to fit Cp
cti = false;   % write cantera .cti format output file (only if refit)
nasa7 = true;  % write nasa 7 (Chemkin II) .dat format output file  (only if refit)
nasa9 = false;  % write nasa 9  .dat format output file  (only if refit)
ck_files = 0;
n9_files = 0;
cti_files = 0;
Ra = gasconstant;
T_zero = 298.15;
P_zero = oneatm;
% Cantera mechanism with thermo data of interest
%mechanism='keromnes2012';
%mechanism='Hong2011';
%mechanism='OH-Ruscic-Burcat';
%mechanism='OH-IUPAC';
%mechanism='OH-partition-SDT';
mechanism ='GRI30-x';
%mechanism = 'JetSurf2';
%mechanism = 'chem';
%mechanism = 'Mevel2018';
%mechanism = 'Mevel2015';
%mechanism = 'hexaneFull'
%mechanism = 'hexanePartial';
%mechanism = 'hexaneReduced';
%mechanism = 'sandiego20161214_H2only';
dir_NAME = mechanism;
if (file)
    status = mkdir(dir_NAME);
end
if (silent)
    log_file = [dir_NAME,'/refit_log.txt'];
    diary(log_file)
end
mech = [mechanism '.cti'];
gas = Solution(mech);
% get species name and temperature limits from text file created by user
if (file)
    fileID = fopen('BadSpecies.txt','r');
    format = '%s';
    S = textscan(fileID, format);
    fclose(fileID);
    Lines = char([S{1}]);
    [nl, nc] = size(Lines);
    nsp = nl/4;
    i = 0;
    sps = struct('name',[],'Tmin',[],'Tmid',[],'Tmax',[]);
    for k=1:4:nl
        i = i + 1;
        sps(i).name = strtrim(Lines(k,:));
        sps(i).Tmin = Lines(k+1,:);
        sps(i).Tmid = Lines(k+2,:);
        sps(i).Tmax = Lines(k+3,:);
    end
    display(['Read ',num2str(nsp),' species from BadSpecies.txt']);
else
    % or only do a single species specified below
    nsp = 1;
    display(['Mechanism file: ',mechanism,'\n species: ']);
    All_species = speciesNames(gas);
    nTotal = nSpecies(gas);
    for i = 1:nTotal
        display(char(All_species(i)));
    end
end

%% loop through all the species in the file (nsp> 1) or input species (nsp = 1)
for j = 1:nsp
    if (file)
        Tmins = str2double(sps(j).Tmin);
        Tmids = str2double(sps(j).Tmid);
        Tmaxs = str2double(sps(j).Tmax);
        species = sps(j).name;
    else
        species  = input('Species ','s');
        Tmins  = input('Minimum temperature ');
        Tmids  = input('Midpoint temperature ');
        Tmaxs  = input('Msximum  temperature ');
    end
    Tmin = Tmins-rem(Tmins,10);  % round down to decade to eliminate issues in refit
    Tmid = Tmids;
    Tmax = Tmaxs;
    T = Tmin:10:Tmax;
    imax = length(T);
    CP_R = zeros(1,imax); H_RT = zeros(1,imax);S_R=zeros(1,imax);
    i_species = speciesIndex(gas,species);
    X = [species,':1.0'];
    ne = nElements(gas);
    Natoms = 0;
    for k =1:ne
        el(k) = elementName(gas,k);    %name of elements in gas object
        at(k) = nAtoms(gas,species,elementIndex(gas,el(k))); %number of atoms of each element in the species
        Natoms = Natoms + at(k);
    end
    els = char(el);  %convert to character string from cell array
    for i=1:imax
        set(gas, 'T', T(i), 'P', oneatm, 'X', X);
        cp_all = cp_R(gas);
        h_rt_all = enthalpies_RT(gas);
        s_r_all = entropies_R(gas);
        CP_R(i) = cp_all(i_species);
        H_RT(i) = h_rt_all(i_species);
        S_R(i) = s_r_all(i_species);
    end
    set(gas, 'T', T_zero, 'P', oneatm, 'X', X);
    w = molecularWeights(gas);
    m = w(i_species);
    cp_all = cp_R(gas);
    CP_R_zero = cp_all(i_species);
    h_rt_all = enthalpies_RT(gas);
    H_RT_zero = h_rt_all(i_species);
    DeltaH_f0 = H_RT_zero*Ra*T_zero;
    s_r_all = entropies_R(gas);
    S_R_zero = s_r_all(i_species);
    if (file)
        display (['**** species number ',num2str(j),' in input file ']);
    end
    display(['Mechanism: ',mech]);
    display(['Species: ',species]);
    display(['Total number of atoms ', num2str(Natoms)]);
    display(['Standard state values:']);
    display(['Cp/R                ', num2str(CP_R_zero)]);
    display(['H/RT                ', num2str(H_RT_zero)]);
    display(['S/R                 ', num2str(S_R_zero)]);
    display(['Delta H_f0 (kJ/mol) ', num2str(DeltaH_f0/1E6)]);
    display(['S_0         (J/mol) ', num2str(S_R_zero*Ra/1E3)]);
    display(['Minimum fit temperature (K) = ',num2str(Tmins)]);
    set(gas, 'T', Tmins, 'P', oneatm, 'X', X);
    cp_all = cp_R(gas);
    CP_R_min = cp_all(i_species);
    display(['Cp/R  (Tmin)        ', num2str(CP_R_min)]);
    display(['Maximum fit temperature (K) = ',num2str(Tmaxs)]);
    set(gas, 'T', Tmaxs, 'P', oneatm, 'X', X);
    cp_all = cp_R(gas);
    CP_R_max = cp_all(i_species);
    display(['Cp/R  (Tmax)        ', num2str(CP_R_max)]);
    display(['Midpoint fit temperature (K) = ',num2str(Tmids)]);
    if (Natoms == 1)
        cplimit = 5./2.;     % atom
    elseif (Natoms == 2)
        cplimit = 9./2.;     % diatomic
    elseif (Natoms > 2)
        if (~silent)
            geo = input('Linear (L) or nonlinear (N) molecule ?','s');
        else
            geo = 'N';
        end
        if (strcmp(geo,'n') || strcmp(geo,'N'))
            cplimit = (6*Natoms - 4)/2 ;   % nonlinear multiatomic
            display(['geometry nonlinear']);
        elseif (strcmp(geo,'l') || strcmp(geo,'L'))
            cplimit = (6*Natoms - 3)/2 ;   % linear multiatomic
            display(['geometry linear']);
        else
            cplimit = (6*Natoms - 4)/2 ;   % default nonlinear
            display(['geometry nonlinear']);
        end
    end 
    display(['Cp/R  (T_infty)        ', num2str(cplimit)]);
    %
    % Wang extrapolation method for T > Tmaxs
    %     for i = 1:imax
    %         if (T(i) >= Tmaxs)
    %             CPI = 1./cplimit  -(1./cplimit - 1./CP_R_max)*Tmaxs/T(i);
    %             CP_R(i) = 1/CPI;
    %         end
    %     end
    %% plot thermo data vs temperature
    if (~silent)
        figure('Name','Cantera thermo' );
        maxy = max(H_RT(:));
        miny = min(H_RT(:));
        maxx = max(T(:));
        minx = min(T(:));
        fontsize=12;
        axis([minx maxx miny maxy]);
        hold on;
        set(gca,'FontSize',fontsize,'LineWidth',2);
        plot(T(:),H_RT(:));
        xlabel('T (K)','FontSize',fontsize);
        ylabel('h/RT','FontSize',fontsize);
        title(['Cantera data ', mechanism, ' species ', species],'FontSize',fontsize);
        hold off
        
        % plot thermo data vs temperature
        figure('Name','Cantera thermo' );
        maxy = max(S_R(:));
        miny = min(S_R(:));
        maxx = max(T(:));
        minx = min(T(:));
        fontsize=12;
        axis([minx maxx miny maxy]);
        hold on;
        set(gca,'FontSize',fontsize,'LineWidth',2);
        plot(T(:),S_R(:));
        xlabel('T (K)','FontSize',fontsize);
        ylabel('S/R','FontSize',fontsize);
        title(['Cantera data ', mechanism, ' species ', species],'FontSize',fontsize);
        hold off
        
        figure('Name','Cantera thermo' );
        maxy = max(CP_R(:));
        miny = min(CP_R(:));
        maxx = max(T(:));
        minx = min(T(:));
        fontsize=12;
        axis([minx maxx miny maxy]);
        hold on
        set(gca,'FontSize',fontsize,'LineWidth',2);
        plot(T(:),CP_R(:));
        xlabel('T (K)','FontSize',fontsize);
        ylabel('C_p/R','FontSize',fontsize);
        title(['Cantera data ', mechanism, ' species ', species],'FontSize',fontsize);
        hold off
        display(['Paused, CR to continue']);
        pause;
        close all;
    end
    %% refit the data
    if (refit)
        if (~silent)
            str = input('Midpoint temperature for new fit ? ','s'); % for best results this should not be original value if Cp jumps at Tmid
        else
            str = '';
        end
        if isempty(str)
            Tmid = Tmids  + 200.;
        else
            Tmid = str2num(str);
            if ((Tmid == Tmids) | (Tmid < Tmin) | (Tmid > Tmax))
                Tmid = Tmids  + 200.;
            end
        end
        % find specific heat derivative at midpoint
        index_mid = max(find(T(:)<=Tmid));   %find closest temperature to selected value
        Tmid = T(index_mid:index_mid);
        display(['Adjusted temperature for midpoint of fit = ',num2str(Tmid)]);
        display(['Adjusted temperature for minimum of fit = ',num2str(Tmin)]);
        display(['Refitting...']);
        z = [T; CP_R; H_RT; S_R];
        DCPRDT = (CP_R(index_mid+1)-CP_R(index_mid-1))/(T(index_mid+1) - T(index_mid-1));
        [x, cp_fit, index_mid, a_0_6, a_1_6, s_fit,  a_0_7, a_1_7, h_fit] = poly_cp(z,Tmin,Tmid,Tmax,CP_R_zero,DCPRDT,S_R_zero,H_RT_zero);
        %% Check for jumps and maximum differences with original data
        cp_jump = cp_fit(index_mid+1)-cp_fit(index_mid);
        h_jump = h_fit(index_mid+1)-h_fit(index_mid);
        s_jump = h_fit(index_mid+1)-h_fit(index_mid);
        display(['Result of optimization after adjustment']);
        display(['Jump at Tmid in Cp/R ', num2str(cp_jump)]);
        display(['Jump at Tmid in H/RT ', num2str(h_jump)]);
        display(['Jump at Tmid  in S/R ', num2str(s_jump)]);
        %Compute maximum errors (remove duplicate value at Tmid)
        cp_fit(index_mid:index_mid) = [];
        h_fit(index_mid:index_mid) = [];
        s_fit(index_mid:index_mid) = [];
        cp_diff = cp_fit'-CP_R;
        h_diff = h_fit'- H_RT;
        s_diff = s_fit'-S_R;
        cp_max_diff = max(abs(cp_diff));
        h_max_diff = max(abs(h_diff));
        s_max_diff = max(abs(s_diff));
        display(['Maximum difference in Cp/R ', num2str(cp_max_diff)]);
        display(['Maximum difference in H/RT ', num2str(h_max_diff)]);
        display(['Maximum difference in S/R ', num2str(s_max_diff)]);
        if (~silent)
            % plot differences
            figure('Name','Original vs refit data');
            fontsize=12;
            hold on;
            set(gca,'FontSize',fontsize,'LineWidth',2);
            plot(T(:),cp_diff(:));
            plot(T(:), h_diff(:));
            plot(T(:), s_diff(:));
            xlabel('T (K)','FontSize',fontsize);
            ylabel('Differences','FontSize',fontsize);
            legend('Cp/R', 'H/RT', 'S/R','Location','southeast');
            title(['Cantera data ', mechanism, ' species ', species],'FontSize',fontsize);
            hold off
            
            % plot thermo data vs temperature
            figure('Name','New fit of original data' );
            maxy = max(cp_fit(:));
            miny = min(cp_fit(:));
            maxx = max(T(:));
            minx = min(T(:));
            fontsize=12;
            axis([minx maxx miny maxy]);
            hold on
            set(gca,'FontSize',fontsize,'LineWidth',2);
            plot(T(:),cp_fit(:));
            plot(T(1:20:imax),CP_R(1:20:imax),'ok');
            xlabel('T (K)','FontSize',fontsize);
            ylabel('C_p/R','FontSize',fontsize);
            title(['Cantera data ', mechanism, ' species ', species],'FontSize',fontsize);
            legend('New fit','Old Fit','Location','southeast');
            hold off
            
            % plot thermo data vs temperature
            figure('Name','New fit of original data' );
            maxy = max(h_fit(:));
            miny = min(h_fit(:));
            maxx = max(T(:));
            minx = min(T(:));
            fontsize=12;
            axis([minx maxx miny maxy]);
            hold on;
            set(gca,'FontSize',fontsize,'LineWidth',2);
            plot(T(:),h_fit(:));
            plot(T(1:20:imax),H_RT(1:20:imax),'ok');
            xlabel('T (K)','FontSize',fontsize);
            ylabel('h/RT','FontSize',fontsize);
            title(['Cantera data ', mechanism, ' species ', species],'FontSize',fontsize);
            legend('New fit','Old Fit','Location','southeast');
            hold off
            
            % plot thermo data vs temperature
            figure('Name','New fit of original data' );
            maxy = max(s_fit(:));
            miny = min(s_fit(:));
            maxx = max(T(:));
            minx = min(T(:));
            fontsize=12;
            axis([minx maxx miny maxy]);
            hold on;
            set(gca,'FontSize',fontsize,'LineWidth',2);
            plot(T(:),s_fit(:));
            plot(T(1:20:imax),S_R(1:20:imax),'ok');
            xlabel('T (K)','FontSize',fontsize);
            ylabel('S/R','FontSize',fontsize);
            title(['Cantera data ', mechanism, ' species ', species],'FontSize',fontsize);
            legend('New fit','Old Fit','Location','southeast');
            hold off
            display(['Paused, CR to continue']);
            pause;
            close all;
        end
        %% create coefficient in correct order for NASA-7 format
        %The NASA-7 polynomials are of the form:
        %Cp/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4
        %H/RT = a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T
        %S/R  = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7
        % a = [x(1:5)',a_0_6,a_0_7,x(6:10)',a_1_6,a_1_7];
        
        if (~silent)
            %% Adjust heat of formation?
            display(['Delta H_f0 (kJ/mol) ', num2str(DeltaH_f0/1E6)]);
            tr = input('New value of heat of formation (kJ/mol) ? ','s');
            if isempty(str)
                delta_a6 = 0;
            else
                h_new = str2num(str);
                delta_a6 = (h_new*1e6-DeltaH_f0)/Ra;
            end
            a_0_6 = a_0_6+delta_a6;
            a_1_6 = a_1_6+delta_a6;
            %% Adjust standard state entropy?
            display(['S_0         (J/mol) ', num2str(S_R_zero*Ra/1E3)]);
            str = input('New value of standard state entropy (J/mol) ? ','s');
            if isempty(str)
                delta_a7 = 0;
            else
                s_new = str2num(str);
                delta_a7 = (s_new*1e3-S_R_zero*Ra)/Ra;
            end
            a_0_7 = a_0_7 + delta_a7;
            a_1_7 = a_1_7 + delta_a7;
        end
        a = [x(1:5)',a_0_6,a_0_7,x(6:10)',a_1_6,a_1_7];
        %% Cantera .cti format output file
        if (cti)
            display(['Writing Cantera format .cti file for species ',species]);
            format = ' %c:%d ';
            atoms = '';
            for k = 1:ne
                if (at(k) > 0)
                    str = sprintf(format,els(k),at(k));
                    atoms = strcat(atoms,str);
                end
            end
            Cantera_output_text = ['name = "',species,'",\n',...
                '        atoms = "',atoms, '",' ];
            species_clean = strrep(species,'*','star');
            outputfilename =[dir_NAME,'/',species_clean,'-refit','.cti'];
            outputfile = fopen(outputfilename,'w');
            fprintf(outputfile, ['species(', Cantera_output_text,'\n', ...
                '        thermo=(NASA([%.3f, %.3f],\n',...
                '                     ['], Tmin, Tmid);
            for i = 1:6
                fprintf(outputfile, '%+.8E, ', a(i));
            end
            fprintf(outputfile, '%+.8E]),\n', a(7));
            fprintf(outputfile, ['                NASA([%.3f, %.3f],\n',...
                '                     ['], Tmid, Tmax);
            for i = 8:13
                fprintf(outputfile, '%+.8E, ',a(i));
            end
            fprintf(outputfile, '%+.8E])', a(14));
            fprintf(outputfile, [ '),\n',...
                '        note = ''Refit from ',mechanism,' by JES ', datestr(now),''')']);
            outputfile = fopen(outputfilename, 'r');
            outputstring = fscanf(outputfile,'%c' );
            fclose(outputfile);
            cti_files = cti_files + 1;
        end
        %% NASA (Chemkin II) format output file
        if (nasa7)
            display(['Writing NASA (Chemkin II) format .dat file for species ',species]);
            format = '%2c%3i';
            atoms = '';
            for k = 1:ne
                if (at(k) > 0)
                    str = sprintf(format,els(k),at(k));
                    atoms = strcat(atoms,str);
                end
            end
            % pad species names and atom fields to fixed length
            name_field = species;
            name_char_number = size(species,2);
            number_spaces = 18 - name_char_number;
            for i = 1:number_spaces
                name_field =[name_field,' '];
            end
            c_size = size(atoms,2);
            n_spaces = 20 -c_size;
            for i = 1:n_spaces
                atoms =[atoms,' '];
            end
            %NASA format is particular about columns alignment
            species_clean = strrep(species,'*','star');
            outputfilename =[dir_NAME,'/',species_clean,'.dat'];
            outputfile = fopen(outputfilename, 'w');
            %            fprintf(outputfile,'THERMO\n%10.3f%10.3f%10.2f\n',Tmin, Tmax,Tmid);
            space = '      ';
            fprintf(outputfile, '%s%s%s%s%10.2f%10.2f%8.2f%s1\n',name_field,' SDT18',atoms,'G',Tmin, Tmax,Tmid,space);
            for i = 8:12
                fprintf(outputfile, '%+.8E', a(i));
            end
            fprintf(outputfile, '    2\n');
            for i = 13:14
                fprintf(outputfile, '%+.8E', a(i));
            end
            for i = 1:3
                fprintf(outputfile, '%+.8E', a(i));
            end
            fprintf(outputfile, '    3\n');
            for i = 4:7
                fprintf(outputfile, '%+.8E', a(i));
            end
            fprintf(outputfile,'                   4\n');
            %            fprintf(outputfile,'END\n');
            fclose(outputfile);
            ck_files = ck_files + 1;
        end
        %% NASA 9 term (1996) format output file
        %
        %Record     Constants     Format     Column
        %1  Species name or formula A24 1 to 24
        %   Comments (data source) A56 25-80
        %2  Number of T intervals I2 2
        %   Optional identification code A6 4-9
        %   Chemical formulas, symbols, and numbers 5(A2,F6.2) 11-50
        %   Zero for gas and nonzero for condensed phases I1 52
        %   Molecular weight F13.5 53-65
        %   Heat of formation at 298.15 K, J/mol F13.5 66-80
        %3  Temperature range 2F10.3 2-21
        %   Number of coefficients for Cp0/R I1 23
        %   T exponents in empirical equation for Cp0/R 8F5.1 24-63
        %   {H0(298.15)-H0(0)}, J/mol F15.3 66-80
        %4  First five coefficients for Cp0/R 5D16.8 1-80
        %5  Last two coefficients for Cp0/R 3D16.8 1-48
        %   Integration constants b1 and b2 2D16.8 49-80
        %... Repeat 3, 4, and 5 for each interval
        %
        if (nasa9)
            display(['Writing NASA-9 format .dat file for species ',species]);
            format = '%2c%6.2f';
            atoms = '';
            for k = 1:5
                if (k <= ne)
                    if (at(k) > 0)
                        str = sprintf(format,strjust(els(k),'left'),at(k));
                        atoms = strcat(atoms,str);
                    else
                        str = sprintf(format,[' '],0.0);
                        atoms = strcat(atoms,str);
                    end
                else
                    str = sprintf(format,[' '],0.0);
                    atoms = strcat(atoms,str);
                end
            end
            % pad species names and atom fields to fixed length
            name_field = species;
            name_char_number = size(species,2);
            number_spaces = 24 - name_char_number;
            for i = 1:number_spaces
                name_field =[name_field,' '];
            end
            c_size = size(atoms,2);
            n_spaces = 40 -c_size;
            for i = 1:n_spaces
                atoms =[atoms,' '];
            end
            %NASA format is particular about columns alignment
            species_clean = strrep(species,'*','star');
            outputfilename =[dir_NAME,'/',species_clean,'-partition-NASA9','.dat'];
            outputfile = fopen(outputfilename, 'w');
            %            fprintf(outputfile,'thermo nasa9\n);
            space = '      ';
            % record 1
            fprintf(outputfile, '%s%s%s%s\n',name_field,'Refit from ',mechanism,' by JES ', datestr(now));
            % record 2
            fprintf(outputfile,' %1d %s %s %1d%13.5f%13.5f\n',2,space,atoms,0,m,DeltaH_f0);
            % record 3
            dummy = 0;  % need to set this to  H(298.15) - H(0)
            fprintf(outputfile,' %10.3f%10.3f %1d%5.1f%5.1f%5.1f%5.1f%5.1f%5.1f%5.1f%5.1f  %15.3f\n', Tmin, Tmid,7,-2.,-1.,0.,1.,2.,3.,4.,0.,dummy);
            %record 4
            fprintf(outputfile, '%16.8E%16.8E%16.8E%16.8E%16.8E\n', 0,0,a(1),a(2),a(3));
            %Record 5
            fprintf(outputfile, '%16.8E%16.8E%16.8E%16.8E%16.8E\n', a(4),a(5),0.0,a(6),a(7));
            fprintf(outputfile,' %10.3f%10.3f %1d%5.1f%5.1f%5.1f%5.1f%5.1f%5.1f%5.1f%5.1f  %15.3f\n', Tmid, Tmax,7,-2.,-1.,0.,1.,2.,3.,4.,0.,dummy);
            %record 4
            fprintf(outputfile, '%16.8E%16.8E%16.8E%16.8E%16.8E\n', 0,0,a(8),a(9),a(10));
            %Record 5
            fprintf(outputfile, '%16.8E%16.8E%16.8E%16.8E%16.8E\n', a(11),a(12),0.0,a(13),a(14));
            fclose(outputfile);
            n9_files =n9_files + 1;
        end
    end
end
display(['Wrote ',num2str(ck_files), ' NASA 7 (chemkin) data files']);
display(['Wrote ',num2str(n9_files), ' NASA 9 data files']);
display(['Wrote ',num2str(cti_files),' cti (cantera) data files']);

if (silent)
    diary off
end
    
function [x, cp_fit, index_mid, a_0_6, a_1_6, s_fit, a_0_7, a_1_7, h_fit] = poly_cp(z, Tmin, Tmid, Tmax,CP_R_zero,DCPRDT,S_R_zero,H_RT_zero)
% poly_cp.m - A function that fits a two segment polynomial to specific
% heat, use NASA 7 term format from SP-272. Gordon and McBride 1971.
% Coefficients a6 and a7 for enthalpy and entropy are also computed.
%
%Eq. 90: cp0/R=a1+a2T+a3*T^2+a4*T^3+a5*T^4 
%Eq. 91: H0_T/(RT)=a1+a2/2*T+a3/3*T^2+a4/4*T^3+a5/5*T^4+a6/T 
%Eq. 92: S0_T/R=a1*ln(T)+a2*T+a3/2*T^2+a4/3*T^3+a5/4*T^4+a7
%
% This program fits the specific heat only then uses standard state data
% and the analytical integration to obtain the additional coefficients for
% entropy and enthalpy.  The specific heat fit constrains the derivative to
% vanish at the lowest temperature, the value and derivative are continuous
% at the knot at T_mid and the value of specific heat at 298.15 K is
% constrained to be the specified value of the standard state. The value of
% the entropy and entropy are set to standard state values at 298.15 K. The
% entropy and enthalpy values at the midpoints are continuous.  This
% procedure produces high quality fits as long as the input entropy and enthalpy
% are thermodynamically consistent with the specific heat data.  This will
% always be the case if the thermodynamic data have been derived from
% partition functions.   
%
% Matlab R2018a
% JES 7/14/2018
%
T_zero = 298.15;
T = z(1,:);
cp_over_R = z(2,:);
h_over_RT = z(3,:);
s_over_R = z(4,:);
index_max = length(T);
Tmin = min(T);
Tmax = max(T);
index_mid = max(find(T(:)<=Tmid));   %find closest temperature in vector
Tmid = T(index_mid:index_mid);
%
% set up coefficient matrix
C = zeros(index_max+1,10);
for i = 1:5
    for k =1:index_mid
        C(k,i) = T(k)^(i-1);
    end
end
for i = 1:5
    for k = index_mid:index_max
        C(k+1,i+5) = T(k)^(i-1);
    end
end
% initial guess for coefficients
for i = 1:10
    x0(i)  = 1.;
end
% right hand side
for k = 1:index_mid
    d(k) = cp_over_R(k);
end
for k = index_mid:index_max
    d(k+1) = cp_over_R(k);
end
% constraint values 
for i=1:5
    Aeq(1,i) = T_zero^(i-1);   %fix standard state value for poly 1
    Aeq(1,i+5) = 0.0;
    %
    Aeq(2,i) = T(index_mid)^(i-1);   %fix midpoint value  for poly 1
    Aeq(2,i+5) = 0.0;
    %
    Aeq(3,i) = 0.;                    %fix midpoint value  for poly 2
    Aeq(3,i+5) = T(index_mid)^(i-1);
    %
    Aeq(4,i) = (i-1)*T(index_mid)^(i-2);  %fix derivative for poly 1
    Aeq(4,i+5) = 0.0;
    %
    Aeq(5,i) = 0.;                       %fix derivative for poly 2                 
    Aeq(5,i+5) = (i-1)*T(index_mid)^(i-2);
    %
%    Aeq(6,i) = (i-1)*T(1)^(i-2);    % sero derivative at lowest temperature
%    Aeq(6,i+5) = 0.0;
end
beq(1) = CP_R_zero;  % match standard state, polynomial 1
beq(2) = cp_over_R(index_mid);  % continuity at knot, polynomial 1
beq(3) = cp_over_R(index_mid);  % continuity at knot, polynomial 2
beq(4) = DCPRDT;   % match derivative at knot, polynomial 1
beq(5) = DCPRDT;   % match derivative at knot, polynomial 2
%beq(6) = 0.;        % zero derivative at lowest temperature in poly 1
% least squares fit, diagnostics can be turned on if necessary
%options = optimoptions('lsqlin','Algorithm','active-set','Diagnostics','on');
options = optimoptions('lsqlin','Algorithm','interior-point','Diagnostics','off');
x = lsqlin(C,d,[],[],Aeq,beq,[],[],x0,options);
% evaluate polynomial for specific heat
cp_fit = C*x;
% value at standard state
% check fit at standard state and eliminate any jump at midpoint
display(['Result of initial optimization']);
sum = 0.;
for i = 1:5
    sum = sum + x(i)*T_zero^(i-1);
end
delta_cp = CP_R_zero-sum;
display(['Difference in CP at T_zero = ',num2str(delta_cp)]);
% check for jump at T mid
cp_jump = cp_fit(index_mid+1)-cp_fit(index_mid);
display(['Jump in CP at midpoint = ',num2str(cp_jump)]);
% remove jump
x(6) = x(6) - cp_jump;
% update initial guess and re-run optimization
x1 = x;
x = lsqlin(C,d,[],[],Aeq,beq,[],[],x1,options);
x2 = x;
x = lsqlin(C,d,[],[],Aeq,beq,[],[],x2,options);
display(['Result of re-running optimization']);
%evaluate fit again
cp_fit = C*x;
sum = 0.;
for i = 1:5
    sum = sum + x(i)*T_zero^(i-1);
end
delta_cp = CP_R_zero-sum;
display(['Difference at T_zero in CP = ',num2str(delta_cp)]);
% check for jump at T mid
cp_jump = cp_fit(index_mid+1)-cp_fit(index_mid);
display(['Jump in CP at midpoint = ',num2str(cp_jump)]);
%
% evaluate the enthalpy fit
%
C2 = zeros(index_max+1,10);
for i = 1:5
    for k =1:index_mid
        C2(k,i) = T(k)^(i-1)/i;
    end
end
for i = 1:5
    for k = index_mid:index_max
        C2(k+1,i+5) = T(k)^(i-1)/i;
    end
end
% fix poly 1 constant by matching enthalpy at 298.15 
sum = 0.;
for i = 1:5
   sum = sum + x(i)*T_zero^(i-1)/i;
end
a_0_6 = T_zero*(H_RT_zero - sum);
% fix poly 2 constant by matching entropy at knot 
sum1 = 0.;
for i = 1:5
   sum1 = sum1 + x(i)*T(index_mid)^(i-1)/i;
end
sum1 = sum1 + a_0_6/T(index_mid);
% match polynomials at knot
sum2 = 0.;
for i = 1:5
   sum2 = sum2 + x(5+i)*T(index_mid)^(i-1)/i;
end
a_1_6 = T(index_mid)*(sum1 - sum2);
% evaluate corrected fit
h_fit=C2*x; 
for k = 1:index_mid
    h_fit(k) = h_fit(k) + a_0_6/T(k);
end
for k = index_mid:index_max
    h_fit(k+1) = h_fit(k+1) + a_1_6/T(k);
end
%
% evaluate entropy fit
% 
C1 = zeros(index_max+1,10);
for i = 1:5
    for k =1:index_mid
        if (i==1)
            C1(k,i) = log(T(k));
        else
            C1(k,i) = T(k)^(i-1)/(i-1);
        end
    end
end
for i = 1:5
    for k = index_mid:index_max
        if (i==1)
            C1(k+1,i+5) = log(T(k));
        else
            C1(k+1,i+5) = T(k)^(i-1)/(i-1);
        end
    end
end
% fix poly 1 constant by matching entropy at 300 K, gives lower error than
% matching interpolated value at T_zero
s_fit=C1*x; 
index_300 = find(T(:)==300);
S_R_300 = s_over_R(index_300);
a_0_7 = S_R_300 - s_fit(index_300) ;
% fix poly 2 constant by matching entropy at knot 
a_1_7 =(s_fit(index_mid) + a_0_7) - s_fit(index_mid+1);
% evaluate corrected fit
for k = 1:index_mid
    s_fit(k) = s_fit(k) + a_0_7;
end
for k = index_mid:index_max
    s_fit(k+1) = s_fit(k+1) + a_1_7;
end
end
