function [pdf,err, options] = cov2mmtpdf(R,dt,u,def,Nstart,hg,options)
%COV2MMTPDF Joint density of Maximum, minimum and period. 
%
% CALL  [pdf, err, options] = cov2mmtpdf(R,dt,u,def,Nstart,hg,options)
%
% pdf     = calculated pdf size Nx x Ntime
% err     = error estimate   
% options = requested and actual rindoptions used in integration.
% R       = [R0,R1,R2,R3,R4] column vectors with autocovariance and its
%          derivatives, i.e., Ri (i=1:4) are vectors with the 1'st to 4'th
%          derivatives of R0.  size Ntime x Nr+1
% dt      = time spacing between covariance samples, i.e., 
%           between R0(1),R0(2).
% u       = crossing level
% def     = integer defining pdf calculated:
%           0 : maximum and the following minimum. (M,m) (default)
%           1 : level v separated Maximum and minimum   (M,m)_v 
%           2 : maximum, minimum and period between (M,m,TMm)
%           3 : level v separated Maximum, minimum and period
%               between (M,m,TMm)_v
%           4 : level v separated Maximum, minimum and the period
%               from Max to level v-down-crossing (M,m,TMd)_v.
%           5 : level v separated Maximum, minimum and the period from 
%               level v-down-crossing to min. (M,m,Tdm)_v 
% Nstart  = index to where to start calculation, i.e., t0 = t(Nstart)
% hg      = vector of amplitudes length Nx or 0
% options = rind options structure defining the integration parameters 
%  
% COV2MMTPDF computes joint density of the  maximum and the following    
% minimum or level u separated maxima and minima + period/wavelength    
%
% For DEF = 0,1 : (Maxima, Minima and period/wavelength)
%         = 2,3 : (Level v separated Maxima and Minima and
%                  period/wavelength between them)
%  
% If Nx==1 then the conditional  density for  period/wavelength between
% Maxima and Minima given the Max and Min is returned
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
% Y=  X'(t2)..X'(ts)..X'(tn-1)||X''(t1) X''(tn)|| X'(t1) X'(tn)  X(t1) X(tn) 
% = [       Xt                   Xd                    Xc            ]  
%                                                                       
% Nt = tn-2, Nd = 2, Nc = 4
%                                                                       
% Xt= contains Nt time points in the indicator function                 
% Xd=    "     Nd    derivatives in Jacobian
% Xc=    "     Nc    variables to condition on                          
%                                                                       
% There are 3 (NI=4) regions with constant barriers:                    
% (indI(1)=0);     for i\in (indI(1),indI(2)]    Y(i)<0.                
% (indI(2)=Nt)  ;  for i\in (indI(2)+1,indI(3)], Y(i)<0 (deriv. X''(t1)) 
% (indI(3)=Nt+1);  for i\in (indI(3)+1,indI(4)], Y(i)>0 (deriv. X''(tn))  
% 
%
% For DEF = 4,5 (Level v separated Maxima and Minima and
%                period/wavelength from Max to crossing)
%
% If Nx==1 then the conditional joint density for  period/wavelength
% between Maxima, Minima and Max to level v crossing given the Max and
% the min is returned
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Y=  X'(t2)..X'(ts)..X'(tn-1)||X''(t1) X''(tn) X'(ts)|| X'(t1) X'(tn)  X(t1) X(tn) X(ts)
% = [       Xt                      Xd                     Xc            ]
%                                                                         
% Nt = tn-2, Nd = 3, Nc = 5
%                                                                         
% Xt= contains Nt time points in the indicator function                      
% Xd=    "     Nd    derivatives                                             
% Xc=    "     Nc    variables to condition on                               
%                                                                            
% There are 4 (NI=5) regions with constant barriers:                        
% (indI(1)=0);     for i\in (indI(1),indI(2)]    Y(i)<0.                    
% (indI(2)=Nt)  ;  for i\in (indI(2)+1,indI(3)], Y(i)<0 (deriv. X''(t1))    
% (indI(3)=Nt+1);  for i\in (indI(3)+1,indI(4)], Y(i)>0 (deriv. X''(tn))
% (indI(4)=Nt+2);  for i\in (indI(4)+1,indI(5)], Y(i)<0 (deriv. X'(ts))    
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  
% History
% Revised pab 22Nov2003
%  -new inputarguments  
%  -added err to output  
%Revised pab 03.04.2003
% -translated from f90 to matlab
%Revised pab 22.04.2000
% - added mean separated min/max + (Tdm, TMd) period distributions
% - added scis 

R0 = R(:,1);
R1 = R(:,2);
R2 = R(:,3);
R3 = R(:,4);
R4 = R(:,5);

Ntime = length(R0);

Nx0  = max(1,length(hg));
Nx1  = Nx0;
%Nx0 = Nx1; % just plain Mm
if (def>1) 
  Nx1 = Nx0/2;
 %  Nx0 = 2*Nx1;   % level v separated max2min densities wanted     
end
%disp(sprintf('def = %d',def))

% ***** The bound 'infinity' is set to 100*sigma *****
XdInf = 100.d0*sqrt(R4(1));
XtInf = 100.d0*sqrt(-R2(1));

Nc = 4;
NI = 4;
Nd = 2; 
Mb = 1; 
Nj = 0;

Nstart = max(2,Nstart); 

symmetry = 0;
isOdd = mod(Nx1,2);
if (def<=1) %THEN % just plain Mm 
   Nx = Nx1*(Nx1-1)/2;
   IJ = (Nx1+isOdd)/2;         
   if (hg(1)+hg(Nx1)==0 &  (hg(IJ)==0 |hg(IJ)+hg(IJ+1)==0) ) 
     symmetry=0;
     disp(' Integration region symmetric')
     % May save Nx1-isOdd integrations in each time step 
     % This is not implemented yet.
     %Nx = Nx1*(Nx1-1)/2-Nx1+isOdd
   end

   % CC = normalizing constant = 1/ expected number of zero-up-crossings of X' 
   CC = 2*pi*sqrt(-R2(1)/R4(1)); 
   %  XcScale = log(CC)
   XcScale = log(2*pi*sqrt(-R2(1)/R4(1)));
else
   % level u separated Mm
   Nx = (Nx1-1)*(Nx1-1);
   if ( abs(u)<=eps & hg(1)+hg(Nx1+1)==0 & (hg(Nx1)+hg(2*Nx1)==0) )
     symmetry=0;
     disp(' Integration region symmetric')
     % Not implemented for DEF <= 3
     %IF (DEF.LE.3) Nx = (Nx1-1)*(Nx1-2)/2 
   end

   if (def>3) %THEN
      Nstart = max(Nstart,3);
      Nc = 5;
      NI = 5; 
      Nd = 3;
   end %ENDIF
   %CC= normalizing constant= 1/ expected number of u-up-crossings of X
   CC = 2*pi*sqrt(-R0(1)/R2(1))*exp(0.5D0*u*u/R0(1)); 
   XcScale = log(2*pi*sqrt(-R0(1)/R2(1))) + 0.5D0*u*u/R0(1);
end %ENDIF

options = setfield(options,'xcscale',XcScale);
opt0 = struct2cell(options);
%opt0 = opt0(1:10);
%seed = [];
%opt0 =  {SCIS,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS,seed,NIT1};



if (Nx>1)
   if ((def==0 | def==2)) %  (M,m) or (M,m)v distribution wanted
      asize = [Nx1,Nx1];
   else           
     % (M,m,TMm), (M,m,TMm)v  (M,m,TMd)v or (M,M,Tdm)v distributions wanted 
      asize = [Nx1,Nx1,Ntime];
   end 
elseif (def>3) %
   % Conditional distribution for (TMd,TMm)v or (Tdm,TMm)v given (M,m)  wanted 
   asize = [1,Ntime,Ntime];
else
  % Conditional distribution for  (TMm) or (TMm)v given (M,m) wanted
   asize = [1,1,Ntime];
end
   
% Initialization
%~~~~~~~~~~~~~~~~~
pdf = zeros(asize);
err = pdf;
BIG  = zeros(Ntime+Nc+1,Ntime+Nc+1);
ex   = zeros(1, Ntime + Nc + 1);          

fxind = zeros(Nx,1);
xc    = zeros(Nc,Nx);
   
indI  = zeros(1,NI);
a_up  = zeros(1,NI-1);
a_lo  = zeros(1,NI-1); 


% INFIN =  INTEGER, array of integration limits flags:  size 1 x Nb   (in)
    %            if INFIN(I) < 0, Ith limits are (-infinity, infinity);
    %            if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)];
    %            if INFIN(I) = 1, Ith limits are [Hlo(I), infinity);
    %            if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)].
INFIN = repmat(0,1,NI-1);
INFIN(3)  = 1;
a_up(1,3) = +XdInf;
a_lo(1,1:2) = [-XtInf -XdInf];
if (def>3)
   a_lo(1,4) = -XtInf;   
end

IJ = 0 ;
if (def<=1) % THEN     % Max2min and period/wavelength
   for I=2:Nx1
      J = IJ+I-1;
      xc(3,IJ+1:J) =  hg(I);
      xc(4,IJ+1:J) =  hg(1:I-1).';
      IJ = J;
   end %do
else
   % Level u separated Max2min
   xc(Nc,:) = u;
   % Hg(1) = Hg(Nx1+1)= u => start do loop at I=2 since by definition we must have:  minimum<u-level<Maximum
   for i=2:Nx1                
      J = IJ+Nx1-1;
      xc(3,IJ+1:J) =  hg(i);              % Max > u
      xc(4,IJ+1:J) =  hg(Nx1+2:2*Nx1).';    % Min < u
      IJ = J;
   end %do   
end %IF
if (def <=3)
   h11 = fwaitbar(0,[],sprintf('Please wait ...(start at: %s)',datestr(now)));

   for Ntd = Nstart:Ntime 
      %Ntd=tn
      Ntdc = Ntd+Nc;
      Nt   = Ntd-Nd;
      indI(2) = Nt;
      indI(3) = Nt+1;
      indI(4) = Ntd;
      % positive wave period 
      BIG(1:Ntdc,1:Ntdc) = covinput(BIG(1:Ntdc,1:Ntdc),R0,R1,R2,R3,R4,Ntd,0);
      [fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),a_lo,a_up,indI, ...
                            xc,Nt,opt0{:});
      
      
      %fxind  = CC*rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),xc,Nt,NIT1,...
      %speed1,indI,a_lo,a_up);

      
      if (Nx<2) %THEN
         % Density of TMm given the Max and the Min. Note that the
         % density is not scaled to unity
         pdf(1,1,Ntd) = fxind(1);
         err(1,1,Ntd) = err0(1).^2;
         %GOTO 100
      else
      
         IJ = 0;
         switch def
         case {-2,-1,0},  % joint density of (Ac,At),(M,m_rfc) or (M,m).
            for i = 2:Nx1 
               J = IJ+i-1;
               pdf(1:i-1,i,1) = pdf(1:i-1,i,1)+fxind(IJ+1:J).'*dt;%*CC;
               err(1:i-1,i,1) = err(1:i-1,i,1)+(err0(IJ+1:J).'*dt).^2;
               IJ = J;
            end %do 
         case  {1}  % joint density of (M,m,TMm)
            for  i = 2:Nx1
               J = IJ+i-1;
               pdf(1:i-1,i,Ntd) = fxind(IJ+1:J).';%*CC
               err(1:i-1,i,Ntd) = (err0(IJ+1:J).').^2;%*CC
               IJ = J;
            end %do
         case {2},  % joint density of level v separated (M,m)v
            for  i = 2:Nx1
               J = IJ+Nx1-1;
               pdf(2:Nx1,i,1) = pdf(2:Nx1,i,1)+fxind(IJ+1:J).'*dt;%*CC;
               err(2:Nx1,i,1) = err(2:Nx1,i,1)+(err0(IJ+1:J).'*dt).^2;
               IJ = J;
            end %do 
         case {3} 
            % joint density of level v separated (M,m,TMm)v
            for  i = 2:Nx1 
               J = IJ+Nx1-1;
               pdf(2:Nx1,i,Ntd) = pdf(2:Nx1,i,Ntd)+fxind(IJ+1:J).';%*CC;
               err(2:Nx1,i,Ntd) = err(2:Nx1,i,Ntd)+(err0(IJ+1:J).').^2;
               IJ = J;
            end %do 
         end % SELECT   
      end %ENDIF
      waitTxt = sprintf('%s Ready: %d of %d',datestr(now),Ntd,Ntime);
      fwaitbar(Ntd/Ntime,h11,waitTxt);
      
   end %do
   close(h11);
   err = sqrt(err);
   %  goto 800
else
   %200 continue
   waitTxt = sprintf('Please wait ...(start at: %s)',datestr(now));
   h11 = fwaitbar(0,[],waitTxt);
   tnold= -1;
   for tn = Nstart:Ntime,
      Ntd  = tn+1;
      Ntdc = Ntd + Nc;
      Nt   = Ntd - Nd;
      indI(2) = Nt;
      indI(3) = Nt + 1;
      indI(4) = Nt + 2;
      indI(5) = Ntd;
      
      if (~symmetry) %IF (SYMMETRY) GOTO 300
        
        for ts = 2:tn-1,
          % positive wave period
          BIG(1:Ntdc,1:Ntdc) = covinput(BIG(1:Ntdc,1:Ntdc),R0,R1,R2,R3, ...
                                        R4,tn,ts,tnold); 
          
          [fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc), ...
                              a_lo,a_up,indI,xc,Nt,opt0{:}); 
          
          %tnold = tn;
          switch (def),
           case {3,4} 
            if (Nx==1) %THEN
                       % Joint density (TMd,TMm) given the Max and the min. 
                       % Note the density is not scaled to unity
                       pdf(1,ts,tn) = fxind(1);%*CC
                       err(1,ts,tn) = err0(1).^2;%*CC
            else
              % 4,  gives level u separated Max2min and wave period
              % from Max to the crossing of level u (M,m,TMd).
              IJ = 0;
              for  i = 2:Nx1,
                J = IJ+Nx1-1; 
                pdf(2:Nx1,i,ts) = pdf(2:Nx1,i,ts)+ fxind(IJ+1:J).'*dt;
                err(2:Nx1,i,ts) = err(2:Nx1,i,ts)+ (err0(IJ+1:J).'*dt).^2;
                IJ = J;
              end %do 
            end
           case  {5}
            if (Nx==1) 
              % Joint density (Tdm,TMm) given the Max and the min. 
              % Note the density is not scaled to unity
              pdf(1,tn-ts+1,tn) = fxind(1); %*CC
              err(1,tn-ts+1,tn) = err0(1).^2;
            else  
              % 5,  gives level u separated Max2min and wave period from 
              % the crossing of level u to the min (M,m,Tdm).
                  
              IJ = 0;
              for  i = 2:Nx1 
                J = IJ+Nx1-1;
                pdf(2:Nx1,i,tn-ts+1)=pdf(2:Nx1,i,tn-ts+1) + fxind(IJ+1:J).'*dt;%*CC;
                err(2:Nx1,i,tn-ts+1)=err(2:Nx1,i,tn-ts+1) + (err0(IJ+1:J).'*dt).^2;
                IJ = J;
              end %do 
            end
          end % SELECT
        end%         enddo
      else % exploit symmetry
         %300   Symmetry
         for ts = 2:floor(Ntd/2)     
           % Using the symmetry since U = 0 and the transformation is
           % linear.
           % positive wave period
            BIG(1:Ntdc,1:Ntdc) = covinput(BIG(1:Ntdc,1:Ntdc),R0,R1,R2,R3, ...
                                          R4,tn,ts,tnold); 
            
            [fxind,err0] = rind(Big(1:Ntdc,1:Ntdc),ex,a_lo,a_up,indI, ...
                                  xc,Nt,opt0{:}); 
            
            %tnold = tn;
            if (Nx==1) % THEN
               % Joint density of (TMd,TMm),(Tdm,TMm) given the max and
               % the min.
               % Note that the density is not scaled to unity
               pdf(1,ts,tn) = fxind(1);%*CC;
               err(1,ts,tn) = err0(1).^2;
               if (ts<tn-ts+1) %THEN
                  pdf(1,tn-ts+1,tn) = fxind(1);%*CC;
                  err(1,tn-ts+1,tn) = err0(1).^2;
               end
               %GOTO 350
            else
              IJ = 0 ;
              switch (def)
               case {4} 
                
                % 4,  gives level u separated Max2min and wave period from 
                % Max to the crossing of level u (M,m,TMd).
                for  i = 2:Nx1 
                  J = IJ+Nx1-1;
                  pdf(2:Nx1,i,ts) = pdf(2:Nx1,i,ts) + fxind(IJ+1:J)*dt;%*CC;
                  err(2:Nx1,i,ts) = err(2:Nx1,i,ts) + (err0(IJ+1:J)*dt).^2;
                  if (ts.LT.tn-ts+1) 
                    % exploiting the symmetry
                    pdf(i,2:Nx1,tn-ts+1) = pdf(i,2:Nx1,tn-ts+1)+fxind(IJ+1:J)*dt;%*CC 
                    err(i,2:Nx1,tn-ts+1) = err(i,2:Nx1,tn-ts+1)+(err0(IJ+1:J)*dt).^2;
                  end
                  IJ = J;
                end %do 
               case {5}
                  % 5,   gives level u separated Max2min and wave period
                  % from the crossing of level u to min (M,m,Tdm).
                  for  i = 2:Nx1, 
                     J = IJ+Nx1-1;
                     
                     pdf(2:Nx1,i,tn-ts+1)=pdf(2:Nx1,i,tn-ts+1)+fxind(IJ+1:J)*dt; 
                     err(2:Nx1,i,tn-ts+1)=err(2:Nx1,i,tn-ts+1)+(err0(IJ+1:J)*dt).^2; 
                     if (ts.LT.tn-ts+1) 
                       %*CC; % exploiting the symmetry
                        pdf(i,2:Nx1,ts) = pdf(i,2:Nx1,ts)+ fxind(IJ+1:J)*dt;
                        err(i,2:Nx1,ts) = err(i,2:Nx1,ts)+ (err0(IJ+1:J)*dt).^2;
                     end %ENDIF
                     IJ = J;
                  end %do 
               end %END SELECT 
            end
            %350     
         end %do
      end
      waitTxt = sprintf('%s Ready: %d of %d',datestr(now),tn,Ntime);
      fwaitbar(tn/Ntime,h11,waitTxt);
      
      %400     print *,'Ready: ',tn,' of ',Ntime
   end %do
   close(h11);
   err = sqrt(err);
end % if      

%Nx1,size(pdf) def  Ntime
if (Nx>1)% THEN
   IJ = 1;
   if (def>2 | def ==1)
      IJ = Ntime;
   end
   pdf = pdf(1:Nx1,1:Nx1,1:IJ);
   err = err(1:Nx1,1:Nx1,1:IJ);
else
   IJ = 1;
   if (def>3)
      IJ = Ntime;
   end
   pdf = squeeze(pdf(1,1:IJ,1:Ntime));
   err = squeeze(err(1,1:IJ,1:Ntime));
end
return
      
   
function BIG = covinput(BIG,R0,R1,R2,R3,R4,tn,ts,tnold)
%COVINPUT Sets up the covariance matrix  
%
% CALL BIG = covinput(BIG, R0,R1,R2,R3,R4,tn,ts)
%   
%  BIG = covariance matrix for X = [Xt,Xd,Xc] in spec2mmtpdf problems.
%     
% The order of the variables in the covariance matrix
% are organized as follows: 
% for  ts <= 1:
%    X'(t2)..X'(ts),...,X'(tn-1) X''(t1),X''(tn)  X'(t1),X'(tn),X(t1),X(tn) 
% = [          Xt               |      Xd       |          Xc             ]
%
% for ts > =2:
%    X'(t2)..X'(ts),...,X'(tn-1) X''(t1),X''(tn) X'(ts)  X'(t1),X'(tn),X(t1),X(tn) X(ts) 
% = [          Xt               |      Xd               |          Xc             ]
%
% where 
%
% Xt= time points in the indicator function
% Xd= derivatives
% Xc=variables to condition on
 
% Computations of all covariances follows simple rules: Cov(X(t),X(s)) = r(t,s),
% then  Cov(X'(t),X(s))=dr(t,s)/dt.  Now for stationary X(t) we have
% a function r(tau) such that Cov(X(t),X(s))=r(s-t) (or r(t-s) will give the same result).
%
% Consequently  Cov(X'(t),X(s))    = -r'(s-t)    = -sign(s-t)*r'(|s-t|)
%               Cov(X'(t),X'(s))   = -r''(s-t)   = -r''(|s-t|)
%               Cov(X''(t),X'(s))  =  r'''(s-t)  = sign(s-t)*r'''(|s-t|)
%               Cov(X''(t),X(s))   =  r''(s-t)   = r''(|s-t|)
% Cov(X''(t),X''(s)) =  r''''(s-t) = r''''(|s-t|)

if nargin<9|isempty(tnold) tnold = -1; end

    
if (ts>1) % THEN
   shft = 1;
   N    = tn + 5 + shft;
   %Cov(Xt,Xc)
   %for
    i = 1:tn-2;
    %j = abs(i+1-ts)
    %BIG(i,N)  = -sign(R1(j+1),R1(j+1)*dble(ts-i-1)) %cov(X'(ti+1),X(ts)) 
    j = i+1-ts;
    tau = abs(j)+1;
    %BIG(i,N)  = abs(R1(tau)).*sign(R1(tau).*j.');
    BIG(i,N)  = R1(tau).*sign(j.');
   %end
   %Cov(Xc)
   BIG(N         ,N) =  R0(1);       % cov(X(ts),X(ts))
   BIG(tn+shft+3 ,N) =  R0(ts);      % cov(X(t1),X(ts))
   BIG(tn+shft+4 ,N) =  R0(tn-ts+1); % cov(X(tn),X(ts))
   BIG(tn+shft+1 ,N) = -R1(ts);      % cov(X'(t1),X(ts))
   BIG(tn+shft+2 ,N) =  R1(tn-ts+1); % cov(X'(tn),X(ts))
   %Cov(Xd,Xc)
   BIG(tn-1 ,N) =  R2(ts);      %cov(X''(t1),X(ts))     
   BIG(tn   ,N) =  R2(tn-ts+1); %cov(X''(tn),X(ts))
   
   %ADD a level u crossing  at ts
   
   %Cov(Xt,Xd)
   %for
      i = 1:tn-2;
      j = abs(i+1-ts);
      BIG(i,tn+shft)  = -R2(j+1); %cov(X'(ti+1),X'(ts)) 
   %end
   %Cov(Xd)  
   BIG(tn+shft,tn+shft) = -R2(1);  %cov(X'(ts),X'(ts))
   BIG(tn-1   ,tn+shft) =  R3(ts); %cov(X''(t1),X'(ts))
   BIG(tn     ,tn+shft) = -R3(tn-ts+1);  %cov(X''(tn),X'(ts))

   %Cov(Xd,Xc)
   BIG(tn+shft ,N       ) =  0.d0;        %cov(X'(ts),X(ts)) 
   BIG(tn+shft,tn+shft+1) = -R2(ts);      % cov(X'(ts),X'(t1))
   BIG(tn+shft,tn+shft+2) = -R2(tn-ts+1); % cov(X'(ts),X'(tn))
   BIG(tn+shft,tn+shft+3) =  R1(ts);      % cov(X'(ts),X(t1))
   BIG(tn+shft,tn+shft+4) = -R1(tn-ts+1); % cov(X'(ts),X(tn))
        
  if (tnold == tn)
     % A previous call to covinput with tn==tnold has been made
     % need only to update  row and column N and tn+1 of big:
     return
     % make lower triangular part equal to upper and then return
     for j=1:tn+shft
        BIG(N,j)      = BIG(j,N)
        BIG(tn+shft,j) = BIG(j,tn+shft)   
     end
     for j=tn+shft+1:N-1
        BIG(N,j) = BIG(j,N)
        BIG(j,tn+shft) = BIG(tn+shft,j)   
     end
     return
  end %if
  tnold = tn;
else
   N = tn+4;
   shft = 0;
end %if
          
if (tn>2)      
   %for
   i=1:tn-2;
   %cov(Xt)
   %   for j=i:tn-2
   %     BIG(i,j) = -R2(j-i+1)              % cov(X'(ti+1),X'(tj+1))
   %  end %do
   
   BIG(i,i) = toeplitz(-R2(i));  % cov(Xt) =   % cov(X'(ti+1),X'(tj+1))


   %cov(Xt,Xc)
   BIG(i      ,tn+shft+1) = -R2(i+1);         %cov(X'(ti+1),X'(t1))  
   BIG(tn-1-i ,tn+shft+2) = -R2(i+1);         %cov(X'(ti+1),X'(tn))  
   BIG(i      ,tn+shft+3) =  R1(i+1);         %cov(X'(ti+1),X(t1))  
   BIG(tn-1-i ,tn+shft+4) = -R1(i+1);         %cov(X'(ti+1),X(tn))  
  
  %Cov(Xt,Xd)
  BIG(i,tn-1)       = R3(i+1);          %cov(X'(ti+1),X''(t1))  
  BIG(tn-1-i,tn)    =-R3(i+1);          %cov(X'(ti+1),X''(tn)) 
  %end %do
end      
%cov(Xd)
BIG(tn-1  ,tn-1  ) = R4(1);
BIG(tn-1  ,tn    ) = R4(tn);     %cov(X''(t1),X''(tn))
BIG(tn    ,tn    ) = R4(1);

%cov(Xc)
BIG(tn+shft+3 ,tn+shft+3) = R0(1);        % cov(X(t1),X(t1))
BIG(tn+shft+3 ,tn+shft+4) = R0(tn);       % cov(X(t1),X(tn))
BIG(tn+shft+1 ,tn+shft+3) = 0.d0 ;        % cov(X(t1),X'(t1))
BIG(tn+shft+2 ,tn+shft+3) = R1(tn);       % cov(X(t1),X'(tn))
BIG(tn+shft+4 ,tn+shft+4) = R0(1) ;       % cov(X(tn),X(tn))
BIG(tn+shft+1 ,tn+shft+4) =-R1(tn);       % cov(X(tn),X'(t1))
BIG(tn+shft+2 ,tn+shft+4) = 0.d0;         % cov(X(tn),X'(tn)) 
BIG(tn+shft+1 ,tn+shft+1) =-R2(1);        % cov(X'(t1),X'(t1))
BIG(tn+shft+1 ,tn+shft+2) =-R2(tn);       % cov(X'(t1),X'(tn))
BIG(tn+shft+2 ,tn+shft+2) =-R2(1);        % cov(X'(tn),X'(tn))
%Xc=X(t1),X(tn),X'(t1),X'(tn) 
%Xd=X''(t1),X''(tn)
%cov(Xd,Xc)
BIG(tn-1  ,tn+shft+3) = R2(1);           %cov(X''(t1),X(t1))     
BIG(tn-1  ,tn+shft+4) = R2(tn);          %cov(X''(t1),X(tn))     
BIG(tn-1  ,tn+shft+1) = 0.d0  ;          %cov(X''(t1),X'(t1))     
BIG(tn-1  ,tn+shft+2) = R3(tn);          %cov(X''(t1),X'(tn))     
BIG(tn    ,tn+shft+3) = R2(tn);          %cov(X''(tn),X(t1))     
BIG(tn    ,tn+shft+4) = R2(1);           %cov(X''(tn),X(tn))     
BIG(tn    ,tn+shft+1) =-R3(tn);          %cov(X''(tn),X'(t1))     
BIG(tn    ,tn+shft+2) = 0.d0;            %cov(X''(tn),X'(tn))
      
     
% make lower triangular part equal to upper 
%for j=1:N-1
%   for i=j+1:N
%      BIG(i,j) = BIG(j,i)
%   end %do
%end %do
lp      = find(tril(ones(size(BIG)),-1)); % indices to lower triangular part
BIGT    = BIG.';
BIG(lp) = BIGT(lp);
return
% END  SUBROUTINE COV_INPUT