function [fxind, tid] = rind(BIG,Ex,xc,Nt1,NIT1,speed1,indI,B_lo,B_up)
%RIND computes  E[Jacobian*Indicator|Condition ]*f_{Xc}(xc(:,ix)) 
% where
%     "Indicator" = I{ H_lo(i) < X(i) < H_up(i), i=1:N_t+N_d }
%     "Jacobian"  = J(X(Nt+1),...,X(Nt+Nd+Nc)), special case is 
%     "Jacobian"  = |X(Nt+1)*...*X(Nt+Nd)|=|Xd(1)*Xd(2)..Xd(Nd)|
%     "condition" = Xc=xc(:,ix),  ix=1,...,Nx.
%     X = [Xt; Xd; Xc], a stochastic vector of Multivariate Gaussian 
%         variables where Xt,Xd and Xc have the length Nt, Nd and Nc,
%         respectively. (Recommended limitations Nx,Nt<=100, Nd<=6 and Nc<=10) 
%
% CALL [E, tid] = rind(S,m,xc,[Nt Nj],NIT,speed,indI,B_lo,B_up);
%
%      tid = time taken for the fortran executable to integrate
%        E = expectation/density as explained above size 1 x Nx
%        S = Covariance matrix of X=[Xt;Xd;Xc] size Ntdc x Ntdc (Ntdc=Nt+Nd+Nc)
%        m = the expectation of X=[Xt;Xd;Xc]   size N x 1
%       xc = values to condition on            size Nc x Nx
%       Nt = size of Xt
%       Nj = # of varibles of Xt integrated directly (default Nj=0)
%      NIT = 0,1,2..., dimension of numerical integration (default NIT=1)
%               -1 Integrate all by SADAPT for Ndim<9 and by KRBVRC otherwise 
%               -2 Integrate all by SADAPT by Genz (1992) (Fast)
%               -3 Integrate all by KRBVRC by Genz (1993) (Fast)
%               -4 Integrate all by KROBOV by Genz (1992) (Fast)
%               -5 Integrate all by RCRUDE by Genz (1992)
%               -6 Integrate all by RLHCRUDE using MLHD and center of cell 
%               -7 Integrate all by RLHCRUDE using  LHD and center of cell
%               -8 Integrate all by RLHCRUDE using MLHD and random point within
%               the cell 
%               -9 Integrate all by RLHCRUDE using LHD and random point within the
%               cell
%    speed = defines accuraccy of calculations by choosing different 
%               parameters, possible values: 1,2...,9 (9 fastest, default 4).
%     indI = vector of indices to the different barriers in the  
%            indicator function,  length NI, where   NI = Nb+1 
%            (NB! restriction  indI(1)=0, indI(NI)=Nt+Nd )
%B_lo,B_up = Lower and upper barriers used to compute the integration 
%            limits, H_lo and H_up, respectively. size  Mb x Nb 
%  
%  This routine also uses global data to transport the inputs:  
%  NB!  INITDATA may be called  before calling RIND
%  to initialize some global data
%
%  If  Mb<Nc+1 then B_lo(Mb+1:Nc+1,:) is assumed to be zero. The relation 
%  to the integration limits  are as follows
%
%              H_lo(i)=B_lo(1,j)+B_lo(2:Nc+1,j).'*xc(:,ix), 
%              H_up(i)=B_up(1,j)+B_up(2:Nc+1,j).'*xc(:,ix), 
%
%  where i=indI(j-1)+1:indI(j), j=2:NI, ix=1:Nx
%      
% Examples:% Compute the probability that X1<0,X2<0,X3<0,X4<0,X5<0,
%          % Xi are zeromean Gaussian variables with variances one
%          % and correlations Cov(X(i),X(i+1))=0.3 and zero otherwise:
%          % indI=[0 5], and barriers B_lo=[-inf 0], B_lo=[0  inf]     
%          % gives H_lo = [-inf -inf -inf -inf -inf]  H_lo = [0 0 0 0 0] 
%
%   N = 5; rho=0.3; NIT=3; Nt=N; indI=[0 N];
%   B_lo=-10; B_up=0; m=1.2*ones(N,1);
%   Sc=(ones(N)-eye(N))*rho+eye(N);
%   E = rind(Sc,m,[],Nt,NIT,[],indI,B_lo,B_up) % exact prob. 0.00195
%
% See also: initdata

% References
% Podgorski et al. (1999)
% "Exact distributions for apparent waves in irregular seas"
% Ocean Engineering 
%
%
% Tested on:  Matlab 5.3,DIGITAL UNIX Fortran90 compiler
% History:
% revised ir 10.9 2001 File closed after call to Fortran (last fclose)
% revised jr 16.2 2001  (i) speed added in the call of this routine
%                      (ii) The first four lines after the declaration of 
%                           global variables added (call of initdata, speed)
% revised ir 6.2.2001 Adapted to MATLAB 6
% revised jr 1.2.2001 The definition of Nb, Mb changed again. 
% revised ir 15.11.2000 A bug fixed in the definition of Nb Mb.
% revised ir 10.11.2000 Compilation with rind60.f
% revised jr 09.11.2000 The call to init_data replaced with initdata
%% revised pab 23.05.2000 replaced the matlab call with a call to the F90
%         program RINDD.exe
% revised by Per A. Brodtkorb 14.05.1999
%       - No limitations on size of the inputs
%       - enabled recursive calls
%       - fixed some bugs
%       - added some additonal checks
%       - updated to fortran90 
% by  Igor Rychlik 29.10.1998 (PROGRAM RIND11 --- Version 1.0)

global Nt Nj Nd Nc Nx Ntd Ntdc Nstoc Ndet indXt SQ Mb NI SCIS N_int
global NIT index1 xedni fxcEPSS  H_lo H_up EPSS EPS2 xCutOff SP 


if nargin<7|isempty(speed1)
  speed1=4;
end			    
initdata(speed1);

if isempty(SCIS), SCIS=0;end
if nargin<5|isempty(NIT1)
  NIT=2;
elseif NIT1<0
  SCIS=min(abs(NIT1),10);
  NIT=1;
else
  NIT = NIT1;
  SCIS = 0;
end

Nt=Nt1(1);
if length(Nt1)==2
 Nj=Nt1(2);
else
 Nj=0;
end
Nj=min(Nj,max(Nt,0)); % make sure Nj<Nt


%if Nt>1
%  [U,S,V] = svd(BIG(1:Nt,1:Nt));
%end

Ntdc=size(BIG,1);
if isempty(xc)
  Nc=0;Nx=1;
  xc=1;
else
  [Nc, Nx]=size(xc);  
end
[Mb, Nb]=size(B_lo);
%[Nb, Mb]=size(B_lo);     % The change in November 2000
NI=length(indI);
Nd=Ntdc-Nt-Nc;
Ntd=Nt+Nd;


if 1,% call fortran procedure instead

disp('   Writing data.')
filename=['BIG.in'];
if exist(filename)
  delete(filename)
end
fid=fopen(filename,'wt');
for ix=1:Ntdc
  fprintf(fid,'%12.10f \n',BIG(ix,:));
end
fclose(fid);

filename=['Ex.in'];
if exist(filename)
  delete(filename)
end
fid=fopen(filename,'wt');
fprintf(fid,'%12.10f \n',Ex);
fclose(fid);

filename=['xc.in'];
if exist(filename)
  delete(filename)
end
fid=fopen(filename,'wt');
fprintf(fid,'%12.10f \n',xc);
fclose(fid);

filename=['indI.in'];
if exist(filename)
  delete(filename)
end
fid=fopen(filename,'wt');
fprintf(fid,'%2.0f \n',indI);
fclose(fid);

filename=['B_lo.in'];
if exist(filename)
  delete(filename)
end
fid=fopen(filename,'wt');
fprintf(fid,'%12.10f \n',B_lo');
fclose(fid);

filename=['B_up.in'];
fid=fopen(filename,'wt');
if exist(filename)
  delete(filename)
end
fid=fopen(filename,'wt');
fprintf(fid,'%12.10f \n',B_up');
fclose(fid);


global speed rateLHD XSPLT RELEPS
if isempty(rateLHD), rateLHD=20;end
if isempty(XSPLT), XSPLT=1.5;end
if isempty(RELEPS), RELEPS=EPSS;end
if exist('sizeinfo.in'), delete sizeinfo.in, end
fid=fopen('sizeinfo.in','wt');
fprintf(fid,'%2.0f \n', speed);
fprintf(fid,'%2.0f \n', Nt);
fprintf(fid,'%2.0f \n', Nd);
fprintf(fid,'%2.0f \n', Nc);
fprintf(fid,'%2.0f \n', NI);
fprintf(fid,'%2.0f \n', Mb);
fprintf(fid,'%2.0f \n', Nx);
fprintf(fid,'%2.0f \n', NIT);
fprintf(fid,'%2.0f \n', Nj);
fprintf(fid,'%2.0f \n', rand*10^9); % put seed  resets it to a different state each time.
fprintf(fid,'%2.0f \n', SCIS);
fprintf(fid,'%2.0f \n', rateLHD);
fprintf(fid,'%12.10f \n', XSPLT);
fprintf(fid,'%12.10f \n', EPSS);    % absolute error
fprintf(fid,'%12.10f \n', EPS2); % 
fprintf(fid,'%12.10f \n', xCutOff); % truncation value
fprintf(fid,'%12.10f \n', RELEPS);  % relativ error
fprintf(fid,'%2.0f \n', N_int(1));
fprintf(fid,'%2.0f \n', 1); % minQnr
fprintf(fid,'%2.0f \n', N_int(1)); % Le2Qnr
fclose(fid);

disp('   Starting Fortran executable.')

tic;
dos([wafoexepath, 'rindd70.exe']);
tid=toc;
fxind=load('rind.out');
return
end