%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Constrained Least Squares Formulation %%%%
%%%%%%%% For Thermodynamic Properties %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% By Jack Ziegler Jan. 2008
% jackalak@caltech.edu
%
% This demo interpolates the enthalpy, pressure independent entropy,
% and specific heat at constant pressure.  It outputs the 
% coefficients for the piecewise functions for use with the
% CHEMKIN or Cantera software (using the ck2cti program
% to convert the CHEMKIN DATA to a Cantera Mechanism (.cti file). 
% See the Shock and Detonation toolbox documentation, and also,
% see the CHEMKIN and Cantera Manuals for more information
%
% This demo also creates output for a cantera .cti file, one may then
% directly cut and paste this data into an existing or new .cti file
%
%
% 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)
%
% Run this file and input as directed at the matlab prompt
% This Demo is set up for fiting 2-butenal, 
% with the twobutenal.m file as additional input.

clc;
close all;
clear all;

% This must match the input matrix (DATA) being used
Tmin = 200; % in Kelvin
Tmax = 5000;
Tmid = 2500; 

% Get the DATA matrix from this script file
% This data was simply copied and pasted from the output of the statistical
% mechanics based thermodynamic package LINGRAF in to the twobutenal.m file
% The data is defined at temperatures every 100 K from Tmin to Tmax
run twobutenal % here twobutenal is the input data script

% Note that the below are case sensitive!
substance = 'CH3CHCHCHO'; %substance name to apear in CHEMKIN II output
                                    %the CHEMKIN output file appears as substance_output
%  Modify the below string for your desired species name and its compostion
% "For the Cantera Mechanism File (.cti)"
Cantera_output_text = ['name = "CH3CHCHCHO",\n\t', 'atoms = " C:4 H:6 O:1 ",' ]; 

%  Standard state entropy and enthapy of formation at approximately 298 K                                   
so_over_R_298 =  40.06;  % so/R at 298K (From Wang 2000)
ho_over_RT_298 = -109.7/8.314472 *1000/298.15; % ho/RT at 298K (enthalpy of formation) 
% http://webbook.nist.gov/cgi/cbook.cgi?ID=C4170303&Units=SI&Mask=1#Thermo-Gas                                   

                                    
% type the composition here or use the below user prompt    
% Note that these formats follow the old CHEMKIN standard
% see http://www.galcit.caltech.edu/EDL/public/formats/chemkin.html

%composition = 'C   4H   6O   1    ';   % type in the composition, here there are 4 carbon
                                               % atoms and 7 hydrogen here 
                                               %(do not change the number of total characters
                                               % in this string (including spaces!!!!)
                                               
element_types = input('input number of element types ') ;
if element_types > 3
    disp('more than 3 elements is not supported');
end
element = input('first element (in CAPS) ', 's');
number = input('multiplicity ', 's') ;
composition  = [element , '   ', number];
if element_types > 1
    element = input('second element (in CAPS) ', 's');
    number = input('multiplicity ', 's') ;
    composition  = [composition,element , '   ', number];
end
if element_types > 2
    element = input('third element (in CAPS) ', 's');
    number = input('multiplicity ', 's') ;
    composition  = [composition,element , '   ', number];
end
comp_char_number = size(composition,2);
number_spaces = 19 - comp_char_number;
for i = 1:number_spaces
    composition =[composition, ' '];
end

% See interp_thermo.m
interp_thermo(substance, composition, Cantera_output_text, DATA, Tmin, Tmid, Tmax, so_over_R_298,ho_over_RT_298)