Figure 9.36 (page 558):

Fit of LC measurement versus time for reduced model; early time measurements have been added.

Code for Figure 9.36

Text of the GNU GPL.

main.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

global vrdiclo k1k2ratio vro teaf flowrate tflow vro timeout idx ...
      lcmeas fcn x xdot sx sxdot x0 sx0 tcrit

%%
%% global variables                     
%%
%% vrdiclo                      initial moles of dichloro
%% k1k2ratio                    k1/k2
%% vro	                        initial reactor volume  cu dm
%% teef	                        TEA feed concentration kgmol/cu dm
%% teaden                       TEA density kg/cu dm
%% tflow, flowrate              measured flowrate at time tflow
%%                              read from "flow.dat" containing time
%%                              (min) and flowrate (kg/min) converted to
%%                              (cu dm/min)
%% tlc, lcmeas                  lc measurement at time tlc (min)
%%
%% optimal values
vrdiclo     =  2.3497;
k1k2ratio   =  2;
%%
%%
%%
vro       = 2370;
teaf      = 0.00721;
teaden    = 0.728;
flow = load ('flow.dat')
lc = load ('lc.dat')
tflow  = [0. ; flow(:,1)];
flowrate = [0. ; flow(:,2)./teaden];
tlc    = lc(:,1);
lcmeas = lc(:,2);
%%
%% Simulated the lc measurement from the beginning of the experiment.
%% These data were saved by hand in lcsim.dat with 
%% k1k2ratio=2, vrdiclo=2.3497;
%% 2nd column is without noise, 3rd column is with noise.
%%
%% prepend simulated data
%%
lcsim = load ('lcsim.dat')
tlc = [lcsim(:,1); tlc];
lcmeas = [lcsim(:,3); lcmeas];


tcrit = flow(:,1);

ntimes    = 80;
tlin = linspace (0, tflow(end), ntimes)';
timeout = sort (unique ([tlin; tflow; tlc]));

%% I can't think of a way to do this without a loop.
idx = zeros (size (tlc));
for i = 1:length (tlc)
  %% This assignment will fail if find returns more than one element,
  %% but that's OK, because each element of tlc should only appear once
  %% in timeout.
  idx(i) = find (timeout == tlc(i));
end

%% sanity check
if (any (timeout (idx) ~= tlc))
  error ('extracting index vector for tlc');
end



% function retval = gradient(theta)
%   global vrdiclo k1 k2 vro teaf flowrate tflow xsmall idx fcn sx0 ...
%       timeout tcrit x xdot sx sxdot sens lcmeas resid grad Hgn ...
%       lcpred
%   %%
%   %% npsol should have just called model with same value of theta and
%   %% model has commputed grad.
%   %% So pull the gradient out of the global variables
%   theta;
%   retval = grad;
% endfunction

%%
%% perform optimization
%%
theta0 = [vrdiclo; k1k2ratio];
% lb     = [0; 0];
% large  = 1e20;
% ub     = [large;large];

% more off;
% npsol_options ('major print level', 10);
% fcnprec = 1e-5;  % absolute tolerance for ssq error small
% npsol_options ('function precision', fcnprec);
% npsol_options ('optimality_tolerance', 1e-3);

% [theta, obj, info, lambda] = npsol (theta0, 'model')
fcn = 'res_reduced';
%phi = model(theta0)
x0 = [vro; 0];
sx0 = zeros(2,2);
ddasac_options ('initial step size', 1e-6);
[x, xdot, sx, sxdot] = ddasac (fcn, x0, theta0, sx0, timeout, tcrit);
%%
%% express the other states in terms of vr and eps2
%%
vr   = x(:,1);
eps2 = x(:,2);
badded = (vr-vro)*teaf;
eps1 = badded - eps2;
ca = (vrdiclo-eps1)./vr;
cb = badded-(eps1+eps2);  % sanity check, should be zero
cc = (eps1 - eps2)./vr;
ccd = eps2./vr;
lcpred = cc./(cc +2*ccd);
%% avoid 0/0 at first two time points;
lcpred(1:2) = 1;
%plot (tlc, lcmeas, '*', timeout, lcpred)
table1 = [timeout lcpred];
table2 = [tlc lcmeas];
save -ascii bvsm_red_ex_pred.dat table1;
save -ascii bvsm_red_ex_data.dat table2;
%%
%%compute objective, phi, and gradient, grad, and Gauss-Newton Hessian
%%approximation, Hgn, for the sum squares prediction error
%%
lc_at_tlc = lcpred(idx);
resid = lcmeas - lc_at_tlc;
phi  = resid'*resid;
%% compute lc sensitivity; set to zero until B is added
dlcdeps2 = -2./badded;
dlcdeps2( (badded == 0)) = 0;
sx1 = nth(sx,1);
lc1 = sx1(:,2).*dlcdeps2;
sx2 = nth(sx,2);
lc2 = sx2(:,2).*dlcdeps2;

%%
%%compute objective, phi, and gradient, grad, and Gauss-Newton Hessian
%%approximation, Hgn, for the sum squares prediction error
%%
lc1_at_tlc = lc1(idx);
lc2_at_tlc = lc2(idx);

sens = [lc1_at_tlc lc2_at_tlc];

grad = -2 * sens' * resid;
Hgn  = 2 * sens'*sens;
ndata = size(tlc,1);
np = 2;
alpha = 0.95;
Fstat = np*finv(alpha,np,ndata-np);
%%
%% fix me, need the F-statistic here, approximate with chisq
%%
%chisq = 5.99;  % alpha=0.95, np=2;
a     = 2*resid'*resid/(ndata-np)*Fstat; 
conf  = sqrt(diag(inv(Hgn))*a)
[xx, yy, major, minor, bbox] = ellipse (Hgn, a, 100, theta0);
tmp = [xx, yy];
save -ascii bvsm_red_ex-outline.dat tmp;
save -ascii bvsm_red_ex-bbox.dat bbox;

plot (table1(:,1), table1(:,2), table2(:,1), table2(:,2), '+');
title ('Figure 9.36')

res_reduced.m

function xdot = res_reduced(x,time,theta)
  global vro teaf flowrate tflow vro
  %%	Passed variable                  Local equivalent
  %%
  %%	x(1)                             reactor contents volume
  %%	x(2)                             extent of second reaction
  %%	
  vrdiclo = theta(1);
  k1k2ratio = theta(2);
  vr = x(1);
  if ( vr <= 0.0 )
    error ('(res_reduced).... vr is zero, time= %d',time);
  end
  eps2   = x(2);
  badded = (vr-vro)*teaf;
  eps1   = badded - eps2;
  %%
  %% look up flowrate at current time 
  %%
  [index, is_switching_time] = find_time_index (tflow, time);
  fteaf   = flowrate(index) * teaf;
  %%
  %%	calculate xdots
  %%
  xdot(1) = flowrate(index);
  if (badded == 0.)
    xdot(2) = 0.;
  else
    xdot(2) = fteaf/( 1 + k1k2ratio*(vrdiclo - (badded-eps2) )   ...
	     /(badded-2*eps2) );
  end

model.m

function phi = model(theta)
  global vrdiclo k1k2ratio vro teaf flowrate tflow vro timeout idx ...
      lcmeas fcn x xdot sx sxdot x0 sx0 tcrit
  vrdiclo = theta(1);
  k1k2ratio = theta(2);
  %%
  %% find the solution of the reduced model at timeout
  %%
  x0 = [vro; 0];
  sx0 = zeros(2,2);
  [x, xdot, sx, sxdot] = ddasac (fcn, x0, theta, sx0, timeout, tcrit)
  %%
  %% express the other states in terms of vr and eps2
  %%
  vr   = x(:,1);
  eps2 = x(:,2);
  badded = (vr-vro)*teaf;
  eps1 = badded - eps2;
  ca = (vrdiclo-eps1)./vr;
  cb = badded-(eps1+eps2);  % sanity check, should be zero
  cc = (eps1 - eps2)./vr;
  ccd = eps2./vr;
  lcpred = cc./(cc +2*ccd);
  lc_at_tlc = lcpred(idx);
  %%table = [timeout vr ca cb cc ccd lcpred eps1 eps2];
  %%
  %%compute objective, phi, and gradient, grad, and Gauss-Newton Hessian
  %%approximation, Hgn, for the sum squares prediction error
  %%
  resid = lcmeas - lc_at_tlc;
  theta
  phi  = resid'*resid;

../../util/common/ellipse.m

%% Copyright (C) 2001, James B. Rawlings and John W. Eaton
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

%% [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)
%%
%% Given a 2x2 matrix, generate ellipse data for plotting.  The
%% arguments N and SHIFT are optional.  If N is an empty matrix, a
%% default value of 100 is used.

function [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)

  if (nargin < 3)
    n = 100;
  end

  if (isempty (n))
    n = 100;
  end

  if (nargin < 4)
    shift = [0, 0];
  end

  ss = size (shift);

  if (any (ss ~= [1, 2]))
    if (ss == [2, 1])
      shift = shift';
    else
      error ('shift must be a 2-element row vector');
    end
  end

  if (nargin > 1)

    [v, l] = eig (amat / level);

    dl = diag(l);
    if (any (imag (dl)) || any (dl <= 0))
      error ('ellipse: amat must be positive definite');
    end

    %% Generate contour data.

    a = 1 / sqrt (l(1,1));
    b = 1 / sqrt (l(2,2));

    t = linspace (0, 2*pi, n)';

    xt = a * cos (t);
    yt = b * sin (t);

    %% Rotate the contours.

    ra = atan2 (v(2,1), v(1,1));

    cos_ra = cos (ra);
    sin_ra = sin (ra);

    x = xt * cos_ra - yt * sin_ra + shift(1);
    y = xt * sin_ra + yt * cos_ra + shift(2);

    %% Endpoints of the major and minor axes.

    minor = (v * diag ([a, b]))';
    major = minor;

    major(2,:) = -major(1,:);
    minor(1,:) = -minor(2,:);

    t = [1; 1] * shift;

    major = major + t;
    minor = minor + t;

    %% Bounding box for the ellipse using magic formula.

    ainv = inv (amat);
    xbox = sqrt (level * ainv(1,1));
    ybox = sqrt (level * ainv(2,2));

    bbox = [xbox ybox; xbox -ybox; -xbox -ybox; -xbox ybox; xbox ybox];

    t = [1; 1; 1; 1; 1] * shift;
    bbox = bbox + t;

  else
    error ('usage: ellipse (amat, level, n, shift)');
  end

%% endfunction

Data files

bvsm_red-bbox.dat

 2.35678015e+00 2.28802081e+00
 2.35678015e+00 1.71197919e+00
 2.34261985e+00 1.71197919e+00
 2.34261985e+00 2.28802081e+00
 2.35678015e+00 2.28802081e+00

bvsm_red-outline.dat

 2.35505856e+00 1.71197920e+00
 2.35534126e+00 1.71256454e+00
 2.35560126e+00 1.71430728e+00
 2.35583749e+00 1.71720041e+00
 2.35604901e+00 1.72123226e+00
 2.35623496e+00 1.72638662e+00
 2.35639460e+00 1.73264272e+00
 2.35652729e+00 1.73997538e+00
 2.35663248e+00 1.74835506e+00
 2.35670975e+00 1.75774803e+00
 2.35675880e+00 1.76811646e+00
 2.35677943e+00 1.77941861e+00
 2.35677155e+00 1.79160896e+00
 2.35673520e+00 1.80463842e+00
 2.35667052e+00 1.81845454e+00
 2.35657777e+00 1.83300168e+00
 2.35645733e+00 1.84822126e+00
 2.35630967e+00 1.86405200e+00
 2.35613540e+00 1.88043015e+00
 2.35593522e+00 1.89728977e+00
 2.35570994e+00 1.91456297e+00
 2.35546045e+00 1.93218019e+00
 2.35518776e+00 1.95007050e+00
 2.35489298e+00 1.96816186e+00
 2.35457729e+00 1.98638142e+00
 2.35424196e+00 2.00465581e+00
 2.35388834e+00 2.02291146e+00
 2.35351786e+00 2.04107485e+00
 2.35313200e+00 2.05907285e+00
 2.35273233e+00 2.07683298e+00
 2.35232044e+00 2.09428374e+00
 2.35189800e+00 2.11135484e+00
 2.35146671e+00 2.12797756e+00
 2.35102831e+00 2.14408496e+00
 2.35058456e+00 2.15961218e+00
 2.35013725e+00 2.17449670e+00
 2.34968817e+00 2.18867858e+00
 2.34923915e+00 2.20210072e+00
 2.34879198e+00 2.21470907e+00
 2.34834846e+00 2.22645286e+00
 2.34791039e+00 2.23728481e+00
 2.34747952e+00 2.24716129e+00
 2.34705760e+00 2.25604255e+00
 2.34664631e+00 2.26389281e+00
 2.34624733e+00 2.27068047e+00
 2.34586224e+00 2.27637819e+00
 2.34549261e+00 2.28096303e+00
 2.34513992e+00 2.28441654e+00
 2.34480559e+00 2.28672480e+00
 2.34449097e+00 2.28787852e+00
 2.34419732e+00 2.28787306e+00
 2.34392583e+00 2.28670843e+00
 2.34367759e+00 2.28438933e+00
 2.34345360e+00 2.28092510e+00
 2.34325477e+00 2.27632968e+00
 2.34308189e+00 2.27062157e+00
 2.34293565e+00 2.26382377e+00
 2.34281665e+00 2.25596365e+00
 2.34272537e+00 2.24707284e+00
 2.34266218e+00 2.23718717e+00
 2.34262732e+00 2.22634642e+00
 2.34262094e+00 2.21459426e+00
 2.34264307e+00 2.20197800e+00
 2.34269361e+00 2.18854845e+00
 2.34277237e+00 2.17435968e+00
 2.34287902e+00 2.15946882e+00
 2.34301313e+00 2.14393584e+00
 2.34317418e+00 2.12782328e+00
 2.34336149e+00 2.11119602e+00
 2.34357434e+00 2.09412102e+00
 2.34381184e+00 2.07666702e+00
 2.34407306e+00 2.05890431e+00
 2.34435694e+00 2.04090442e+00
 2.34466233e+00 2.02273981e+00
 2.34498800e+00 2.00448364e+00
 2.34533265e+00 1.98620942e+00
 2.34569488e+00 1.96799073e+00
 2.34607324e+00 1.94990093e+00
 2.34646621e+00 1.93201286e+00
 2.34687219e+00 1.91439854e+00
 2.34728957e+00 1.89712892e+00
 2.34771664e+00 1.88027352e+00
 2.34815171e+00 1.86390022e+00
 2.34859301e+00 1.84807494e+00
 2.34903876e+00 1.83286141e+00
 2.34948718e+00 1.81832090e+00
 2.34993646e+00 1.80451193e+00
 2.35038478e+00 1.79149013e+00
 2.35083035e+00 1.77930793e+00
 2.35127137e+00 1.76801437e+00
 2.35170605e+00 1.75765494e+00
 2.35213266e+00 1.74827135e+00
 2.35254948e+00 1.73990137e+00
 2.35295482e+00 1.73257873e+00
 2.35334705e+00 1.72633289e+00
 2.35372460e+00 1.72118902e+00
 2.35408595e+00 1.71716782e+00
 2.35442963e+00 1.71428548e+00
 2.35475427e+00 1.71255362e+00
 2.35505856e+00 1.71197920e+00

bvsm_red_data.dat

 4.14000000e+02 1.70500000e-01
 4.24000000e+02 1.60400000e-01
 4.34000000e+02 1.32500000e-01
 4.44000000e+02 1.08400000e-01
 4.93000000e+02 1.01400000e-01
 5.03000000e+02 1.04200000e-01
 5.13000000e+02 1.04500000e-01
 5.23000000e+02 9.70000000e-02
 5.33000000e+02 8.24000000e-02
 5.82000000e+02 6.78000000e-02
 5.92000000e+02 6.97000000e-02
 6.02000000e+02 7.02000000e-02
 6.12000000e+02 6.92000000e-02
 6.21000000e+02 6.84000000e-02
 6.31000000e+02 6.75000000e-02
 6.41000000e+02 6.73000000e-02
 6.51000000e+02 6.73000000e-02
 6.61000000e+02 6.54000000e-02
 6.71000000e+02 5.77000000e-02
 6.81000000e+02 5.71000000e-02
 6.90000000e+02 5.63000000e-02
 7.00000000e+02 4.66000000e-02
 7.10000000e+02 4.55000000e-02
 7.20000000e+02 4.49000000e-02
 7.30000000e+02 4.56000000e-02
 7.40000000e+02 4.50000000e-02
 7.50000000e+02 4.34000000e-02
 7.59000000e+02 4.28000000e-02
 7.69000000e+02 4.26000000e-02
 7.79000000e+02 4.09000000e-02
 7.89000000e+02 3.99000000e-02
 7.99000000e+02 3.97000000e-02
 8.48000000e+02 3.94000000e-02
 8.58000000e+02 3.74000000e-02
 8.68000000e+02 3.71000000e-02

bvsm_red_pred.dat

 0.00000000e+00 1.00000000e+00
 9.00000000e+00 9.99999798e-01
 1.10000000e+01 9.99610513e-01
 1.90000000e+01 9.98062790e-01
 2.20000000e+01 9.90392655e-01
 2.90000000e+01 9.72458200e-01
 3.30000000e+01 9.62406626e-01
 3.90000000e+01 9.47297428e-01
 4.40000000e+01 9.34904852e-01
 4.90000000e+01 9.22487195e-01
 5.50000000e+01 9.07844745e-01
 5.90000000e+01 8.98063606e-01
 6.60000000e+01 8.81312885e-01
 6.90000000e+01 8.74120140e-01
 7.70000000e+01 8.55329392e-01
 7.90000000e+01 8.50622926e-01
 8.80000000e+01 8.40490459e-01
 8.90000000e+01 8.39363608e-01
 9.90000000e+01 8.16778643e-01
 1.09000000e+02 7.88947246e-01
 1.10000000e+02 7.86218186e-01
 1.19000000e+02 7.61605720e-01
 1.21000000e+02 7.56279021e-01
 1.29000000e+02 7.34930066e-01
 1.32000000e+02 7.27048027e-01
 1.39000000e+02 7.08621353e-01
 1.43000000e+02 6.98296577e-01
 1.49000000e+02 6.82781039e-01
 1.54000000e+02 6.70073222e-01
 1.59000000e+02 6.57342993e-01
 1.65000000e+02 6.42326041e-01
 1.69000000e+02 6.32297920e-01
 1.76000000e+02 6.15117757e-01
 1.79000000e+02 6.07742907e-01
 1.87000000e+02 5.88369493e-01
 1.89000000e+02 5.83518891e-01
 1.98000000e+02 5.62072947e-01
 1.99000000e+02 5.59686547e-01
 2.09000000e+02 5.36272358e-01
 2.19000000e+02 5.13078028e-01
 2.20000000e+02 5.10186146e-01
 2.29000000e+02 4.84117324e-01
 2.31000000e+02 4.78381037e-01
 2.39000000e+02 4.55400975e-01
 2.42000000e+02 4.46868449e-01
 2.49000000e+02 4.26931068e-01
 2.53000000e+02 4.15723098e-01
 2.59000000e+02 3.98889511e-01
 2.64000000e+02 3.85166344e-01
 2.69000000e+02 3.71427300e-01
 2.75000000e+02 3.55130088e-01
 2.79000000e+02 3.44254082e-01
 2.86000000e+02 3.25430625e-01
 2.89000000e+02 3.17356258e-01
 2.97000000e+02 2.96225991e-01
 2.99000000e+02 2.90939561e-01
 3.08000000e+02 2.67676672e-01
 3.09000000e+02 2.65090413e-01
 3.19000000e+02 2.39609837e-01
 3.29000000e+02 2.14446830e-01
 3.30000000e+02 2.13759438e-01
 3.39000000e+02 2.07572868e-01
 3.41000000e+02 2.05184542e-01
 3.49000000e+02 1.95630916e-01
 3.52000000e+02 1.87645768e-01
 3.59000000e+02 1.69015942e-01
 3.63000000e+02 1.66210223e-01
 3.69000000e+02 1.62002003e-01
 3.74000000e+02 1.61991702e-01
 3.79000000e+02 1.61981409e-01
 3.85000000e+02 1.61975310e-01
 3.89000000e+02 1.61971236e-01
 3.96000000e+02 1.61962967e-01
 3.99000000e+02 1.61959433e-01
 4.07000000e+02 1.61950861e-01
 4.09000000e+02 1.61948752e-01
 4.14000000e+02 1.61944754e-01
 4.18000000e+02 1.61941567e-01
 4.19000000e+02 1.61940752e-01
 4.24000000e+02 1.59006944e-01
 4.29000000e+02 1.56073392e-01
 4.34000000e+02 1.43198714e-01
 4.39000000e+02 1.30331460e-01
 4.40000000e+02 1.28314253e-01
 4.44000000e+02 1.20248119e-01
 4.49000000e+02 1.10172146e-01
 4.51000000e+02 1.07893904e-01
 4.59000000e+02 9.87858974e-02
 4.62000000e+02 9.87800639e-02
 4.69000000e+02 9.87664651e-02
 4.73000000e+02 9.87589127e-02
 4.79000000e+02 9.87476144e-02
 4.84000000e+02 9.87389839e-02
 4.89000000e+02 9.87303773e-02
 4.93000000e+02 9.87261085e-02
 4.95000000e+02 9.87239741e-02
 4.99000000e+02 9.87197124e-02
 5.03000000e+02 9.87156333e-02
 5.06000000e+02 9.87125939e-02
 5.09000000e+02 9.87095603e-02
 5.13000000e+02 9.87063856e-02
 5.17000000e+02 9.87032110e-02
 5.19000000e+02 9.87016182e-02
 5.23000000e+02 9.66844232e-02
 5.28000000e+02 9.41635576e-02
 5.29000000e+02 9.36594580e-02
 5.33000000e+02 8.69281373e-02
 5.39000000e+02 7.68418458e-02
 5.49000000e+02 7.66794336e-02
 5.50000000e+02 7.66631956e-02
 5.59000000e+02 7.65170649e-02
 5.61000000e+02 7.64921957e-02
 5.69000000e+02 7.63928571e-02
 5.72000000e+02 7.63408816e-02
 5.79000000e+02 7.62195788e-02
 5.82000000e+02 7.29707135e-02
 5.83000000e+02 7.18881165e-02
 5.89000000e+02 6.53965064e-02
 5.92000000e+02 6.54165410e-02
 5.94000000e+02 6.54298895e-02
 5.99000000e+02 6.54632291e-02
 6.02000000e+02 6.54571191e-02
 6.05000000e+02 6.54509894e-02
 6.09000000e+02 6.54428200e-02
 6.12000000e+02 6.54432288e-02
 6.16000000e+02 6.54437738e-02
 6.19000000e+02 6.54441831e-02
 6.21000000e+02 6.54499038e-02
 6.27000000e+02 6.54670662e-02
 6.29000000e+02 6.54727987e-02
 6.31000000e+02 6.54926853e-02
 6.38000000e+02 6.55622898e-02
 6.39000000e+02 6.55722237e-02
 6.41000000e+02 6.55725257e-02
 6.49000000e+02 6.55736338e-02
 6.51000000e+02 6.55826236e-02
 6.59000000e+02 6.56185678e-02
 6.60000000e+02 6.52411395e-02
 6.61000000e+02 6.48637273e-02
 6.69000000e+02 6.18453056e-02
 6.71000000e+02 6.05433152e-02
 6.79000000e+02 5.53385281e-02
 6.81000000e+02 5.53170103e-02
 6.82000000e+02 5.53062698e-02
 6.89000000e+02 5.52311164e-02
 6.90000000e+02 5.52225365e-02
 6.93000000e+02 5.51968410e-02
 6.99000000e+02 5.51454183e-02
 7.00000000e+02 5.40459597e-02
 7.04000000e+02 4.96506455e-02
 7.09000000e+02 4.41625253e-02
 7.10000000e+02 4.41617964e-02
 7.15000000e+02 4.41581339e-02
 7.19000000e+02 4.41552084e-02
 7.20000000e+02 4.41545303e-02
 7.26000000e+02 4.41504612e-02
 7.29000000e+02 4.41484745e-02
 7.30000000e+02 4.41482142e-02
 7.37000000e+02 4.41465046e-02
 7.39000000e+02 4.41460314e-02
 7.40000000e+02 4.41458957e-02
 7.48000000e+02 4.41448095e-02
 7.49000000e+02 4.41446645e-02
 7.50000000e+02 4.41439320e-02
 7.59000000e+02 4.41373371e-02
 7.69000000e+02 4.18653306e-02
 7.70000000e+02 4.18638962e-02
 7.79000000e+02 4.18509648e-02
 7.81000000e+02 4.15590507e-02
 7.89000000e+02 4.03914875e-02
 7.92000000e+02 4.03858704e-02
 7.99000000e+02 4.03727815e-02
 8.03000000e+02 4.03653023e-02
 8.09000000e+02 4.03540835e-02
 8.14000000e+02 4.03447347e-02
 8.19000000e+02 4.03353912e-02
 8.25000000e+02 4.03251358e-02
 8.29000000e+02 4.03183098e-02
 8.36000000e+02 4.03074825e-02
 8.39000000e+02 4.03028584e-02
 8.47000000e+02 4.02930981e-02
 8.48000000e+02 4.02918785e-02
 8.49000000e+02 4.02906693e-02
 8.58000000e+02 4.02816436e-02
 8.59000000e+02 4.02806321e-02
 8.68000000e+02 4.02706316e-02
 8.69000000e+02 4.02695243e-02

lc.dat

    414        0.17050
    424        0.16040
    434        0.13250
    444        0.10840
    493        0.10140
    503        0.10420
    513        0.10450
    523        0.09700
    533        0.08240
    582        0.06780
    592        0.06970
    602        0.07020
    612        0.06920
    621        0.06840
    631        0.06750
    641        0.06730
    651        0.06730
    661        0.06540
    671        0.05770
    681        0.05710
    690        0.05630
    700        0.04660
    710        0.04550
    720        0.04490
    730        0.04560
    740        0.04500
    750        0.04340
    759        0.04280
    769        0.04260
    779        0.04090
    789        0.03990
    799        0.03970
    848        0.03940
    858        0.03740
    868        0.03710

flow.dat

   9.000000      9.7817838E-02   4.786000    
   19.00000       1.290601       14.21800    
   29.00000       1.262779       13.99800    
   39.00000       1.240016       13.81800    
   49.00000       1.215736       13.62600    
   59.00000       1.187156       13.40000    
   69.00000       1.160599       13.19000    
   79.00000      0.5545961       8.398001    
   89.00000       1.109508       12.78600    
   99.00000       1.361925       14.78200    
   109.0000       1.332333       14.54800    
   119.0000       1.294648       14.25000    
   129.0000       1.271885       14.07000    
   139.0000       1.244569       13.85400    
   149.0000       1.220794       13.66600    
   159.0000       1.197778       13.48400    
   169.0000       1.170463       13.26800    
   179.0000       1.150988       13.11400    
   189.0000       1.128983       12.94000    
   199.0000       1.105968       12.75800    
   209.0000       1.092563       12.65200    
   219.0000       1.360155       14.76800    
   229.0000       1.344473       14.64400    
   239.0000       1.329045       14.52200    
   249.0000       1.305524       14.33600    
   259.0000       1.275426       14.09800    
   269.0000       1.259239       13.97000    
   279.0000       1.244063       13.85000    
   289.0000       1.219783       13.65800    
   299.0000       1.191961       13.43800    
   309.0000       1.173751       13.29400    
   319.0000       1.158322       13.17200    
   329.0000      0.3163429       6.514000    
   339.0000      0.5495375       8.358001    
   349.0000       1.224841       13.69800    
   359.0000      0.3228684       6.565600    
   369.0000      9.4845297E-04   4.020000    
   379.0000      4.6794413E-04   4.016201    
   389.0000      5.4383278E-04   4.016800    
   399.0000      4.9326418E-04   4.016400    
   409.0000      3.6680698E-04   4.015400    
   419.0000      0.2701340       6.148601    
   429.0000       1.185891       13.39000    
   439.0000      0.9299334       11.36600    
   449.0000      0.5259147       8.171202    
   459.0000      8.9793204E-04   4.019600    
   469.0000      8.7256433E-04   4.019400    
   479.0000      7.9672335E-04   4.018800    
   489.0000      4.9319270E-04   4.016400    
   499.0000      4.6799183E-04   4.016200    
   509.0000      3.6678315E-04   4.015400    
   519.0000      0.2330809       5.855600    
   529.0000      0.7784327       10.16800    
   539.0000      7.5244429E-03   4.072000    
   549.0000      7.5244186E-03   4.072000    
   559.0000      5.7539940E-03   4.058000    
   569.0000      8.0302954E-03   4.076001    
   579.0000      0.5019882       7.982001    
   589.0000     -3.0983209E-03   3.988000    
   599.0000      9.4850064E-04   4.020000    
   609.0000     -6.3252446E-05   4.012000    
   619.0000     -1.3278484E-03   4.002000    
   629.0000     -4.6158792E-03   3.976000    
   639.0000     -6.3204767E-05   4.012000    
   649.0000     -2.0865917E-03   3.996000    
   659.0000      0.1752122       5.398000    
   669.0000      0.3024322       6.404000    
   679.0000      4.9952744E-03   4.052001    
   689.0000      3.9835931E-03   4.044001    
   699.0000      0.5114223       8.056601    
   709.0000      3.4151078E-04   4.015200    
   719.0000      3.1619071E-04   4.015000    
   729.0000      1.1386871E-04   4.013400    
   739.0000      6.3300133E-05   4.013000    
   749.0000      3.4151078E-04   4.015200    
   759.0000      0.1059619       4.850400    
   769.0000      6.7019463E-04   4.017800    
   779.0000      6.8099439E-02   4.551000    
   789.0000      8.7261200E-04   4.019400    
   799.0000      9.7372534E-04   4.020200    
   809.0000      8.7256433E-04   4.019400    
   819.0000      7.9669955E-04   4.018800    
   829.0000      7.2083471E-04   4.018200    
   839.0000      5.6912901E-04   4.017000    
   849.0000      4.6801567E-04   4.016200    
   859.0000      5.1856041E-04   4.016600    
   869.0000      4.4267176E-04   4.016000    

lcsim.dat

 9 0.999999811053562 1.01722259342699
 11 0.999587658944779 0.99385414630478
 19 0.997925636447302 1.0063033273586
 22 0.989776227519767 0.973524345205085
 29 0.97073863708603 0.951277459498521
 33 0.96009379707585 0.955916127193075
 39 0.944126388142496 0.933176658461481
 44 0.931059995373657 0.912129208525589
 49 0.917993531731551 0.924765345958655
 55 0.902620794363343 0.929423720426881
 59 0.892372164807407 0.895315152189819
 66 0.874858939641787 0.878874365681482
 69 0.867353293988678 0.887610914800141
 77 0.847785922726378 0.848974596028909
 79 0.842894103481712 0.848985755994263
 88 0.83237481405589 0.829011613033741
 89 0.831205888518126 0.828704133459361
 99 0.807823351175219 0.808888562620313
 109 0.779121310864896 0.786381406938047
 110 0.77631345417554 0.776252149913632
 119 0.751042820774906 0.756749889814251
 121 0.745585969175354 0.732134962762848
 129 0.723758566435087 0.723820015722546
 132 0.715717201841988 0.725510391009964
 139 0.696954020548539 0.697223175175415
 143 0.68646246764278 0.67349142047977
 149 0.670725076033402 0.668328166140843
 154 0.657861133646678 0.662317372691344
 159 0.644997187699813 0.657514340724487
 165 0.629851520352139 0.621581418804898
 169 0.61975441688066 0.626779273090415
 176 0.602487407389775 0.585884236756459
 179 0.595087253950168 0.597729143403102
 187 0.575681893808323 0.580350965840298
 189 0.570830568007055 0.572046209133092
 198 0.549416907668722 0.541541840673101
 199 0.547037560114431 0.548612747410106
 209 0.523729613746048 0.538541249955535
 219 0.500704189037897 0.513368588184931
 220 0.497837706149487 0.50676216119614
 229 0.472039332645793 0.477205527084727
 231 0.466372465615925 0.488838142877278
 239 0.443705011167859 0.437219925680493
 242 0.435302251160462 0.432986828416903
 249 0.415695826849098 0.416071875495965
 253 0.404690424332785 0.401912894973921
 259 0.388182323384162 0.39285382966983
 264 0.37474272150017 0.378923137070906
 269 0.361303120788802 0.367484106358756
 275 0.345380258106417 0.338508666896052
 279 0.33476502597149 0.333439320473033
 286 0.316412252844235 0.313210164070985
 289 0.308546738752642 0.302731294760027
 297 0.287981499998565 0.284780383973117
 299 0.282840191132878 0.274816315777158
 308 0.260231996792122 0.273443081396385
 309 0.257719953569755 0.251919518622741
 319 0.232983524949703 0.231578003318462
 329 0.208572150752592 0.225291260764646
 330 0.207905468886441 0.22729540104205
 339 0.201905298645643 0.206577487404493
 341 0.199589040875074 0.210143334269163
 349 0.190323859821764 0.195116540264097
 352 0.182579938375238 0.187682610356096
 359 0.164510700875708 0.151127109700629
 363 0.161788969441419 0.158451310797697
 369 0.15770634207812 0.162343283034236
 374 0.157696338778412 0.149316304212486
 379 0.157686348622417 0.164821037715053
 385 0.15768043157753 0.15766516626921
 389 0.157676482590873 0.161593237909514
 396 0.157668465879349 0.170786320474533
 399 0.157665027844229 0.153198631274023