%   he0.m  Post-process Fluent calculations of the Barlow & Carter Flame He0
%   Radial profiles from two Fluent calculations are compared to the
%   experimental data.

FLA  = 'radial_base.dat';  % name of 1st Fluent output file (red)
KFLA = 1;  %  =1 for equilibrium calculation,  =2 for flamelet calculation
FLB  =  'radial_base.dat';  % name of 2nd Fluent output file  (blue)
KFLB = 1;  %  =1 for equilibrium calculation,  =2 for flamelet calculation


dmm = 3.75;        % jet diameter (mm)
d   = dmm .*1e-3;  % jet diamter (m)
Lod = 180.;        % L/d - L = visual length
L   = Lod .* d;
xoL = ['0   '; '1/16'; '1/8 '; '1/4 '; '3/8 '; '1/2 '; '5/8 '; '3/4 '; '1   '; '1.39'];
%                                  1 2 3 4 5  6   7 8
% columns in expt. velocity data: nr x y z u varu v	varv	w	varw	k

[vny,v]=get_vel(dmm);

% columns in the scalar data
%  1     2     3    4    5    6     7    8       9     10    11     12
% r(mm) Fblgr Frms T(K) Trms YO2  YO2rms YN2   YN2rms  YH2  YH2rms YH2O  

[sny,s]=get_scalar(dmm);

%  13      14   15    16     17       18       19
% YH2Orms YOH YOHrms  YNO   YNOrms   XNO(ppm) XNOrms

% columns of Fluent data for equilibrium calculation:
% 1 2 3 4 5 6  7  8  9  10 11  12 13  14
% n x y U T k eps H2 N2 O2 H2O OH f var(f)

% columns of Fluent data for flamelet calculation:
% 1 2 3 4 5 6  7  9  12  8 11  10 13  14
% n x y U T k eps H2 N2 O2 H2O OH f var(f)

convert = [1:14]';
convert(8)=9; convert(9)=12; convert(10)=8; convert(12)=10;  % equil index to flamelet index

% 
%  replaced in fl_unpack with:
% 6-sqrt(k)  7-chi  14-rms(f)
kx=2; ky=3; kU=4; kT=5; ktke=6; kH2O=11; kOH=12; kf=13; kfp=14;

faa   = load(FLA); %  first Fluent case
[ax any fa]=fl_unpack( faa, 2, 3, L, d );

fbb   = load(FLB); %  second Fluent case
[bx bny fb]=fl_unpack( fbb, 2, 3, L, d );

rmax = [5 10 15 20 25 30 40 50 50]';  % r/d ranges

figure(1)  %  mean velocity profiles
for k = 1:9
    kk=k+1;
    na =any(kk);
    nb =bny(kk);
    nv = vny(kk);
    subplot(3,3,k)
    plot( fa(1:na,ky,kk), fa(1:na,kU,kk), 'r-', ...
          fb(1:nb,ky,kk), fb(1:nb,kU,kk), 'b--', ...
          v(1:nv,1,kk), v(1:nv,5,kk), 'go' )
    xlim([0 rmax(k)]);
    legend(['x/L=' xoL(kk,:)]);
    if k>6; xlabel('r/d'); end
    if mod(k,3)==1; ylabel('U (m/s)'); end
end

figure(2)  %  sqrt(k) vs. u'
for k = 1:9
    kk=k+1;
    na = any(kk);
    nb = bny(kk);
    nv = vny(kk);
    subplot(3,3,k)
    plot( fa(1:na,ky,kk), fa(1:na,ktke,kk), 'r-', ...
          fb(1:nb,ky,kk), fb(1:nb,ktke,kk), 'b--', ...
          v(1:nv,1,kk), v(1:nv,6,kk), 'go' )
    xlim([0 rmax(k)]);
    legend(['x/L=' xoL(kk,:)]);
    if k>6; xlabel('r/d'); end
    if mod(k,3)==1; ylabel('rms u (m/s)'); end
end

nfig=7;
for ifig=3:nfig
    if ifig==3; kab=kf; ks=2; ylab='mean \xi'; end 
    if ifig==4; kab=kfp; ks=3; ylab='rms \xi'; end 
    if ifig==5; kab=kT; ks=4; ylab='mean T (K)'; end 
    if ifig==6; kab=kH2O; ks=12; ylab='Y_{H_2O}'; end 
    if ifig==7; kab=kOH; ks=14; ylab='Y_{OH}'; end 
    figure(ifig)  
    for k = 1:9
        kk=k+1;
        na = any(kk);
        nb = bny(kk);
        ns = sny(kk);
        ka = kab; if KFLA==2; ka=convert(kab); end
        kb = kab; if KFLB==2; kb=convert(kab); end
        subplot(3,3,k)
        plot( fa(1:na,ky,kk), fa(1:na,ka,kk), 'r-', ...
            fb(1:nb,ky,kk), fb(1:nb,kb,kk), 'b--', ...
            s(1:ns,1,kk), s(1:ns,ks,kk), 'go' )
        xlim([0 rmax(k)]);
        legend(['x/L=' xoL(kk,:)]);
        if k>6; xlabel('r/d'); end
        if mod(k,3)==1; ylabel(ylab); end
    end
end