Home > wafo > trgauss > spec2mmtpdf.m

spec2mmtpdf

PURPOSE ^

Calculates joint density of Maximum, minimum and period.

SYNOPSIS ^

f = spec2mmtpdf(spec,utc,def,paramt,paramu,options,bound)

DESCRIPTION ^

 SPEC2MMTPDF Calculates joint density of Maximum, minimum and period. 
            
   CALL:  f   = spec2mmtpdf(S,u,def,paramt,paramu,options);  
   
     f    = pdf (density structure) of crests (trough) heights  
     S    = spectral density structure  
     u    = reference level (default the most frequently crossed level). 
    def   = string defining density returned 
             'Mm'   : maximum and the following minimum. (M,m) (default) 
            'rfc'   : maximum and the rainflow minimum height. 
            'AcAt'  : (crest,trough) heights.  
            'vMm'   : level v separated Maximum and minimum   (M,m)_v  
            'MmTMm' : maximum, minimum and period between (M,m,TMm) 
            'vMmTMm': level v separated Maximum, minimum and period 
                      between (M,m,TMm)_v 
            'MmTMd' : level v separated Maximum, minimum and the period 
                      from Max to level v-down-crossing (M,m,TMd)_v. 
            'MmTdm' : level v separated Maximum, minimum and the period from  
                      level v-down-crossing to min. (M,m,Tdm)_v 
            NB! All 'T' above can be replaced by 'L' to get wave length 
            instead.    
  paramt  = [0 tn Nt] defines discretization of half period: tn is the 
            longest period considered while Nt is the number of points, 
            i.e. (Nt-1)/tn is the sampling frequnecy. paramt= [0 10 51] 
            implies that the halfperiods are considered at 51 linearly 
            spaced points in the interval [0,10], i.e. sampling frequency 
            is 5 Hz.   
  paramu  = [u v N] defines discretization of maxima and minima ranges:  
            u is the lowest minimum considered, v the highest maximum and N  
            is the number of levels (u,v) included.  
  options = rind-options structure containing optional parameters 
            controlling the performance of the integration.  
            See rindoptset for details. 
    []    = default values are used.  
  
  SPEC2MMTPDF calculates densities of wave characteristics in a 
  stationary Gaussian transform process X(t) where  
   Y(t) = g(X(t)) (Y zero-mean Gaussian with spectrum given in input spec).  
   The tr. g can be estimated using lc2tr, dat2tr, hermitetr or ochitr. 
   
  Examples:% The joint density of zero separated Max2min cycles in time (a); 
           % in space (b); AcAt in time for nonlinear sea model (c):  
  
     Hm0=7;Tp=11; 
     S = jonswap(4*pi/Tp,[Hm0 Tp]);  
     Sk = spec2spec(S,'k1d');      
     L0 = spec2mom(S,1);  
     paramu = [sqrt(L0)*[-4 4] 41];  
     ft = spec2mmtpdf(S,0,'vmm',[],paramu); pdfplot(ft)         % a) 
     fs = spec2mmtpdf(Sk,0,'vmm');  figure, pdfplot(fs)         % b) 
     [sk, ku, me]=spec2skew(S);  
     g = hermitetr([],[sqrt(L0) sk ku me]); 
     Snorm=S; Snorm.S=S.S/L0; Snorm.tr=g; 
     ftg=spec2mmtpdf(Snorm,0,'AcAt',[],paramu); pdfplot(ftg)    % c) 
  
  See also  rindoptset, dat2tr, datastructures, wavedef, perioddef

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function f = spec2mmtpdf(spec,utc,def,paramt,paramu,options,bound)  
0002 %SPEC2MMTPDF Calculates joint density of Maximum, minimum and period. 
0003 %           
0004 %  CALL:  f   = spec2mmtpdf(S,u,def,paramt,paramu,options);  
0005 %  
0006 %    f    = pdf (density structure) of crests (trough) heights  
0007 %    S    = spectral density structure  
0008 %    u    = reference level (default the most frequently crossed level). 
0009 %   def   = string defining density returned 
0010 %            'Mm'   : maximum and the following minimum. (M,m) (default) 
0011 %           'rfc'   : maximum and the rainflow minimum height. 
0012 %           'AcAt'  : (crest,trough) heights.  
0013 %           'vMm'   : level v separated Maximum and minimum   (M,m)_v  
0014 %           'MmTMm' : maximum, minimum and period between (M,m,TMm) 
0015 %           'vMmTMm': level v separated Maximum, minimum and period 
0016 %                     between (M,m,TMm)_v 
0017 %           'MmTMd' : level v separated Maximum, minimum and the period 
0018 %                     from Max to level v-down-crossing (M,m,TMd)_v. 
0019 %           'MmTdm' : level v separated Maximum, minimum and the period from  
0020 %                     level v-down-crossing to min. (M,m,Tdm)_v 
0021 %           NB! All 'T' above can be replaced by 'L' to get wave length 
0022 %           instead.    
0023 % paramt  = [0 tn Nt] defines discretization of half period: tn is the 
0024 %           longest period considered while Nt is the number of points, 
0025 %           i.e. (Nt-1)/tn is the sampling frequnecy. paramt= [0 10 51] 
0026 %           implies that the halfperiods are considered at 51 linearly 
0027 %           spaced points in the interval [0,10], i.e. sampling frequency 
0028 %           is 5 Hz.   
0029 % paramu  = [u v N] defines discretization of maxima and minima ranges:  
0030 %           u is the lowest minimum considered, v the highest maximum and N  
0031 %           is the number of levels (u,v) included.  
0032 % options = rind-options structure containing optional parameters 
0033 %           controlling the performance of the integration.  
0034 %           See rindoptset for details. 
0035 %   []    = default values are used.  
0036 % 
0037 % SPEC2MMTPDF calculates densities of wave characteristics in a 
0038 % stationary Gaussian transform process X(t) where  
0039 %  Y(t) = g(X(t)) (Y zero-mean Gaussian with spectrum given in input spec).  
0040 %  The tr. g can be estimated using lc2tr, dat2tr, hermitetr or ochitr. 
0041 %  
0042 % Examples:% The joint density of zero separated Max2min cycles in time (a); 
0043 %          % in space (b); AcAt in time for nonlinear sea model (c):  
0044 % 
0045 %    Hm0=7;Tp=11; 
0046 %    S = jonswap(4*pi/Tp,[Hm0 Tp]);  
0047 %    Sk = spec2spec(S,'k1d');      
0048 %    L0 = spec2mom(S,1);  
0049 %    paramu = [sqrt(L0)*[-4 4] 41];  
0050 %    ft = spec2mmtpdf(S,0,'vmm',[],paramu); pdfplot(ft)         % a) 
0051 %    fs = spec2mmtpdf(Sk,0,'vmm');  figure, pdfplot(fs)         % b) 
0052 %    [sk, ku, me]=spec2skew(S);  
0053 %    g = hermitetr([],[sqrt(L0) sk ku me]); 
0054 %    Snorm=S; Snorm.S=S.S/L0; Snorm.tr=g; 
0055 %    ftg=spec2mmtpdf(Snorm,0,'AcAt',[],paramu); pdfplot(ftg)    % c) 
0056 % 
0057 % See also  rindoptset, dat2tr, datastructures, wavedef, perioddef 
0058   
0059  
0060 %       bound = 0 the distribution is approximated (default)  
0061 %             = 1 the upper and lower bounds for the distribution are computed. 
0062  
0063  
0064 %               -1 Integrate all by SADAPT for Ndim<9 and by KRBVRC otherwise  
0065 %               -2 Integrate all by SADAPT by Genz (1992) (Fast) 
0066 %               -3 Integrate all by KRBVRC by Genz (1993) (Fast) 
0067 %               -4 Integrate all by KROBOV by Genz (1992) (Fast) 
0068 %               -5 Integrate all by RCRUDE by Genz (1992) 
0069 %               -6 Integrate all by RLHCRUDE using MLHD and center of cell  
0070 %               -7 Integrate all by RLHCRUDE using  LHD and center of cell 
0071 %               -8 Integrate all by RLHCRUDE using MLHD and random point within the cell 
0072 %               -9 Integrate all by RLHCRUDE using LHD and random point within the cell 
0073  
0074 % Tested on : matlab 5.3  
0075 % History: by I. Rychlik 01.10.1998 with name minmax.m, new name: spec2cmat  
0076 % bounds by I.R. 02.01.2000.  
0077 % Revised by pab 09.05.2000  
0078 %  - added level u separated Max and min + period distributions  
0079 %  - new name: spec2cmat -> spec2mmtpdf  
0080 %  - also Made call with directional spectrum possible. 
0081 % revised by IR removing error in transformation 29 VI 2000 
0082 % revised by IR  changed default value for paramt and some ch. of 
0083 % help 2 VII 2000 
0084 % revised by IR, addopted for Matlab 6, 31 III 2001  
0085 % revised pab 10 April 2003 
0086 %  -added matlab version, i.e., a call to sp2mmtpdf   
0087 %revised pab 22Nov 2003 
0088 % removed  nit and speed from input and replaced it  
0089 %  with a rindoptions structure   
0090 %    nit  =  0,...,9. Dimension of numerical integration  (default 2). 
0091 %           -1,-2,... different important sampling type integrations. 
0092 %   speed = defines accuracy of calculations, by choosing different  
0093 %               parameters, possible values: 1,2...,9 (9 fastest, default 
0094 %               4). 
0095    
0096 defaultSpeed   = 4; 
0097 defaultoptions = rindoptset('speed',defaultSpeed,'nit',2,'method',0); 
0098 defaultoptions.speed = []; 
0099    
0100 if ((nargin==1) & (nargout <= 1) &  isequal(spec,'defaults')), 
0101   f = defaultoptions; 
0102   return 
0103 end  
0104 error(nargchk(1,7,nargin)) 
0105 startTime = clock;  
0106  
0107 ftype = freqtype(spec); 
0108 if nargin<3|isempty(def)  
0109   defnr=0;  
0110   if strcmp(ftype,'k') 
0111     def='L'; % distributions in space are default    
0112   else 
0113     def='T'; % distributions in time are default    
0114   end 
0115 else 
0116   switch lower(def)  
0117    case 'acat',    
0118     defnr =-2; % (Ac,At) 
0119    case 'rfc', 
0120     defnr =-1; % (M,m_rfc) 
0121    case 'mm', 
0122     defnr = 0; % max2min. (M,m) (default) 
0123    case {'mmtmm','mmlmm'}, 
0124     defnr = 1; % max2min and period inbetween (M,m,TMm) 
0125    case {'vmm'}, 
0126     defnr = 2; % level v separated Max2min   (Mv,mv) 
0127    case {'vmmtmm',  'vmmlmm'}, 
0128     defnr = 3; % level v separated Max2min and period inbetween (Mu,mu,TMm) 
0129    case {'mmtmd','vmmtmd','mmlmd','vmmlmd'}, 
0130     defnr = 4; % level v separated Max2min and period from Max to level v-down-crossing. 
0131    case {'mmtdm','vmmtdm','mmldm','vmmldm'},  
0132     defnr = 5; % level v separated Max2min and period from level v-down-crossing to min. 
0133     otherwise, error('Unknown def')  
0134   end  
0135   if defnr>=3|defnr==1, 
0136     def = upper(def(end-2)); % Store the kind of distribution that is 
0137                              % wanted: Period or wavelength 
0138   elseif strcmp(ftype,'k') 
0139     def='L'; % distributions in space are default    
0140   else 
0141     def='T'; % distributions in time are default    
0142   end 
0143 end                  
0144  
0145  
0146 switch upper(def(1)) 
0147   case {'L'},  spec = spec2spec(spec,'k1d') ;  ptxt='space';   
0148   case {'T'},  spec = spec2spec(spec,'freq');  ptxt='time';  
0149   otherwise,  error('Unknown def') 
0150 end 
0151  
0152 [S, xl4, L0, L2, L4, L1] = wnormspec(spec);  
0153 A    = sqrt(L0/L2);  
0154   
0155 if isfield(spec,'tr')  
0156    g=spec.tr;  
0157 else  
0158    g=[];  
0159 end  
0160 if isempty(g)  
0161   g  = [sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)']; 
0162 end 
0163  
0164 if nargin<6|isempty(options) 
0165   options = defaultoptions; 
0166 else 
0167   options = rindoptset(defaultoptions,options); 
0168 end 
0169  
0170 if nargin<7|isempty(bound)  
0171     bound=0;  
0172 end 
0173  
0174 if nargin<2|isempty(utc) 
0175   utc_d = gaus2dat([0 0],g); % most frequent crossed level  
0176   utc   = utc_d(1,2); 
0177 end 
0178  
0179 % transform reference level into Gaussian level 
0180 uu = dat2gaus([0. utc],g); 
0181 u  = uu(2); 
0182 disp(['The level u for Gaussian process = ' num2str(u)]) 
0183           
0184 if options.method>0  
0185    bound=0;  
0186 end  
0187 if nargin<4|isempty(paramt)  
0188   paramt = [0 5*pi*sqrt(L2/L4), 41];   
0189 end  
0190 if nargin<5|isempty(paramu)   
0191   paramu=[-5*sqrt(L0) 5*sqrt(L0) 41];  
0192 end  
0193  
0194 t0     = paramt(1);  
0195 tn     = paramt(2);  
0196 Ntime  = paramt(3); 
0197  
0198 t      = levels([0 tn/A Ntime]); % normalized times  
0199 Nstart = 1 + round(t0/tn*(Ntime-1)); % the starting point to evaluate  
0200  
0201  
0202 Nx = paramu(3); 
0203 if (defnr>1) 
0204   paramu(1) = max(0,paramu(1)); 
0205   if (paramu(2)<0),  
0206     error('Discretization levels must be larger than zero'),  
0207   end   
0208 end 
0209    
0210 %Transform amplitudes to Gaussian levels:     
0211 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
0212 h   = levels(paramu);  
0213 h   = reshape(h,Nx,1);  
0214 der = ones(Nx,1); 
0215  
0216 if defnr>1 % level v separated Max2min densities 
0217   der1 = der; 
0218   hg   = tranproc([utc+h der],g); 
0219   der  = abs(hg(:,2));  
0220   hg   = hg(:,1); % Gaussian levels above u 
0221   hg1  = tranproc([utc-h der1],g); 
0222   der1 = abs(hg1(:,2));  
0223   hg   = [hg;hg1(:,1)]; % Gaussian levels below u 
0224 else   % Max2min densities 
0225  
0226   hg  = tranproc([h der],g);  
0227   der = abs(hg(:,2));  
0228   der1= der; 
0229   hg  = hg(:,1); % Gaussian level  
0230 end 
0231  
0232 dt      = t(2)-t(1) ; 
0233  
0234 % Calculating covariances 
0235 %~~~~~~~~~~~~~~~~~~~~~~~~~ 
0236 nr = 4; % number of derivatives 
0237 R = spec2cov2(S,nr,Ntime-1,dt); 
0238   
0239 %NB!!! the spec2XXpdf.exe programmes is very sensitive to how you interpolate   
0240 %      the covariances, especially where the process is very dependent  
0241 %      and the covariance matrix is nearly singular. (i.e. for small t  
0242 %      and high levels of u if Tc and low levels of u if Tt)  
0243 %     The best is to interpolate the spectrum linearly so that S.S>=0 
0244 %     This makes sure that the covariance matrix is positive 
0245 %     semi-definitt, since the circulant spectrum are the eigenvalues of 
0246 %     the circulant covariance matrix. 
0247  
0248  
0249 callFortran = 0; %options.method<0;   
0250 if callFortran, % call fortran 
0251    ftmp = callSp2mmtpdfexe(R,dt,u,defnr,Nstart,hg,bound,options); 
0252    err = repmat(nan,size(ftmp)); 
0253 else 
0254   [ftmp,err,options] = cov2mmtpdf(R,dt,u,defnr,Nstart,hg,options);    
0255 end 
0256  
0257 f       = createpdf;  
0258 if isfield(S,'note') 
0259   f.note  = S.note; 
0260 end 
0261 f.options = options; 
0262 %f.u     = utc; 
0263  
0264 f.elapsedTime = etime(clock,startTime); 
0265  
0266 if Nx>2 
0267   f.labx{1}='Max [m]';  
0268   f.x{1}=h;  
0269   switch defnr  
0270    case -2, f.title = sprintf('Joint density of (Ac,At) in %s', ptxt);  
0271    case -1, f.title = sprintf('Joint density of (M,m_{rfc}) in %s', ptxt);  
0272    case  0, f.title = sprintf('Joint density of (M,m) in %s', ptxt);   
0273    case  1,  
0274     f.title = sprintf('Joint density of (M,m,%sMm) in %s', def(1), ptxt); 
0275    case  2,  
0276     f.title = sprintf('Joint density of (M,m)_{v=%2.5g} in %s',utc, ptxt);    
0277    case  3,  
0278     f.title = sprintf('Joint density of (M,m,%sMm)_{v=%2.5g} in %s',... 
0279               def(1),utc, ptxt);  
0280    case  4,  
0281     f.title = sprintf('Joint density of (M,m,%sMd)_{v=%2.5g} in %s',... 
0282               def(1),utc, ptxt);  
0283    case  5,  
0284     f.title = sprintf('Joint density of (M,m,%sdm)_{v=%2.5g} in %s',... 
0285               def(1),utc, ptxt);  
0286     otherwise, error('Unknown def')  
0287   end  
0288 else   
0289   f.note = [f.note 'Density is not scaled to unity']; 
0290   switch defnr  
0291    case {-2,-1,0,1},  
0292     f.title = sprintf('Density of (%sMm, M = %2.5g, m = %2.5g)',... 
0293               def(1),h(2),h(1));  
0294    case {2,3},  
0295     f.title = sprintf('Density of (%sMm, M = %2.5g, m = %2.5g)_{v=%2.5g}',... 
0296               def(1),h(2),-h(2),utc);  
0297    case 4,      
0298     f.title = sprintf('Density of (%sMd, %sMm, M = %2.5g, m = %2.5g)_{v=%2.5g}',... 
0299               def(1),def(1),h(2),-h(2),utc);  
0300    case 5,      
0301     f.title = sprintf('Density of (%sdm, %sMm, M = %2.5g, m = %2.5g)_{v=%2.5g}',... 
0302               def(1),def(1),h(2),-h(2),utc);  
0303    otherwise, error('Unknown def')  
0304   end 
0305 end 
0306  
0307  
0308  
0309 dh = h(2)-h(1); 
0310 if bound<0.5 
0311  if Nx>2  % amplitude distributions wanted 
0312    f.x{2}    = h; 
0313    f.labx{2} = 'min [m]'; 
0314    if defnr>1|defnr==-2 ,f.u  = utc;end   % save level u 
0315     
0316    if defnr>2 | defnr==1 
0317      ftmp = reshape(ftmp,Nx,Nx,Ntime); 
0318      err  = reshape(err,Nx,Nx,Ntime); 
0319      der0 = der1(:,ones(1,Nx)).*(der(:,ones(1,Nx)).'); 
0320      for ix =1:Ntime 
0321        ftmp(:,:,ix) = ftmp(:,:,ix).*der0;% dh*dh; 
0322        err(:,:,ix)  = err(:,:,ix).*der0; 
0323      end 
0324       
0325      ftmp   = ftmp/A; 
0326      err    = err/A; 
0327      f.x{3} = t(:)*A; 
0328      if strcmpi(def(1),'t') 
0329        f.labx{3} = 'period [sec]'; 
0330      else 
0331        f.labx{3} = 'wave length [m]'; 
0332      end 
0333    else 
0334      ftmp  = reshape(ftmp,Nx,Nx); 
0335      err   = reshape(err,Nx,Nx); 
0336      ftmp  = ftmp.*der(:,ones(1,Nx)).*(der(:,ones(1,Nx)).'); 
0337      err   = err.*der(:,ones(1,Nx)).*(der(:,ones(1,Nx)).'); 
0338      if (defnr==-1) 
0339        ftmp0 = fliplr(mctp2rfc(fliplr(ftmp))); 
0340        err  = abs(ftmp0-fliplr(mctp2rfc(fliplr(ftmp+err)))); 
0341        ftmp = ftmp0; 
0342      elseif (defnr==-2) 
0343         ftmp0=fliplr(mctp2tc(fliplr(ftmp),utc,paramu))*sqrt(L4*L0)/L2; 
0344     err =abs(ftmp0-fliplr(mctp2tc(fliplr(ftmp+err),utc,paramu))*sqrt(L4*L0)/L2); 
0345         index1=find(f.x{1}>0); 
0346         index2=find(f.x{2}<0); 
0347         ftmp=flipud(ftmp0(index2,index1)); 
0348     err =flipud(err(index2,index1)); 
0349         f.x{1} = f.x{1}(index1); 
0350         f.x{2} = abs(flipud(f.x{2}(index2))); 
0351      end 
0352    end 
0353    f.f = ftmp; 
0354    f.err = err; 
0355  else % Only time or wave length distributions wanted 
0356    f.f = ftmp/A; 
0357    f.err = err/A 
0358    f.x{1}=A*t'; 
0359    if strcmpi(def(1),'t') 
0360      f.labx{1} = 'period [sec]'; 
0361    else 
0362      f.labx{1} = 'wave length [m]'; 
0363    end 
0364    if defnr>3, 
0365      f.f   = reshape(f.f,[Ntime, Ntime]); 
0366      f.err = reshape(f.err,[Ntime, Ntime]); 
0367      f.x{2}= A*t'; 
0368      if strcmpi(def(1),'t') 
0369        f.labx{2} = 'period [sec]'; 
0370      else 
0371         f.labx{2} = 'wave length [m]'; 
0372      end 
0373    end 
0374       
0375  end 
0376   
0377 else    % bound >0.5 
0378   % The lines below must be fixed: (it is not correct yet 
0379   %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
0380   if Nx>2 
0381     ftmp_up=reshape(ftmp(:,1),Nx,Nx);  
0382     ftmp_lo=reshape(ftmp(:,2),Nx,Nx);  
0383     for i=1:Nx  
0384       ftmp_up(:,i)=dh*dh*der(i)*ftmp_up(:,i).*der;%*sqrt(-R(1,6)/R(1,4))/2/pi;  
0385       ftmp_lo(:,i)=dh*dh*der(i)*ftmp_lo(:,i).*der;%sqrt(-R(1,6)/R(1,4))/2/pi; 
0386     end  
0387     if (defnr==0) 
0388       f.f{1}=fliplr(mctp2rfc(fliplr(ftmp_up))); 
0389       f.f{2}=fliplr(mctp2rfc(fliplr(ftmp_lo))); 
0390     end 
0391     if (defnr==1) 
0392       f.f{1}=ftmp_up;  
0393       f.f{2}=ftmp_lo; 
0394     end 
0395     if (defnr==-1) 
0396       f.f{1}=fliplr(mctp2tc(fliplr(ftmp_up))); 
0397       f.f{2}=fliplr(mctp2tc(fliplr(ftmp_lo))); 
0398     end 
0399     f.x{2}=h;  
0400   else 
0401     size(ftmp) 
0402     f.f=ftmp/A;    
0403     f.x{1}=A*t.'; 
0404   end 
0405 end  
0406  
0407 % [f.cl,f.pl]=qlevels(f.f,[10 30 50 70 90 95 99 99.9],f.x{1},f.x{2}); 
0408  
0409 %pdfplot(f)  
0410  
0411 %Test of spec2mmtpdf 
0412 % cd  f:\matlab\matlab\wafo\source\sp2thpdfalan 
0413 % addpath f:\matlab\matlab\wafo ,initwafo, addpath f:\matlab\matlab\graphutil 
0414 % Hm0=7;Tp=11; S = jonswap(4*pi/Tp,[Hm0 Tp]); 
0415 % ft = spec2mmtpdf(S,0,'vMmTMm',[0.3,.4,11],[0 .00005 2]); 
0416  
0417 return % main 
0418  
0419 function dens  = callSp2mmtpdfexe(R,dt,u,defnr,Nstart,hg,bound,options) 
0420 % Write parameters to file 
0421 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
0422   Nx    = max(1,length(hg)); 
0423   if (defnr>1)  
0424     Nx = Nx/2; %level v separated max2min densities wanted      
0425   end 
0426   nr    = size(R,2)-1; 
0427   Ntime = size(R,1); 
0428    
0429    
0430   if exist('h.in'), delete h.in, end  
0431   fid = fopen('h.in','wt');  
0432   fprintf(fid,'%12.10f\n',hg);  
0433    fclose(fid);  
0434    for k=0:nr  
0435       filename=['Cd' int2str(k) '.in'];  
0436       if exist(filename),   
0437     delete(filename), 
0438       end  
0439       fid = fopen(filename,'wt');  
0440       fprintf(fid,'%12.10E \n',R(:,k+1));  
0441       fclose(fid);  
0442    end  
0443    %SCIS=0;  
0444    XSPLT = options.xsplit; 
0445    nit   = options.nit; 
0446    speed = options.speed; 
0447    seed  = options.seed; 
0448    SCIS  = abs(options.method); % method<=0 
0449     
0450   if exist('reflev.in'),  
0451     delete('reflev.in'), 
0452   end  
0453   disp('writing data')  
0454   fid=fopen('reflev.in','wt');  
0455   fprintf(fid,'%2.0f \n',Ntime);  
0456   fprintf(fid,'%2.0f \n',Nstart);  
0457   fprintf(fid,'%2.0f \n',nit);  
0458   fprintf(fid,'%2.0f \n',speed);  
0459   fprintf(fid,'%2.0f \n',SCIS);  
0460   fprintf(fid,'%2.0f \n',seed);  % select a random seed for rind   
0461   fprintf(fid,'%2.0f \n',Nx);  
0462   fprintf(fid,'%12.10E \n',dt); 
0463   fprintf(fid,'%12.10E \n',u); 
0464   fprintf(fid,'%2.0f \n',defnr); % def 
0465   fclose(fid);   
0466    
0467   disp('   Starting Fortran executable.')  
0468   if bound<0.5  
0469      dos([ wafoexepath 'sp2mmt70.exe']); %compiled sp2mmtpdf.f with rind70.f  
0470   else  
0471      dos([ wafoexepath 'sp2mM51.exe']); %compiled spec2tthpdf1.f with rind49.f  
0472   end  
0473   dens = load('dens.out'); 
0474   return 
0475    
0476    
0477   function [pdf,err, options] = cov2mmtpdf(R,dt,u,def,Nstart,hg,options) 
0478 %COV2MMTPDF Joint density of Maximum, minimum and period. 
0479 % 
0480 % CALL  [pdf, err, options] = cov2mmtpdf(R,dt,u,def,Nstart,hg,options) 
0481 % 
0482 % pdf     = calculated pdf size Nx x Ntime 
0483 % err     = error estimate    
0484 % options = requested and actual rindoptions used in integration. 
0485 % R       = [R0,R1,R2,R3,R4] column vectors with autocovariance and its 
0486 %          derivatives, i.e., Ri (i=1:4) are vectors with the 1'st to 4'th 
0487 %          derivatives of R0.  size Ntime x Nr+1 
0488 % dt      = time spacing between covariance samples, i.e.,  
0489 %           between R0(1),R0(2). 
0490 % u       = crossing level 
0491 % def     = integer defining pdf calculated: 
0492 %           0 : maximum and the following minimum. (M,m) (default) 
0493 %           1 : level v separated Maximum and minimum   (M,m)_v  
0494 %           2 : maximum, minimum and period between (M,m,TMm) 
0495 %           3 : level v separated Maximum, minimum and period 
0496 %               between (M,m,TMm)_v 
0497 %           4 : level v separated Maximum, minimum and the period 
0498 %               from Max to level v-down-crossing (M,m,TMd)_v. 
0499 %           5 : level v separated Maximum, minimum and the period from  
0500 %               level v-down-crossing to min. (M,m,Tdm)_v  
0501 % Nstart  = index to where to start calculation, i.e., t0 = t(Nstart) 
0502 % hg      = vector of amplitudes length Nx or 0 
0503 % options = rind options structure defining the integration parameters  
0504 %   
0505 % COV2MMTPDF computes joint density of the  maximum and the following     
0506 % minimum or level u separated maxima and minima + period/wavelength     
0507 % 
0508 % For DEF = 0,1 : (Maxima, Minima and period/wavelength) 
0509 %         = 2,3 : (Level v separated Maxima and Minima and 
0510 %                  period/wavelength between them) 
0511 %   
0512 % If Nx==1 then the conditional  density for  period/wavelength between 
0513 % Maxima and Minima given the Max and Min is returned 
0514 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
0515 % Y=  X'(t2)..X'(ts)..X'(tn-1)||X''(t1) X''(tn)|| X'(t1) X'(tn)  X(t1) X(tn)  
0516 % = [       Xt                   Xd                    Xc            ]   
0517 %                                                                        
0518 % Nt = tn-2, Nd = 2, Nc = 4 
0519 %                                                                        
0520 % Xt= contains Nt time points in the indicator function                  
0521 % Xd=    "     Nd    derivatives in Jacobian 
0522 % Xc=    "     Nc    variables to condition on                           
0523 %                                                                        
0524 % There are 3 (NI=4) regions with constant barriers:                     
0525 % (indI(1)=0);     for i\in (indI(1),indI(2)]    Y(i)<0.                 
0526 % (indI(2)=Nt)  ;  for i\in (indI(2)+1,indI(3)], Y(i)<0 (deriv. X''(t1))  
0527 % (indI(3)=Nt+1);  for i\in (indI(3)+1,indI(4)], Y(i)>0 (deriv. X''(tn))   
0528 %  
0529 % 
0530 % For DEF = 4,5 (Level v separated Maxima and Minima and 
0531 %                period/wavelength from Max to crossing) 
0532 % 
0533 % If Nx==1 then the conditional joint density for  period/wavelength 
0534 % between Maxima, Minima and Max to level v crossing given the Max and 
0535 % the min is returned 
0536 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
0537 % Y=  X'(t2)..X'(ts)..X'(tn-1)||X''(t1) X''(tn) X'(ts)|| X'(t1) X'(tn)  X(t1) X(tn) X(ts) 
0538 % = [       Xt                      Xd                     Xc            ] 
0539 %                                                                          
0540 % Nt = tn-2, Nd = 3, Nc = 5 
0541 %                                                                          
0542 % Xt= contains Nt time points in the indicator function                       
0543 % Xd=    "     Nd    derivatives                                              
0544 % Xc=    "     Nc    variables to condition on                                
0545 %                                                                             
0546 % There are 4 (NI=5) regions with constant barriers:                         
0547 % (indI(1)=0);     for i\in (indI(1),indI(2)]    Y(i)<0.                     
0548 % (indI(2)=Nt)  ;  for i\in (indI(2)+1,indI(3)], Y(i)<0 (deriv. X''(t1))     
0549 % (indI(3)=Nt+1);  for i\in (indI(3)+1,indI(4)], Y(i)>0 (deriv. X''(tn)) 
0550 % (indI(4)=Nt+2);  for i\in (indI(4)+1,indI(5)], Y(i)<0 (deriv. X'(ts))     
0551 %  
0552 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
0553 % 
0554    
0555 % History 
0556 % Revised pab 22Nov2003 
0557 %  -new inputarguments   
0558 %  -added err to output   
0559 %Revised pab 03.04.2003 
0560 % -translated from f90 to matlab 
0561 %Revised pab 22.04.2000 
0562 % - added mean separated min/max + (Tdm, TMd) period distributions 
0563 % - added scis  
0564  
0565 R0 = R(:,1); 
0566 R1 = R(:,2); 
0567 R2 = R(:,3); 
0568 R3 = R(:,4); 
0569 R4 = R(:,5); 
0570  
0571 Ntime = length(R0); 
0572  
0573 Nx0  = max(1,length(hg)); 
0574 Nx1  = Nx0; 
0575 %Nx0 = Nx1; % just plain Mm 
0576 if (def>1)  
0577   Nx1 = Nx0/2; 
0578  %  Nx0 = 2*Nx1;   % level v separated max2min densities wanted      
0579 end 
0580 %disp(sprintf('def = %d',def)) 
0581  
0582 % ***** The bound 'infinity' is set to 100*sigma ***** 
0583 XdInf = 100.d0*sqrt(R4(1)); 
0584 XtInf = 100.d0*sqrt(-R2(1)); 
0585  
0586 Nc = 4; 
0587 NI = 4; 
0588 Nd = 2;  
0589 Mb = 1;  
0590 Nj = 0; 
0591  
0592 Nstart = max(2,Nstart);  
0593  
0594 symmetry = 0; 
0595 isOdd = mod(Nx1,2); 
0596 if (def<=1) %THEN % just plain Mm  
0597    Nx = Nx1*(Nx1-1)/2; 
0598    IJ = (Nx1+isOdd)/2;          
0599    if (hg(1)+hg(Nx1)==0 &  (hg(IJ)==0 |hg(IJ)+hg(IJ+1)==0) )  
0600      symmetry=0; 
0601      disp(' Integration region symmetric') 
0602      % May save Nx1-isOdd integrations in each time step  
0603      % This is not implemented yet. 
0604      %Nx = Nx1*(Nx1-1)/2-Nx1+isOdd 
0605    end 
0606  
0607    % CC = normalizing constant = 1/ expected number of zero-up-crossings of X'  
0608    CC = 2*pi*sqrt(-R2(1)/R4(1));  
0609    %  XcScale = log(CC) 
0610    XcScale = log(2*pi*sqrt(-R2(1)/R4(1))); 
0611 else 
0612    % level u separated Mm 
0613    Nx = (Nx1-1)*(Nx1-1); 
0614    if ( abs(u)<=eps & hg(1)+hg(Nx1+1)==0 & (hg(Nx1)+hg(2*Nx1)==0) ) 
0615      symmetry=0; 
0616      disp(' Integration region symmetric') 
0617      % Not implemented for DEF <= 3 
0618      %IF (DEF.LE.3) Nx = (Nx1-1)*(Nx1-2)/2  
0619    end 
0620  
0621    if (def>3) %THEN 
0622       Nstart = max(Nstart,3); 
0623       Nc = 5; 
0624       NI = 5;  
0625       Nd = 3; 
0626    end %ENDIF 
0627    %CC= normalizing constant= 1/ expected number of u-up-crossings of X 
0628    CC = 2*pi*sqrt(-R0(1)/R2(1))*exp(0.5D0*u*u/R0(1));  
0629    XcScale = log(2*pi*sqrt(-R0(1)/R2(1))) + 0.5D0*u*u/R0(1); 
0630 end %ENDIF 
0631  
0632 options = setfield(options,'xcscale',XcScale); 
0633 opt0 = struct2cell(options); 
0634 %opt0 = opt0(1:10); 
0635 %seed = []; 
0636 %opt0 =  {SCIS,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS,seed,NIT1}; 
0637  
0638  
0639  
0640 if (Nx>1) 
0641    if ((def==0 | def==2)) %  (M,m) or (M,m)v distribution wanted 
0642       asize = [Nx1,Nx1]; 
0643    else            
0644      % (M,m,TMm), (M,m,TMm)v  (M,m,TMd)v or (M,M,Tdm)v distributions wanted  
0645       asize = [Nx1,Nx1,Ntime]; 
0646    end  
0647 elseif (def>3) % 
0648    % Conditional distribution for (TMd,TMm)v or (Tdm,TMm)v given (M,m)  wanted  
0649    asize = [1,Ntime,Ntime]; 
0650 else 
0651   % Conditional distribution for  (TMm) or (TMm)v given (M,m) wanted 
0652    asize = [1,1,Ntime]; 
0653 end 
0654     
0655 % Initialization 
0656 %~~~~~~~~~~~~~~~~~ 
0657 pdf = zeros(asize); 
0658 err = pdf; 
0659 BIG  = zeros(Ntime+Nc+1,Ntime+Nc+1); 
0660 ex   = zeros(1, Ntime + Nc + 1);           
0661  
0662 fxind = zeros(Nx,1); 
0663 xc    = zeros(Nc,Nx); 
0664     
0665 indI  = zeros(1,NI); 
0666 a_up  = zeros(1,NI-1); 
0667 a_lo  = zeros(1,NI-1);  
0668  
0669  
0670 % INFIN =  INTEGER, array of integration limits flags:  size 1 x Nb   (in) 
0671     %            if INFIN(I) < 0, Ith limits are (-infinity, infinity); 
0672     %            if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)]; 
0673     %            if INFIN(I) = 1, Ith limits are [Hlo(I), infinity); 
0674     %            if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)]. 
0675 INFIN = repmat(0,1,NI-1); 
0676 INFIN(3)  = 1; 
0677 a_up(1,3) = +XdInf; 
0678 a_lo(1,1:2) = [-XtInf -XdInf]; 
0679 if (def>3) 
0680    a_lo(1,4) = -XtInf;    
0681 end 
0682  
0683 IJ = 0 ; 
0684 if (def<=1) % THEN     % Max2min and period/wavelength 
0685    for I=2:Nx1 
0686       J = IJ+I-1; 
0687       xc(3,IJ+1:J) =  hg(I); 
0688       xc(4,IJ+1:J) =  hg(1:I-1).'; 
0689       IJ = J; 
0690    end %do 
0691 else 
0692    % Level u separated Max2min 
0693    xc(Nc,:) = u; 
0694    % Hg(1) = Hg(Nx1+1)= u => start do loop at I=2 since by definition we must have:  minimum<u-level<Maximum 
0695    for i=2:Nx1                 
0696       J = IJ+Nx1-1; 
0697       xc(3,IJ+1:J) =  hg(i);              % Max > u 
0698       xc(4,IJ+1:J) =  hg(Nx1+2:2*Nx1).';    % Min < u 
0699       IJ = J; 
0700    end %do    
0701 end %IF 
0702 if (def <=3) 
0703    h11 = fwaitbar(0,[],sprintf('Please wait ...(start at: %s)',datestr(now))); 
0704  
0705    for Ntd = Nstart:Ntime  
0706       %Ntd=tn 
0707       Ntdc = Ntd+Nc; 
0708       Nt   = Ntd-Nd; 
0709       indI(2) = Nt; 
0710       indI(3) = Nt+1; 
0711       indI(4) = Ntd; 
0712       % positive wave period  
0713       BIG(1:Ntdc,1:Ntdc) = covinput(BIG(1:Ntdc,1:Ntdc),R0,R1,R2,R3,R4,Ntd,0); 
0714       [fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),a_lo,a_up,indI, ... 
0715                 xc,Nt,opt0{:}); 
0716        
0717        
0718       %fxind  = CC*rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),xc,Nt,NIT1,... 
0719       %speed1,indI,a_lo,a_up); 
0720  
0721        
0722       if (Nx<2) %THEN 
0723          % Density of TMm given the Max and the Min. Note that the 
0724          % density is not scaled to unity 
0725          pdf(1,1,Ntd) = fxind(1); 
0726      err(1,1,Ntd) = err0(1).^2; 
0727          %GOTO 100 
0728       else 
0729        
0730          IJ = 0; 
0731          switch def 
0732          case {-2,-1,0},  % joint density of (Ac,At),(M,m_rfc) or (M,m). 
0733             for i = 2:Nx1  
0734                J = IJ+i-1; 
0735                pdf(1:i-1,i,1) = pdf(1:i-1,i,1)+fxind(IJ+1:J).'*dt;%*CC; 
0736            err(1:i-1,i,1) = err(1:i-1,i,1)+(err0(IJ+1:J).'*dt).^2; 
0737                IJ = J; 
0738             end %do  
0739          case  {1}  % joint density of (M,m,TMm) 
0740             for  i = 2:Nx1 
0741                J = IJ+i-1; 
0742                pdf(1:i-1,i,Ntd) = fxind(IJ+1:J).';%*CC 
0743            err(1:i-1,i,Ntd) = (err0(IJ+1:J).').^2;%*CC 
0744                IJ = J; 
0745             end %do 
0746          case {2},  % joint density of level v separated (M,m)v 
0747             for  i = 2:Nx1 
0748                J = IJ+Nx1-1; 
0749                pdf(2:Nx1,i,1) = pdf(2:Nx1,i,1)+fxind(IJ+1:J).'*dt;%*CC; 
0750            err(2:Nx1,i,1) = err(2:Nx1,i,1)+(err0(IJ+1:J).'*dt).^2; 
0751                IJ = J; 
0752             end %do  
0753          case {3}  
0754             % joint density of level v separated (M,m,TMm)v 
0755             for  i = 2:Nx1  
0756                J = IJ+Nx1-1; 
0757                pdf(2:Nx1,i,Ntd) = pdf(2:Nx1,i,Ntd)+fxind(IJ+1:J).';%*CC; 
0758            err(2:Nx1,i,Ntd) = err(2:Nx1,i,Ntd)+(err0(IJ+1:J).').^2; 
0759                IJ = J; 
0760             end %do  
0761          end % SELECT    
0762       end %ENDIF 
0763       waitTxt = sprintf('%s Ready: %d of %d',datestr(now),Ntd,Ntime); 
0764       fwaitbar(Ntd/Ntime,h11,waitTxt); 
0765        
0766    end %do 
0767    close(h11); 
0768    err = sqrt(err); 
0769    %  goto 800 
0770 else 
0771    %200 continue 
0772    waitTxt = sprintf('Please wait ...(start at: %s)',datestr(now)); 
0773    h11 = fwaitbar(0,[],waitTxt); 
0774    tnold= -1; 
0775    for tn = Nstart:Ntime, 
0776       Ntd  = tn+1; 
0777       Ntdc = Ntd + Nc; 
0778       Nt   = Ntd - Nd; 
0779       indI(2) = Nt; 
0780       indI(3) = Nt + 1; 
0781       indI(4) = Nt + 2; 
0782       indI(5) = Ntd; 
0783        
0784       if (~symmetry) %IF (SYMMETRY) GOTO 300 
0785          
0786     for ts = 2:tn-1, 
0787       % positive wave period 
0788       BIG(1:Ntdc,1:Ntdc) = covinput(BIG(1:Ntdc,1:Ntdc),R0,R1,R2,R3, ... 
0789                     R4,tn,ts,tnold);  
0790        
0791       [fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc), ... 
0792                   a_lo,a_up,indI,xc,Nt,opt0{:});  
0793        
0794       %tnold = tn; 
0795       switch (def), 
0796        case {3,4}  
0797         if (Nx==1) %THEN 
0798                % Joint density (TMd,TMm) given the Max and the min.  
0799                % Note the density is not scaled to unity 
0800                pdf(1,ts,tn) = fxind(1);%*CC 
0801                err(1,ts,tn) = err0(1).^2;%*CC 
0802         else 
0803           % 4,  gives level u separated Max2min and wave period 
0804           % from Max to the crossing of level u (M,m,TMd). 
0805           IJ = 0; 
0806           for  i = 2:Nx1, 
0807         J = IJ+Nx1-1;  
0808         pdf(2:Nx1,i,ts) = pdf(2:Nx1,i,ts)+ fxind(IJ+1:J).'*dt; 
0809         err(2:Nx1,i,ts) = err(2:Nx1,i,ts)+ (err0(IJ+1:J).'*dt).^2; 
0810         IJ = J; 
0811           end %do  
0812         end 
0813        case  {5} 
0814         if (Nx==1)  
0815           % Joint density (Tdm,TMm) given the Max and the min.  
0816           % Note the density is not scaled to unity 
0817           pdf(1,tn-ts+1,tn) = fxind(1); %*CC 
0818           err(1,tn-ts+1,tn) = err0(1).^2; 
0819         else   
0820           % 5,  gives level u separated Max2min and wave period from  
0821           % the crossing of level u to the min (M,m,Tdm). 
0822                    
0823           IJ = 0; 
0824           for  i = 2:Nx1  
0825         J = IJ+Nx1-1; 
0826         pdf(2:Nx1,i,tn-ts+1)=pdf(2:Nx1,i,tn-ts+1) + fxind(IJ+1:J).'*dt;%*CC; 
0827         err(2:Nx1,i,tn-ts+1)=err(2:Nx1,i,tn-ts+1) + (err0(IJ+1:J).'*dt).^2; 
0828         IJ = J; 
0829           end %do  
0830         end 
0831       end % SELECT 
0832     end%         enddo 
0833       else % exploit symmetry 
0834          %300   Symmetry 
0835          for ts = 2:floor(Ntd/2)      
0836        % Using the symmetry since U = 0 and the transformation is 
0837            % linear. 
0838        % positive wave period 
0839             BIG(1:Ntdc,1:Ntdc) = covinput(BIG(1:Ntdc,1:Ntdc),R0,R1,R2,R3, ... 
0840                       R4,tn,ts,tnold);  
0841          
0842             [fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex,a_lo,a_up,indI, ... 
0843                   xc,Nt,opt0{:});  
0844          
0845             %tnold = tn; 
0846             if (Nx==1) % THEN 
0847                % Joint density of (TMd,TMm),(Tdm,TMm) given the max and 
0848                % the min. 
0849                % Note that the density is not scaled to unity 
0850                pdf(1,ts,tn) = fxind(1);%*CC; 
0851            err(1,ts,tn) = err0(1).^2; 
0852                if (ts<tn-ts+1) %THEN 
0853                   pdf(1,tn-ts+1,tn) = fxind(1);%*CC; 
0854           err(1,tn-ts+1,tn) = err0(1).^2; 
0855                end 
0856                %GOTO 350 
0857             else 
0858           IJ = 0 ; 
0859           switch (def) 
0860                case {4}  
0861          
0862         % 4,  gives level u separated Max2min and wave period from  
0863         % Max to the crossing of level u (M,m,TMd). 
0864         for  i = 2:Nx1  
0865           J = IJ+Nx1-1; 
0866           pdf(2:Nx1,i,ts) = pdf(2:Nx1,i,ts) + fxind(IJ+1:J)*dt;%*CC; 
0867           err(2:Nx1,i,ts) = err(2:Nx1,i,ts) + (err0(IJ+1:J)*dt).^2; 
0868           if (ts.LT.tn-ts+1)  
0869             % exploiting the symmetry 
0870             pdf(i,2:Nx1,tn-ts+1) = pdf(i,2:Nx1,tn-ts+1)+fxind(IJ+1:J)*dt;%*CC  
0871             err(i,2:Nx1,tn-ts+1) = err(i,2:Nx1,tn-ts+1)+(err0(IJ+1:J)*dt).^2; 
0872           end 
0873           IJ = J; 
0874         end %do  
0875                case {5} 
0876                   % 5,   gives level u separated Max2min and wave period 
0877                   % from the crossing of level u to min (M,m,Tdm). 
0878                   for  i = 2:Nx1,  
0879                      J = IJ+Nx1-1; 
0880                       
0881                      pdf(2:Nx1,i,tn-ts+1)=pdf(2:Nx1,i,tn-ts+1)+fxind(IJ+1:J)*dt;  
0882              err(2:Nx1,i,tn-ts+1)=err(2:Nx1,i,tn-ts+1)+(err0(IJ+1:J)*dt).^2;  
0883                      if (ts.LT.tn-ts+1)  
0884                %*CC; % exploiting the symmetry 
0885                         pdf(i,2:Nx1,ts) = pdf(i,2:Nx1,ts)+ fxind(IJ+1:J)*dt; 
0886             err(i,2:Nx1,ts) = err(i,2:Nx1,ts)+ (err0(IJ+1:J)*dt).^2; 
0887                      end %ENDIF 
0888                      IJ = J; 
0889                   end %do  
0890                end %END SELECT  
0891             end 
0892             %350      
0893          end %do 
0894       end 
0895       waitTxt = sprintf('%s Ready: %d of %d',datestr(now),tn,Ntime); 
0896       fwaitbar(tn/Ntime,h11,waitTxt); 
0897        
0898       %400     print *,'Ready: ',tn,' of ',Ntime 
0899    end %do 
0900    close(h11); 
0901    err = sqrt(err); 
0902 end % if       
0903  
0904 %Nx1,size(pdf) def  Ntime 
0905 if (Nx>1)% THEN 
0906    IJ = 1; 
0907    if (def>2 | def ==1) 
0908       IJ = Ntime; 
0909    end 
0910    pdf = pdf(1:Nx1,1:Nx1,1:IJ); 
0911    err = err(1:Nx1,1:Nx1,1:IJ); 
0912 else 
0913    IJ = 1; 
0914    if (def>3) 
0915       IJ = Ntime; 
0916    end 
0917    pdf = squeeze(pdf(1,1:IJ,1:Ntime)); 
0918    err = squeeze(err(1,1:IJ,1:Ntime)); 
0919 end 
0920 return 
0921        
0922     
0923 function BIG = covinput(BIG,R0,R1,R2,R3,R4,tn,ts,tnold) 
0924 %COVINPUT Sets up the covariance matrix   
0925 % 
0926 % CALL BIG = covinput(BIG, R0,R1,R2,R3,R4,tn,ts) 
0927 %    
0928 %  BIG = covariance matrix for X = [Xt,Xd,Xc] in spec2mmtpdf problems. 
0929 %      
0930 % The order of the variables in the covariance matrix 
0931 % are organized as follows:  
0932 % for  ts <= 1: 
0933 %    X'(t2)..X'(ts),...,X'(tn-1) X''(t1),X''(tn)  X'(t1),X'(tn),X(t1),X(tn)  
0934 % = [          Xt               |      Xd       |          Xc             ] 
0935 % 
0936 % for ts > =2: 
0937 %    X'(t2)..X'(ts),...,X'(tn-1) X''(t1),X''(tn) X'(ts)  X'(t1),X'(tn),X(t1),X(tn) X(ts)  
0938 % = [          Xt               |      Xd               |          Xc             ] 
0939 % 
0940 % where  
0941 % 
0942 % Xt= time points in the indicator function 
0943 % Xd= derivatives 
0944 % Xc=variables to condition on 
0945   
0946 % Computations of all covariances follows simple rules: Cov(X(t),X(s)) = r(t,s), 
0947 % then  Cov(X'(t),X(s))=dr(t,s)/dt.  Now for stationary X(t) we have 
0948 % a function r(tau) such that Cov(X(t),X(s))=r(s-t) (or r(t-s) will give the same result). 
0949 % 
0950 % Consequently  Cov(X'(t),X(s))    = -r'(s-t)    = -sign(s-t)*r'(|s-t|) 
0951 %               Cov(X'(t),X'(s))   = -r''(s-t)   = -r''(|s-t|) 
0952 %               Cov(X''(t),X'(s))  =  r'''(s-t)  = sign(s-t)*r'''(|s-t|) 
0953 %               Cov(X''(t),X(s))   =  r''(s-t)   = r''(|s-t|) 
0954 % Cov(X''(t),X''(s)) =  r''''(s-t) = r''''(|s-t|) 
0955  
0956 if nargin<9|isempty(tnold) tnold = -1; end 
0957  
0958      
0959 if (ts>1) % THEN 
0960    shft = 1; 
0961    N    = tn + 5 + shft; 
0962    %Cov(Xt,Xc) 
0963    %for 
0964     i = 1:tn-2; 
0965     %j = abs(i+1-ts) 
0966     %BIG(i,N)  = -sign(R1(j+1),R1(j+1)*dble(ts-i-1)) %cov(X'(ti+1),X(ts))  
0967     j = i+1-ts; 
0968     tau = abs(j)+1; 
0969     %BIG(i,N)  = abs(R1(tau)).*sign(R1(tau).*j.'); 
0970     BIG(i,N)  = R1(tau).*sign(j.'); 
0971    %end 
0972    %Cov(Xc) 
0973    BIG(N         ,N) =  R0(1);       % cov(X(ts),X(ts)) 
0974    BIG(tn+shft+3 ,N) =  R0(ts);      % cov(X(t1),X(ts)) 
0975    BIG(tn+shft+4 ,N) =  R0(tn-ts+1); % cov(X(tn),X(ts)) 
0976    BIG(tn+shft+1 ,N) = -R1(ts);      % cov(X'(t1),X(ts)) 
0977    BIG(tn+shft+2 ,N) =  R1(tn-ts+1); % cov(X'(tn),X(ts)) 
0978    %Cov(Xd,Xc) 
0979    BIG(tn-1 ,N) =  R2(ts);      %cov(X''(t1),X(ts))      
0980    BIG(tn   ,N) =  R2(tn-ts+1); %cov(X''(tn),X(ts)) 
0981     
0982    %ADD a level u crossing  at ts 
0983     
0984    %Cov(Xt,Xd) 
0985    %for 
0986       i = 1:tn-2; 
0987       j = abs(i+1-ts); 
0988       BIG(i,tn+shft)  = -R2(j+1); %cov(X'(ti+1),X'(ts))  
0989    %end 
0990    %Cov(Xd)   
0991    BIG(tn+shft,tn+shft) = -R2(1);  %cov(X'(ts),X'(ts)) 
0992    BIG(tn-1   ,tn+shft) =  R3(ts); %cov(X''(t1),X'(ts)) 
0993    BIG(tn     ,tn+shft) = -R3(tn-ts+1);  %cov(X''(tn),X'(ts)) 
0994  
0995    %Cov(Xd,Xc) 
0996    BIG(tn+shft ,N       ) =  0.d0;        %cov(X'(ts),X(ts))  
0997    BIG(tn+shft,tn+shft+1) = -R2(ts);      % cov(X'(ts),X'(t1)) 
0998    BIG(tn+shft,tn+shft+2) = -R2(tn-ts+1); % cov(X'(ts),X'(tn)) 
0999    BIG(tn+shft,tn+shft+3) =  R1(ts);      % cov(X'(ts),X(t1)) 
1000    BIG(tn+shft,tn+shft+4) = -R1(tn-ts+1); % cov(X'(ts),X(tn)) 
1001          
1002   if (tnold == tn) 
1003      % A previous call to covinput with tn==tnold has been made 
1004      % need only to update  row and column N and tn+1 of big: 
1005      return 
1006      % make lower triangular part equal to upper and then return 
1007      for j=1:tn+shft 
1008         BIG(N,j)      = BIG(j,N) 
1009         BIG(tn+shft,j) = BIG(j,tn+shft)    
1010      end 
1011      for j=tn+shft+1:N-1 
1012         BIG(N,j) = BIG(j,N) 
1013         BIG(j,tn+shft) = BIG(tn+shft,j)    
1014      end 
1015      return 
1016   end %if 
1017   tnold = tn; 
1018 else 
1019    N = tn+4; 
1020    shft = 0; 
1021 end %if 
1022            
1023 if (tn>2)       
1024    %for 
1025    i=1:tn-2; 
1026    %cov(Xt) 
1027    %   for j=i:tn-2 
1028    %     BIG(i,j) = -R2(j-i+1)              % cov(X'(ti+1),X'(tj+1)) 
1029    %  end %do 
1030     
1031    BIG(i,i) = toeplitz(-R2(i));  % cov(Xt) =   % cov(X'(ti+1),X'(tj+1)) 
1032  
1033  
1034    %cov(Xt,Xc) 
1035    BIG(i      ,tn+shft+1) = -R2(i+1);         %cov(X'(ti+1),X'(t1))   
1036    BIG(tn-1-i ,tn+shft+2) = -R2(i+1);         %cov(X'(ti+1),X'(tn))   
1037    BIG(i      ,tn+shft+3) =  R1(i+1);         %cov(X'(ti+1),X(t1))   
1038    BIG(tn-1-i ,tn+shft+4) = -R1(i+1);         %cov(X'(ti+1),X(tn))   
1039    
1040   %Cov(Xt,Xd) 
1041   BIG(i,tn-1)       = R3(i+1);          %cov(X'(ti+1),X''(t1))   
1042   BIG(tn-1-i,tn)    =-R3(i+1);          %cov(X'(ti+1),X''(tn))  
1043   %end %do 
1044 end       
1045 %cov(Xd) 
1046 BIG(tn-1  ,tn-1  ) = R4(1); 
1047 BIG(tn-1  ,tn    ) = R4(tn);     %cov(X''(t1),X''(tn)) 
1048 BIG(tn    ,tn    ) = R4(1); 
1049  
1050 %cov(Xc) 
1051 BIG(tn+shft+3 ,tn+shft+3) = R0(1);        % cov(X(t1),X(t1)) 
1052 BIG(tn+shft+3 ,tn+shft+4) = R0(tn);       % cov(X(t1),X(tn)) 
1053 BIG(tn+shft+1 ,tn+shft+3) = 0.d0 ;        % cov(X(t1),X'(t1)) 
1054 BIG(tn+shft+2 ,tn+shft+3) = R1(tn);       % cov(X(t1),X'(tn)) 
1055 BIG(tn+shft+4 ,tn+shft+4) = R0(1) ;       % cov(X(tn),X(tn)) 
1056 BIG(tn+shft+1 ,tn+shft+4) =-R1(tn);       % cov(X(tn),X'(t1)) 
1057 BIG(tn+shft+2 ,tn+shft+4) = 0.d0;         % cov(X(tn),X'(tn))  
1058 BIG(tn+shft+1 ,tn+shft+1) =-R2(1);        % cov(X'(t1),X'(t1)) 
1059 BIG(tn+shft+1 ,tn+shft+2) =-R2(tn);       % cov(X'(t1),X'(tn)) 
1060 BIG(tn+shft+2 ,tn+shft+2) =-R2(1);        % cov(X'(tn),X'(tn)) 
1061 %Xc=X(t1),X(tn),X'(t1),X'(tn)  
1062 %Xd=X''(t1),X''(tn) 
1063 %cov(Xd,Xc) 
1064 BIG(tn-1  ,tn+shft+3) = R2(1);           %cov(X''(t1),X(t1))      
1065 BIG(tn-1  ,tn+shft+4) = R2(tn);          %cov(X''(t1),X(tn))      
1066 BIG(tn-1  ,tn+shft+1) = 0.d0  ;          %cov(X''(t1),X'(t1))      
1067 BIG(tn-1  ,tn+shft+2) = R3(tn);          %cov(X''(t1),X'(tn))      
1068 BIG(tn    ,tn+shft+3) = R2(tn);          %cov(X''(tn),X(t1))      
1069 BIG(tn    ,tn+shft+4) = R2(1);           %cov(X''(tn),X(tn))      
1070 BIG(tn    ,tn+shft+1) =-R3(tn);          %cov(X''(tn),X'(t1))      
1071 BIG(tn    ,tn+shft+2) = 0.d0;            %cov(X''(tn),X'(tn)) 
1072        
1073       
1074 % make lower triangular part equal to upper  
1075 %for j=1:N-1 
1076 %   for i=j+1:N 
1077 %      BIG(i,j) = BIG(j,i) 
1078 %   end %do 
1079 %end %do 
1080 lp      = find(tril(ones(size(BIG)),-1)); % indices to lower triangular part 
1081 BIGT    = BIG.'; 
1082 BIG(lp) = BIGT(lp); 
1083 return 
1084 % END  SUBROUTINE COV_INPUT 
1085          
1086  
1087  
1088  
1089  
1090  
1091  
1092  
1093  
1094  
1095

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003