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
PDF class constructor | |
Transforms x using the transformation g. | |
returns the frequency type of a Spectral density struct. | |
Fast display of wait bar. | |
Transforms xx using the inverse of transformation g. | |
Calculates discrete levels given the parameter matrix. | |
Rainflow matrix given a Markov matrix of a Markov chain of turning points | |
Frequencies of upcrossing troughs and crests using Markov chain of turning points. | |
Computes multivariate normal expectations | |
Create or alter RIND OPTIONS structure. | |
Computes covariance function and its derivatives, alternative version | |
Transforms between different types of spectra | |
Transforms process X and up to four derivatives | |
Returns the path to executables for the WAFO Toolbox | |
Normalize a spectral density such that m0=m2=1 | |
Current date and time as date vector. | |
Close figure. | |
String representation of date. | |
Delete file or graphics object. | |
Execute DOS command and return result. | |
Display message and abort function. | |
Elapsed time. | |
Check if variables or functions are defined. | |
Close file. | |
Open file. | |
Convert integer to string (Fast version). | |
True if arrays are numerically equal. | |
True if field is in structure array. | |
Convert string to lowercase. | |
Current date and time as date number. | |
Convert number to string. (Fast version) | |
Set structure field contents. | |
Write formatted data to string. | |
Remove singleton dimensions. | |
Compare strings. | |
Compare strings ignoring case. | |
Convert structure array to cell array. | |
Toeplitz matrix. | |
Extract lower triangular part. | |
Convert string to uppercase. |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
Intensity of trough-crest cycles computed from St | |
Intensity of rainflow cycles computed from St |
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
Comments or corrections to the WAFO group