function [Tad] = NASA_TN_D_5452(FA,Pcomb,varargin)
%NASA_TN_D_5452 estimates combustion product temperatures for Natural Gas
%               air.
%
%Reference: Poferl, Davis J., et al, "Thermodynamic and Transport
%Properties of Air and the Combustion Products of Natural Gas and of
%ASTM-A-1 Fuel with Air" (NASA TN D-5452); National Aeronautics and Space
%Administration, Washington, D.C., October 1969)
%
%This is just a crude regression based on experimentation using an old
%version of Matlab.  It shouldn't be used for any production purposes.
%
%Mandatory Inputs:
%
%FA     Ratio of Fuel to that of Air
%Pcom   The pressure of the reactants
%
%Optional Input
%
%conf   The desired confidence interval (default = 0.95).
%
%Output
%
%Tad    A vector: [-conf,"nominal",+conf];
%
%If no output is specified, a formatted string of the same data.

%Built with 2007b.


persistent x X Y B BINT R RINT STATS;

if isempty(X)
    %Build the independent columns of table III:
    x1 = [0;1;2;4;6]/100;%F/A
    x1 = [x1;x1];

    x2 = [ones(5,1)* 3.04e5;ones(5,1)* 10.13e5];%Pressure
    x = [x1 x2];

    %Build the regression model:
    X = theModel(x);
        
end;

if isempty(Y)
    %Build the dependent column (combustion temperature) of table III:
    Y = [298;
        745;
        1130;
        1770;
        2242;
        298;
        745;
        1130;
        1771;
        2262];
end;

%Check for stupid:

if ~((FA >= min(x(:,1)) && FA <= max(x(:,1))))
    error('Fuel-Air Ratio out of regression range.');
end;
if ~((Pcomb >= min(x(:,2)) && Pcomb <= max(x(:,2))))
    error('Reactant Pressure out of regression range.');
end;

%Check for non-default confidence interval:
if nargin == 2
    conf = 0.95;
else
    conf = varargin{1};
    %Check for stupid:
    if isnumeric(conf) && conf > 0 && conf < 1 
    else
        error('Confidence limit must be > 0 and < 1.');
    end;
end;

%Only recalc if required to do so:
if any(isempty(B) || isempty(BINT) || isempty(R) || isempty(RINT) || isempty(STATS) || conf~=0.95)
    %The terms of B are the coefficients for each term of the fitted model,
    %in the same order as given in the X matrix.  For the rest of the
    %variables, see the Matlab help for the regress command.
    [B,BINT,R,RINT,STATS] = regress(Y,X,1-abs(conf)/100);
    %Although the R-square, F, and p values indicate a good fit, the RINT
    %output suggests that it has some odd things going on - small, but odd.
    %It isn't all that great of a model and, for production usage, I'd
    %probably want to track things back to the original model as
    %implemented for the reference.  But this will do for exposition.
    
end;

Y1 = theModel([FA,Pcomb]);

Tad = [sum((BINT(:,1)) .* Y1'),sum(B .* Y1'),sum((BINT(:,2)) .* Y1')];

Tad = round(1e4*Tad)/1e4;%A valid argument could be made for only three significant digits...

%Form the output based on what was requested by the calling workspace:
if nargout == 0
    Tad = sprintf('[%4.0f,%4.0f,%4.0f] Kelvin\n',Tad);
else
    display(sprintf('Temperature results are in Kelvin @ [%0.2f,0.50,%0.2f] Confidence',1-conf,conf));
end;
end

function rslt = theModel(inp)
%The model shape resulted from simple fiddle-farting around for a few
%minutes given inspection of the curve shapes. A MUCH more thorough fit
%would be possible!
rslt = [inp(:,1).^3,inp(:,1).^2,inp(:,1),(inp(:,1).^6).*inp(:,2).^1.5,ones(size(inp,1),1)];

end