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
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:
 createpdf PDF class constructor dat2gaus Transforms x using the transformation g. freqtype returns the frequency type of a Spectral density struct. fwaitbar Fast display of wait bar. gaus2dat Transforms xx using the inverse of transformation g. levels Calculates discrete levels given the parameter matrix. mctp2rfc Rainflow matrix given a Markov matrix of a Markov chain of turning points mctp2tc Frequencies of upcrossing troughs and crests using Markov chain of turning points. rind Computes multivariate normal expectations rindoptset Create or alter RIND OPTIONS structure. spec2cov2 Computes covariance function and its derivatives, alternative version spec2spec Transforms between different types of spectra tranproc Transforms process X and up to four derivatives wafoexepath Returns the path to executables for the WAFO Toolbox wnormspec Normalize a spectral density such that m0=m2=1 clock Current date and time as date vector. close Close figure. datestr String representation of date. delete Delete file or graphics object. dos Execute DOS command and return result. error Display message and abort function. etime Elapsed time. exist Check if variables or functions are defined. fclose Close file. fopen Open file. int2str Convert integer to string (Fast version). isequal True if arrays are numerically equal. isfield True if field is in structure array. lower Convert string to lowercase. now Current date and time as date number. num2str Convert number to string. (Fast version) setfield Set structure field contents. sprintf Write formatted data to string. squeeze Remove singleton dimensions. strcmp Compare strings. strcmpi Compare strings ignoring case. struct2cell Convert structure array to cell array. toeplitz Toeplitz matrix. tril Extract lower triangular part. upper Convert string to uppercase.
This function is called by:
 Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics wafofig10 Intensity of trough-crest cycles computed from St wafofig9 Intensity of rainflow cycles computed from St

## 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
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 %
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
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
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
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