function [tout, y, d, p, err, errmsg] = OsiSim(t,u,ConNo);
%
% OsiSim.m		Simulate OLGA pipeline.   
%
% The calling syntax is:
%   1) Integrate to Tend already stored in OlGA Server (set by OsiSetTend), 
%      with the inputs stored in OLGA Server (set by OsiSetInputData).
%      [t, y, d, p, err, errmsg] = OsiSim;
%   2) Integrate to Tend with the inputs stored in the OLGA Server  
%      (set by OsiSetInputData).
%      [t, y, d, p, err, errmsg] = OsiSim(Tend);
%      Note if Tend is a vector then the results are stored in the outputs 
%      for each time in Tend.
%   3) Integrate to the time steps given in T with the inputs U,
%      [t, y, d, p, err, errmsg] = OsiSim(T, U);
%      NOTE U must contain the same number of rows as T does. 
%      If T is empty then OsiSim integrates to Tend stored in the 
%      OLGA Server using as input the first row of U. 
%
% Copyright (c) 1998 by ABB Corporate Research.
% All Rights Reserved.
%
% $Revision: 1.0,   5/10-1998 by Kjetil Havre
%            1.01, 27/11-1998 by Kjetil Havre, minor update in the return of profile variables$  
% $Revision: 1.1, 20/7-2005 by Espen Storkaas, ABB-EOP, included user-defined connection number$ 

  global OsiIBODYsetProfileVar OsiIBODYsetTrendVar

  tout = []; y = []; d = []; p = []; Nsteps = 1; 
  NrowsT = []; SetTime = 0; SetInput = 0;

  if (nargin >= 1) & (~isempty(t))
     NrowsT = size(t,1);
     SetTime = 1;
  end

  if NrowsT > 1
     for i=1:NrowsT-1
        if t(i) >= t(i+1)
           err = -1;
           errmsg = 'OsiSim error: The time in the vector T must increase';
           disp(errmsg);
           return
        end
     end 
  end

  th = OsiGetTime(ConNo);

  if SetTime & (th >= t(1) )
     err = -1;
     errmsg = 'OsiSim error: current time is greater or equal to T(1)';
     disp(errmsg);
     return;
  end

  if (nargin >= 2) & (~isempty(u))
     SetInput = 1; Nsteps = size(u,1);
     if NrowsT > Nsteps
        u = [u; ones(NrowsT-Nsteps,1)*u(Nsteps,:)];
        Nsteps = NrowsT;
        [Mu,Nu] = size(u);
     end
     if (NrowsT ~= Nsteps) & (NrowsT ~= 0) 
        err = -1;
        errmsg = 'OsiSim error: Number of rows in T must be equal to number of rows in U';
        disp(errmsg);
        return;
     end
     if NrowsT == 0
        Nsteps = 1;
     end
  end
  if nargin<=2
      ConNo=1;
  end
  
  if (NrowsT > Nsteps) & (nargin < 2)
     Nsteps = NrowsT;
  end
  
  err = 0; 
  ID=sprintf('Con%s',num2str(ConNo));
  NproSecM = max(OsiIBODYsetProfileVar);
  NproVar = size(OsiIBODYsetProfileVar,1); 
  NtrdVar = size(getfield(OsiIBODYsetTrendVar,ID),1); 
  if NproVar == 0
     d = [];
  else 
     d = (0:+1:NproSecM-1)';
  end 
  
  %
  % Allocate memory to tout, y, d and p before simulation
  %
  
  tout = zeros(Nsteps,1);
  y    = zeros(Nsteps,NtrdVar);
  d    = zeros(NproSecM,NproVar); 
  p    = zeros(NproSecM*Nsteps,NproVar);
  %
  % Initialize d.
  %
  for k = 1:NproVar
    d(1:OsiIBODYsetProfileVar(k),k) = (1:+1:OsiIBODYsetProfileVar(k))';
  end

  for i = 1:Nsteps,
          
     if SetTime
        [err, errmsg] = OsiSetTend(t(i),ConNo);
     end
     if SetInput & (err == 0)
        [err, errmsg] = OsiSetInputData(u(i,:)',ConNo);
     end


     if err == 0
       [err, msg, ibody, rbody, errmsg] = osmi(ConNo, 'OLGA_SIM_STEP', [], [], []);
       if err == 0   
          [tp, err, errmsg] = OsiGetTime(ConNo);
          if err == 0
             tout(i,1) = tp(1,1);
          else
             disp(errmsg)
             return   
          end
          y(i,:) = ibody';
          
          nrs = (i-1)*NproSecM+1;

          AccV = 1;
          for k = 1:NproVar
             p( nrs:(nrs+OsiIBODYsetProfileVar(k)-1),k ) = rbody( AccV:AccV+OsiIBODYsetProfileVar(k)-1); 
             AccV = AccV+OsiIBODYsetProfileVar(k);
          end
          if (Nsteps >= 10) & (rem(i, Nsteps/40) == 0)
             disp([num2str(i/Nsteps*100), '% of simulation completed.'])
          end          
       else
          disp(errmsg)
          return
       end       
    end    
  end


