SPEC2THPDF Joint density of amplitude and period/wave-length characteristics CALL: f = spec2thpdf(S,u,def,paramt,h,options,plotflag); f = density structure of wave characteristics of half-waves S = spectral density structure u = reference level (default the most frequently crossed level). def = 'Tc', gives crest period, Tc (default). 'Tt', gives trough period, Tt. 'TcAc', gives crest period and wave crest amplitude (Tc,Ac). 'TtAt', gives trough period and wave trough amplitude (Tt,At). 'TcfAc', gives crest front period and wave crest amplitude (Tcf,Ac). 'TtbAt', gives trough back period and wave trough amplitude (Ttb,At). 'TAc', gives minimum of crest front/back period and wave crest amplitude (min(Tcf,Tcb),Ac). 'TAt', gives minimum of trough front/back period and wave trough amplitude (min(Ttf,Ttb),At). NB! All 'T.' above can be replaced by 'L.' to get wave length instead. paramt = [t0 tn Nt] where t0, tn and Nt is the first value, last value and the number of points, respectively, for which the density will be computed. paramt= [5 5 51] implies that the density is computed only for T=5 and using 51 linearly spaced points in the interval [0,5]. h = vector of amplitudes (default levels([0 hmax 31])) note h >= 0 hmax evaluated from spec options = rind-options structure containing optional parameters controlling the performance of the integration. See rindoptset for details. plotflag = plots result if 1, no plot else (default 0) [] = default values are used. SPEC2THPDF calculates densities of wave characteristics of half waves 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 transformation, g, can be estimated or calculated using LC2TR, DAT2TR, HERMITETR or OCHITR. Note that due to symmetry the density of (Tcf,Ac) is equal to (Tcb,Ac). Similarly the density of (Ttf,At) is equal to (Ttb,At). Example: The density of (Tcf,Ac,Tc=5) is computed by opt = rindoptset('speed',5,'nit',3,'method',0); S = jonswap; f = spec2thpdf(S,[],'TcfAc',[5 5 51],[],opt,1); opt1 = rindoptset( 'speed',9,'nit',3,'method',1); paramt = [0 10,51] f = spec2thpdf(S,[],'Tc',paramt,[],opt1); pdfplot(f) See also rindoptset, dat2tr, datastructures, ampdef, perioddef
PDF class constructor | |
Transforms x using the transformation g. | |
Fast display of wait bar. | |
Transforms xx using the inverse of transformation g. | |
Initializes RIND options according to speed. | |
Calculates discrete levels given the parameter matrix. | |
Plot contents of pdf structures | |
Calculates quantile levels which encloses P% of PDF | |
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. | |
Error function. | |
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. | |
Compare strings ignoring case. | |
Convert structure array to cell array. | |
Toeplitz matrix. | |
Extract lower triangular part. | |
Convert string to uppercase. | |
Display warning message; disable or enable warning messages. |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
Probability density distributions (pdf) of wave period, Tt, | |
Joint distribution (pdf) of crest front velocity and wave height: | |
Joint distribution (pdf) of crest front period, Tcf, and crest amplitude, Ac | |
Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac | |
Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac for extremal waves |
0001 function f = spec2thpdf(spec,utc,def,paramt,h,options,plotflag) 0002 %SPEC2THPDF Joint density of amplitude and period/wave-length characteristics 0003 % 0004 % CALL: f = spec2thpdf(S,u,def,paramt,h,options,plotflag); 0005 % 0006 % f = density structure of wave characteristics of half-waves 0007 % S = spectral density structure 0008 % u = reference level (default the most frequently crossed level). 0009 % def = 'Tc', gives crest period, Tc (default). 0010 % 'Tt', gives trough period, Tt. 0011 % 'TcAc', gives crest period and wave crest amplitude (Tc,Ac). 0012 % 'TtAt', gives trough period and wave trough amplitude (Tt,At). 0013 % 'TcfAc', gives crest front period and wave crest amplitude (Tcf,Ac). 0014 % 'TtbAt', gives trough back period and wave trough amplitude (Ttb,At). 0015 % 'TAc', gives minimum of crest front/back period and wave crest 0016 % amplitude (min(Tcf,Tcb),Ac). 0017 % 'TAt', gives minimum of trough front/back period and wave trough 0018 % amplitude (min(Ttf,Ttb),At). 0019 % NB! All 'T.' above can be replaced by 'L.' to get wave 0020 % length instead. 0021 % paramt = [t0 tn Nt] where t0, tn and Nt is the first value, last value 0022 % and the number of points, respectively, for which 0023 % the density will be computed. paramt= [5 5 51] implies 0024 % that the density is computed only for T=5 and 0025 % using 51 linearly spaced points in the interval [0,5]. 0026 % h = vector of amplitudes (default levels([0 hmax 31])) note h >= 0 0027 % hmax evaluated from spec 0028 % options = rind-options structure containing optional parameters 0029 % controlling the performance of the integration. 0030 % See rindoptset for details. 0031 %plotflag = plots result if 1, no plot else (default 0) 0032 % [] = default values are used. 0033 % 0034 % SPEC2THPDF calculates densities of wave characteristics of half waves 0035 % in a stationary Gaussian transform process X(t) where 0036 % Y(t) = g(X(t)) (Y zero-mean Gaussian with spectrum given in input spec). 0037 % The transformation, g, can be estimated or calculated using LC2TR, 0038 % DAT2TR, HERMITETR or OCHITR. 0039 % 0040 % Note that due to symmetry the density of (Tcf,Ac) is equal to (Tcb,Ac). 0041 % Similarly the density of (Ttf,At) is equal to (Ttb,At). 0042 % 0043 % Example: The density of (Tcf,Ac,Tc=5) is computed by 0044 % 0045 % opt = rindoptset('speed',5,'nit',3,'method',0); 0046 % 0047 % S = jonswap; 0048 % f = spec2thpdf(S,[],'TcfAc',[5 5 51],[],opt,1); 0049 % 0050 % opt1 = rindoptset( 'speed',9,'nit',3,'method',1); 0051 % paramt = [0 10,51] 0052 % f = spec2thpdf(S,[],'Tc',paramt,[],opt1); 0053 % pdfplot(f) 0054 % 0055 % See also rindoptset, dat2tr, datastructures, ampdef, perioddef 0056 0057 % Tested on : matlab 5.3,6.X 0058 % History: by I. Rychlik 01.10.1998 with name wave_th1.m 0059 % revised by Per A. Brodtkorb 19.09.1999 0060 % revised by I.R. 30.09.1999, bugs removing. 0061 % revised by pab 21.10.1999 0062 % added option 'TAc', 'TAt' 0063 % updated to new pdf structure 0064 % revised by es 000322. Make call with directional spectrum possible. 0065 % Introduced 'L..' options for def. 0066 % Added plotflag and made some changes in the help 0067 % revised pab 28.03.2000 0068 % - added NIT=-3:-11 0069 % - added XSPLT, rateLHD to reflev.in 0070 % - fixed 2 bugs: 1) contour levels only for 2D pdf's 0071 % 2) title string missed some braces 0072 % revised pab 22.05.2000 0073 % - changed order of negative NIT's ie SCIS 0074 0075 % -1 Integrate all by SADAPT for Ndim<9 and by KRBVRC otherwise 0076 % -2 Integrate all by SADAPT by Genz (1992) (Fast) 0077 % -3 Integrate all by KRBVRC by Genz (1993) (Fast) 0078 % -4 Integrate all by KROBOV by Genz (1992) (Fast) 0079 % -5 Integrate all by RCRUDE by Genz (1992) 0080 % -6 Integrate all by RLHCRUDE using MLHD and center of cell 0081 % -7 Integrate all by RLHCRUDE using LHD and center of cell 0082 % -8 Integrate all by RLHCRUDE using MLHD and random point within the cell 0083 % -9 Integrate all by RLHCRUDE using LHD and random point within the cell 0084 % revised by IR removing error in transformation 29 VI 2000 0085 % revised by IR removing error in normalization when t0=tn defnr=3,-3 1 VII 2000 0086 % revised by IR adapted to MATLAB6 + reshaping h, 8 II 2001 0087 % revised by jr 01.03.28 changed 'h.in' to fid at line 181 0088 % revised pab 03.03.03 0089 % - added sp2thpdf.m 0090 % revised pab 19.11.2003 0091 % -added err to output 0092 % -removed some code and replaced it with a call to spec2cov2 0093 % -removed nit and speed from input/output and replaced it 0094 % with a rindoptions structure 0095 % nit = 0,...,9. Dimension of numerical integration (default 2). 0096 % -1,-2,... different important sampling type integrations. 0097 % speed = defines accuracy of calculations, by choosing different 0098 % parameters, possible values: 1,2...,9 (9 fastest, default 4). 0099 0100 0101 defaultSpeed = 7; 0102 defaultoptions = initoptions(defaultSpeed); 0103 defaultoptions.nit = 2; 0104 defaultoptions.speed = []; 0105 0106 if ((nargin==1) & (nargout <= 1) & isequal(spec,'defaults')), 0107 f = defaultoptions; 0108 return 0109 end 0110 error(nargchk(1,7,nargin)) 0111 startTime = clock; 0112 if nargin<3|isempty(def) 0113 def='tc'; 0114 end 0115 if strcmpi('l',def(1)) 0116 spec = spec2spec(spec,'k1d'); 0117 elseif strcmpi('t',def(1)) 0118 spec = spec2spec(spec,'freq'); 0119 else 0120 error('Unknown def') 0121 end 0122 switch lower(def(2:end)) 0123 case 'c', defnr = 1; 0124 case 't', defnr =-1; 0125 case 'cac', defnr = 2; 0126 case 'tat', defnr =-2; 0127 case {'cfac','cbac'} , defnr = 3; 0128 case {'tbat', 'tfat'}, defnr =-3; 0129 case {'ac'} , defnr = 4; 0130 case {'at'}, defnr =-4; 0131 otherwise, error('Unknown def') 0132 end 0133 0134 0135 [S, xl4, L0, L2, L4, L1] = wnormspec(spec); 0136 A = sqrt(L0/L2); 0137 0138 if isfield(spec,'tr') 0139 g = spec.tr; 0140 else 0141 g = []; 0142 end 0143 if isempty(g) 0144 g = [sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)']; 0145 end 0146 if nargin<2|isempty(utc) 0147 utc_d = gaus2dat([0, 0],g); % most frequently crossed level 0148 utc = utc_d(1,2); 0149 end 0150 0151 % transform reference level into Gaussian level 0152 uu = dat2gaus([0., utc],g); 0153 u = uu(2); 0154 disp(['The level u for Gaussian process = ', num2str(u)]) 0155 0156 if nargin<6|isempty(options) 0157 options = defaultoptions; 0158 else 0159 options = rindoptset(defaultoptions,options); 0160 end 0161 0162 if nargin<7|isempty(plotflag) 0163 plotflag=0; 0164 end 0165 0166 if nargin<4|isempty(paramt) 0167 z2 = u^2/2; 0168 z = -sign(defnr)*u/sqrt(2); 0169 expectedMaxPeriod = 1.5*ceil(2*pi*A*exp(z)*(0.5+erf(z)/2)) 0170 paramt = [0, expectedMaxPeriod, 51]; 0171 end 0172 0173 t0 = paramt(1); 0174 tn = paramt(2); 0175 Ntime = paramt(3); 0176 t = levels([0, tn/A, Ntime]); % normalized times 0177 Nstart = 1 + round(t0/tn*(Ntime-1)); % index to starting point to evaluate 0178 0179 if abs(defnr)<2 0180 nr = 2; 0181 Nx = 1; 0182 hg = []; 0183 else 0184 nr = 4; 0185 if nargin<5|isempty(h) 0186 Nx = 31; 0187 px = gaus2dat([0., u;1, 4.*sign(defnr)],g); 0188 px = abs(px(2,2)-px(1,2)); 0189 h = levels([0, px, Nx]).'; 0190 else 0191 h = sort(abs(h(:))); % make sure values are positive and increasing 0192 Nx = length(h); 0193 h = reshape(h,Nx,1); 0194 end 0195 %Transform amplitudes to Gaussian levels: 0196 der = ones(Nx,1); % dh/dh=1 0197 hg = tranproc([utc+sign(defnr)*h, der],g); 0198 der = abs(hg(:,2)); 0199 hg = hg(:,1); % Gaussian level 0200 0201 end 0202 0203 dt = t(2)-t(1); 0204 % Calculating covariances 0205 %~~~~~~~~~~~~~~~~~~~~~~~~ 0206 R = spec2cov2(S,nr,Ntime-1,dt); 0207 0208 callFortran = 0; %options.method<0; 0209 if callFortran, % call fortran program 0210 dens = callSp2thpdfexe(R,dt,u,defnr,Nstart,hg,options); 0211 err = repmat(nan,size(dens)); 0212 else 0213 [dens,err,options] = cov2thpdf(R,dt,u,defnr,Nstart,hg,options); 0214 end 0215 0216 def1 = upper(def(1)); 0217 0218 if abs(defnr)>1 0219 f = createpdf(2); 0220 f.err = reshape(err/A,Nx,Ntime); 0221 f.f = reshape(dens/A,Nx,Ntime).*der(:,ones(1,Ntime)); 0222 f.x{2} = h; 0223 f.labx{2}='amplitude [m]'; 0224 else 0225 f = createpdf(1); 0226 f.f = dens(:)/A; 0227 f.err = err(:)/A; 0228 end 0229 if isfield(spec,'note') 0230 f.note = spec.note; 0231 end 0232 dt = dt*A; 0233 0234 switch defnr 0235 case 1, Htxt = sprintf('Density of %sc',def1); 0236 case -1, Htxt = sprintf('Density of %st',def1); 0237 case 2, Htxt = sprintf('Joint density of (%sc,Ac)',def1); 0238 case -2, Htxt = sprintf('Joint density of (%st,At)',def1); 0239 case 3, 0240 if (t0==tn) 0241 Htxt = ['Joint density of (',def1,'cf,Ac,',def1,'c=' num2str(tn) ')']; 0242 f.f = f.f/dt; 0243 f.err = f.err/dt; 0244 else 0245 Htxt = ['Joint density of (',def1,'cf,Ac)']; 0246 end 0247 case -3, 0248 if t0==tn 0249 Htxt = ['Joint density of (',def1,'cb,At,',def1,'t=' num2str(tn) ')']; 0250 f.f = f.f/dt; 0251 f.err = f.err/dt; 0252 else 0253 Htxt=['Joint density of (',def1,'cb,At)']; 0254 end 0255 case 4, Htxt = ['Joint density of (min(',def1,'cf,',def1,'cb),Ac)']; 0256 case -4, Htxt = ['Joint density of (min(',def1,'tf.',def1,'tb),At)']; 0257 end 0258 0259 Htxt = [Htxt,sprintf('_{v =%2.5g}',utc)]; 0260 0261 f.title = Htxt; 0262 if strcmpi('t',def(1)) 0263 f.labx{1} = 'period [s]'; 0264 else 0265 f.labx{1} = 'wave length [m]'; 0266 end 0267 0268 f.x{1} = t.'*A; 0269 f.options = options; 0270 f.u = utc; 0271 f.elapsedTime = etime(clock,startTime); 0272 0273 if abs(defnr)>1, 0274 try 0275 [f.cl,f.pl] = qlevels(f.f,[10, 30, 50, 70, 90, 95, 99],f.x{1},f.x{2}); 0276 catch 0277 warning('Unable to calculate contourlevels') 0278 end 0279 end 0280 0281 if plotflag 0282 pdfplot(f) 0283 end 0284 return % Main 0285 0286 function dens = callSp2thpdfexe(R,dt,u,defnr,Nstart,hg,options) 0287 0288 Nx = max(1,length(hg)); 0289 nr = size(R,2)-1; 0290 Ntime = size(R,1); 0291 0292 for k=0:nr 0293 filename=['Cd', int2str(k), '.in']; 0294 0295 if exist(filename) 0296 delete(filename) 0297 end 0298 fid=fopen(filename,'wt'); 0299 fprintf(fid,'%12.10E \n',R(:,k+1)); 0300 fclose(fid); 0301 end 0302 0303 %disp(SCIS) 0304 if ~isempty(hg) 0305 if exist('h.in'), delete h.in, end 0306 fid=fopen('h.in','wt'); 0307 fprintf(fid,'%12.10f\n',hg); 0308 fclose(fid) 0309 end 0310 0311 rateLHD=3; 0312 XSPLT = options.xsplit; 0313 nit = options.nit; 0314 speed = options.speed; 0315 seed = options.seed; 0316 SCIS = abs(options.method); % method<=0 0317 0318 if exist('reflev.in'), delete reflev.in, end 0319 disp('writing data') 0320 fid=fopen('reflev.in','wt'); 0321 fprintf(fid,'%12.10E \n',u); 0322 fprintf(fid,'%2.0f \n',defnr); 0323 fprintf(fid,'%2.0f \n',Ntime); 0324 fprintf(fid,'%2.0f \n',Nstart); 0325 fprintf(fid,'%2.0f \n',nit); 0326 fprintf(fid,'%2.0f \n',speed); 0327 fprintf(fid,'%2.0f \n',SCIS); 0328 fprintf(fid,'%2.0f \n',seed); % select a random seed for rind 0329 fprintf(fid,'%2.0f \n',Nx); 0330 fprintf(fid,'%12.10E \n',dt); 0331 fprintf(fid,'%2.0f \n',rateLHD); 0332 fprintf(fid,'%12.10E \n',XSPLT); 0333 fclose(fid); 0334 disp(' Starting Fortran executable.') 0335 %dos([ wafoexepath 'sp2thpdf70prof.exe']) 0336 dos([ wafoexepath 'sp2thpdf70.exe']); 0337 %dos([ wafoexepath 'sp2thpdfalan.exe']); % using EXINV 0338 %dos([ wafoexepath 'sp2thpdfalan3.exe']); % using EXINV 0339 %dos([ wafoexepath 'sp2thpdfalanEX.exe']); % using EXINV 0340 %dos([ wafoexepath 'sp2thpdfalanFI.exe']); % using FIINV 0341 0342 %dos([ wafoexepath 'sp2thpdfalanOld.exe']); 0343 0344 dens = load('dens.out'); 0345 return 0346 0347 function [pdf,err,options] = cov2thpdf(R,dt,u,def,Nstart,h,options) 0348 %COV2THPDF Joint density of amplitude and period/wave-length 0349 % 0350 % CALL [pdf, err, options] = cov2thpdf(R,dt,u,def,Nstart,h,options) 0351 % 0352 % pdf = calculated pdf size Nx x Ntime 0353 % err = error estimate 0354 % options = requested and actual rindoptions used in integration. 0355 % R = [R0,R1,R2,R3,R4] column vectors with autocovariance and its 0356 % derivatives, i.e., Ri (i=1:4) are vectors with the 1'st to 4'th 0357 % derivatives of R0. size Ntime x Nr+1 0358 % dt = time spacing between covariance samples, i.e., 0359 % between R0(1),R0(2). 0360 % u = crossing level 0361 % def = integer defining pdf calculated: 0362 % 1, gives half wave period, Tc (default). 0363 % -1, gives half wave period, Tt. 0364 % 2, gives half wave period and wave crest amplitude (Tc,Ac). 0365 % -2, gives half wave period and wave trough amplitude (Tt,At). 0366 % 3, gives crest front period and wave crest amplitude (Tcf,Ac). 0367 % -3, gives trough back period and wave trough amplitude (Ttb,At). 0368 % 4, gives minimum of crest front/back period and wave crest 0369 % amplitude (min(Tcf,Tcb),Ac). 0370 % -4, gives minimum of trough front/back period and wave trough 0371 % amplitude (min(Ttf,Ttb),At). 0372 % Nstart = index to where to start calculation, i.e., t0 = t(Nstart) 0373 % h = vector of amplitudes length Nx or 0 0374 % options = rind options structure defining the integration parameters 0375 % 0376 % COV2THPDF program computes density of S_i,Hi,T_i in a gaussian process, 0377 % i.e., quart wavelength (up-crossing to crest) and crest amplitude. 0378 % 0379 % See also rind, rindoptset 0380 0381 %History: 0382 % Revised pab 22Nov2003 0383 % -new inputarguments 0384 % revised Per A. Brodtkorb 03.03.2003 0385 % -converted from fortran 90 to matlab 0386 % -introduced XcScale instead of CC 0387 % revised Per A. Brodtkorb 04.04.2000 0388 % - 0389 % revised Per A. Brodtkorb 23.11.99 0390 % - fixed a bug in calculating pdf for def = +/- 4 0391 % revised Per A. Brodtkorb 03.11.99 0392 % - added def = +/-4 0393 % revised Per A. Brodtkorb 23.09.99 0394 % - minor changes to covinput 0395 % - removed the calculation of the transformation to spec2thpdf.m 0396 % by Igor Rychlik 0397 R0 = R(:,1); 0398 R1 = R(:,2); 0399 R2 = R(:,3); 0400 nr = size(R,2)-1; 0401 if nr <3 0402 R3 = []; 0403 R4 = []; 0404 else 0405 R3 = R(:,4); 0406 R4 = R(:,5); 0407 end 0408 0409 Ntime = length(R0); 0410 Nx = max(1,length(h)); 0411 0412 %C ***** The bound 'infinity' is set to 100*sigma ***** 0413 XdInf = 100.d0*sqrt(-R2(1)); 0414 XtInf = 100.d0*sqrt(R0(1)); 0415 % normalizing constant 0416 %CC=TWPI*SQRT(-R0(1)/R2(1))*exp(u*u/(2.d0*R0(1)) ); 0417 XcScale = log(2*pi*sqrt(-R0(1)/R2(1))) + u*u/(2.d0*R0(1)); 0418 0419 options = setfield(options,'xcscale',XcScale); 0420 0421 pdf = zeros(Nx,Ntime); 0422 err = zeros(Nx,Ntime); 0423 0424 fxind = zeros(Nx,1); 0425 err0 = fxind; 0426 0427 ABSEPS = options.abseps; 0428 0429 opt0 = struct2cell(options); 0430 %opt0 = opt0(1:10); 0431 %opt0 = {SCIS,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS}; 0432 if (abs(def)<2) 0433 NI = 4; 0434 Nd = 2; 0435 Nc = 2; 0436 Mb = 1; 0437 Nx = 1; 0438 BIG = zeros(Ntime+Nc,Ntime+Nc); 0439 xc = zeros(Nc,Nx); 0440 ex = zeros(1,Ntime+Nc); 0441 indI = zeros(1,NI); 0442 INFIN = repmat(2,1,NI-1); 0443 a_up = zeros(Mb,NI-1); 0444 a_lo = zeros(Mb,NI-1); 0445 0446 0447 xc(1,1)=u; 0448 xc(2,1)=u; 0449 % INFIN = INTEGER, array of integration limits flags: size 1 x Nb (in) 0450 % if INFIN(I) < 0, Ith limits are (-infinity, infinity); 0451 % if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)]; 0452 % if INFIN(I) = 1, Ith limits are [Hlo(I), infinity); 0453 % if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)]. 0454 0455 if (def>0) %then 0456 INFIN(1:2) = 1; 0457 INFIN(3) = 0; 0458 a_up(1,1) = u+XtInf ; 0459 a_lo(1,1)= u; 0460 a_up(1,2)= XdInf; 0461 a_lo(1,3)=-XdInf; 0462 else 0463 INFIN(1:2) = 0; 0464 INFIN(3) = 1; 0465 a_up(1,1)=u ; 0466 a_lo(1,1)=u-XtInf; 0467 a_lo(1,2)=-XdInf; 0468 a_up(1,3)= XdInf ; 0469 end %if 0470 %print *,'Nstart',Nstart 0471 Nstart=max(2,Nstart) ; 0472 %print *,'Nstart',Nstart 0473 waitTxt = sprintf('Please wait ...(start at: %s)',datestr(now)); 0474 h11 = fwaitbar(0,[],waitTxt); 0475 0476 for Ntd=Nstart:Ntime, 0477 %CALL COV_INPUT2(BIG,Ntd, R0,R1,R2) 0478 BIG = covinput(BIG,R0,R1,R2,R3,R4,Ntd,-1); % positive wave period 0479 Nt = Ntd-Nd; 0480 indI(2) = Nt; 0481 indI(3) = Nt+1; 0482 indI(4) = Ntd; 0483 Ntdc = Ntd+Nc; 0484 [pdf(Ntd),err(Ntd)] = rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),... 0485 a_lo,a_up,indI,xc,Nt,opt0{:}); 0486 waitTxt = sprintf('%s Ready: %d of %d',datestr(now),Ntd,Ntime); 0487 fwaitbar(Ntd/Ntime,h11,waitTxt); 0488 0489 %disp('sp2thpdf hit a button to continue') 0490 %pause 0491 end 0492 close(h11) 0493 0494 else 0495 XddInf = 100.d0*sqrt(R4(1)); 0496 NI=7; 0497 Nd=3; 0498 Nc=4; 0499 Mb=2; 0500 BIG = zeros(Ntime+Nc+1,Ntime+Nc+1); 0501 xc = zeros(Nc,Nx); 0502 ex = zeros(1,Ntime+Nc+1); 0503 indI = zeros(1,NI); 0504 INFIN = repmat(2,1,NI-1); 0505 a_up = zeros(Mb,NI-1); 0506 a_lo = zeros(Mb,NI-1); 0507 0508 xc(1,1:Nx) = h(1:Nx).'; 0509 xc(2,1:Nx) = u; 0510 xc(3,1:Nx) = u; 0511 xc(4,1:Nx) = 0.d0; 0512 0513 % INFIN = INTEGER, array of integration limits flags: size 1 x Nb (in) 0514 % if INFIN(I) < 0, Ith limits are (-infinity, infinity); 0515 % if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)]; 0516 % if INFIN(I) = 1, Ith limits are [Hlo(I), infinity); 0517 % if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)]. 0518 if (def>0) % then 0519 INFIN(2) = -1; 0520 INFIN(4) = 0; 0521 INFIN(5) = 1; 0522 INFIN(6) = 0; 0523 a_up(2,1) = 1.d0; %*h 0524 a_lo(1,1) = u; 0525 a_up(1,2) = XtInf; % X(ts) is redundant 0526 a_lo(1,2) = -XtInf; 0527 a_up(2,2) = 1.d0; % *h 0528 a_lo(2,2) = 1.d0; % *h 0529 a_up(2,3) = 1.d0; %*h 0530 a_lo(1,3) = u; 0531 0532 a_lo(1,4) = -XddInf; 0533 a_up(1,5) = XdInf; 0534 a_lo(1,6) = -XdInf ; 0535 else %def<0 0536 INFIN(2) = -1; 0537 INFIN(4) = 1; 0538 INFIN(5) = 0; 0539 INFIN(6) = 1; 0540 a_up(1,1) = u ; 0541 a_lo(2,1) = 1.d0; %*h 0542 a_up(1,2) = XtInf; % X(ts) is redundant 0543 a_lo(1,2) = -XtInf; 0544 a_up(2,2) = 1.d0; % *h 0545 a_lo(2,2) = 1.d0; % *h 0546 a_up(1,3) = u; 0547 a_lo(2,3) = 1.d0; %*h 0548 a_up(1,4) = XddInf; 0549 a_lo(1,5) = -XdInf; 0550 a_up(1,6) = XdInf; 0551 end %if 0552 EPSOLD = ABSEPS; 0553 Nstart = max(Nstart,3); 0554 0555 waitTxt = sprintf('Please wait ...(start at: %s)',datestr(now)); 0556 h11 = fwaitbar(0,[],waitTxt); 0557 0558 for tn = Nstart:Ntime, 0559 Ntd = tn+1; 0560 Nt = Ntd-Nd; 0561 Ntdc = Ntd+Nc; 0562 indI(4) = Nt; 0563 indI(5) = Nt+1; 0564 indI(6) = Nt+2; 0565 indI(7) = Ntd; 0566 0567 %opt0{3} =min(sqrt(tn)*EPSOLD*0.5D0,0.1D0); 0568 0569 for ts=2:floor((tn+1)/2.d0), 0570 %print *,'ts,tn' ,ts,tn 0571 BIG(1:Ntdc,1:Ntdc) = covinput(BIG(1:Ntdc,1:Ntdc),R0,R1,R2,R3,R4,tn,ts); % positive wave period 0572 indI(2) = ts-2; 0573 indI(3) = ts-1; 0574 [fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),a_lo, ... 0575 a_up,indI,xc,Nt,opt0{:}); 0576 0577 %disp('sp2thpdf hit a button to continue') 0578 %pause 0579 0580 switch abs(def) 0581 case {2} 0582 % 2, gives half wave period and wave crest amplitude (Tc,Ac). 0583 % -2, gives half wave period and wave trough amplitude (Tt,At). 0584 if (ts == tn-ts+1) % then 0585 pdf(1:Nx,tn) = pdf(1:Nx,tn) + fxind.'*dt; 0586 err(1:Nx,tn) = err(1:Nx,tn) + (err0(:)*dt).^2; 0587 else 0588 pdf(1:Nx,tn) = pdf(1:Nx,tn) + fxind.'*2.d0*dt; 0589 err(1:Nx,tn) = err(1:Nx,tn) + (err0(:)*2.0*dt).^2; 0590 end %if 0591 case {3} 0592 % 3, gives crest front period and wave crest amplitude (Tcf,Ac). 0593 % -3, gives trough back period and wave trough amplitude (Ttb,At). 0594 pdf(1:Nx,ts) = pdf(1:Nx,ts) + fxind.'*dt; 0595 err(1:Nx,ts) = err(1:Nx,ts) + (err0(:)*dt).^2; 0596 if ((ts<tn-ts+1)) 0597 % exploiting the symmetry 0598 pdf(1:Nx,tn-ts+1) = pdf(1:Nx,tn-ts+1) + fxind.'*dt; 0599 err(1:Nx,tn-ts+1) = err(1:Nx,tn-ts+1) + (err0(:)*dt).^2; 0600 end %if 0601 case {4} 0602 % 4, gives minimum of crest front/back period and wave crest 0603 % amplitude (min(Tcf,Tcb),Ac). 0604 % -4, gives minimum of trough front/back period and wave 0605 % trough amplitude (min(Ttf,Ttb),At). 0606 0607 if (ts == tn-ts+1) % then 0608 pdf(1:Nx,ts) = pdf(1:Nx,ts)+fxind.'*dt; 0609 err(1:Nx,ts) = err(1:Nx,ts)+(err0(:)*dt).^2; 0610 else 0611 pdf(1:Nx,ts) = pdf(1:Nx,ts)+fxind.'*2.0*dt; 0612 err(1:Nx,ts) = err(1:Nx,ts)+(err0(:)*2.0*dt).^2; 0613 end %if 0614 end %select 0615 end %do % ts 0616 %print *,'Ready: ',tn,' of ',Ntime, ' ABSEPS = ', ABSEPS 0617 waitTxt = sprintf('%s Ready: %d of %d',datestr(now),tn,Ntime); 0618 fwaitbar(tn/Ntime,h11,waitTxt); 0619 end %do %tn 0620 close(h11) 0621 % err0 is given as 3 standarddeviations of the variability of fxind. 0622 % thus err is given as a variance 0623 err = sqrt(err); % convert to standarddeviations 0624 end %if 0625 return % main 0626 0627 0628 0629 function BIG = covinput(BIG, R0,R1,R2,R3,R4,tn,ts) 0630 %COVINPUT Sets up the covariance matrix 0631 % 0632 % CALL BIG = covinput(BIG, R0,R1,R2,R3,R4,tn,ts) 0633 % 0634 % BIG = covariance matrix for X = [Xt,Xd,Xc] in spec2thpdf problems. 0635 % 0636 % The order of the variables in the covariance matrix 0637 % are organized as follows: 0638 % For ts>1: 0639 % ||X(t2)..X(ts),..X(tn-1)||X''(ts) X'(t1) X'(tn)||X(ts) X(t1) X(tn) X'(ts)|| 0640 % = [Xt Xd Xc] 0641 % 0642 % For ts<=1: 0643 % ||X(t2)..,..X(tn-1)||X'(t1) X'(tn)||X(t1) X(tn)|| 0644 % = [Xt Xd Xc] 0645 0646 % where 0647 % 0648 % Xt = time points in the indicator function 0649 % Xd = derivatives 0650 % Xc = variables to condition on 0651 0652 % Computations of all covariances follows simple rules: Cov(X(t),X(s)) = r(t,s), 0653 % then Cov(X'(t),X(s))=dr(t,s)/dt. Now for stationary X(t) we have 0654 % a function r(tau) such that Cov(X(t),X(s))=r(s-t) (or r(t-s) will give the same result). 0655 % 0656 % Consequently Cov(X'(t),X(s)) = -r'(s-t) = -sign(s-t)*r'(|s-t|) 0657 % Cov(X'(t),X'(s)) = -r''(s-t) = -r''(|s-t|) 0658 % Cov(X''(t),X'(s)) = r'''(s-t) = sign(s-t)*r'''(|s-t|) 0659 % Cov(X''(t),X(s)) = r''(s-t) = r''(|s-t|) 0660 % Cov(X''(t),X''(s)) = r''''(s-t) = r''''(|s-t|) 0661 0662 if nargin<8|isempty(ts) 0663 ts = -1; 0664 end 0665 if (ts<1) % THEN 0666 Ntd1 = tn; 0667 N = Ntd1+2; 0668 shft = 0; % def=1 want only crest period Tc 0669 else 0670 Ntd1 = tn+1; 0671 N = Ntd1+4; 0672 shft = 1; % def=2 or 3 want Tc Ac or Tcf, Ac 0673 end %if 0674 if tn>2 0675 %do 0676 i=1:tn-2; 0677 %cov(Xt) 0678 BIG(i,i) = toeplitz(R0(i));% cov(X(ti+1),X(tj+1)); 0679 0680 %cov(Xt,Xc) 0681 BIG(i ,Ntd1+1+shft) = R0(i+1); %cov(X(ti+1),X(t1)) 0682 BIG(tn-1-i ,Ntd1+2+shft) = R0(i+1); %cov(X(t.. ),X(tn)) 0683 0684 %Cov(Xt,Xd)=cov(X(ti+1),x(tj) 0685 BIG(i,Ntd1-1) =-R1(i+1); %cov(X(ti+1),X' (t1)) 0686 BIG(tn-1-i,Ntd1) = R1(i+1); %cov(X(ti+1),X' (tn)) 0687 %enddo 0688 end 0689 %call echo(big(1:tn,1:tn),tn) 0690 %cov(Xd) 0691 BIG(Ntd1 ,Ntd1 ) = -R2(1); 0692 BIG(Ntd1-1,Ntd1 ) = -R2(tn); %cov(X'(t1),X'(tn)) 0693 BIG(Ntd1-1,Ntd1-1) = -R2(1); 0694 0695 %cov(Xc) 0696 BIG(Ntd1+1+shft,Ntd1+1+shft) = R0(1); % cov(X(t1),X (t1)) 0697 BIG(Ntd1+1+shft,Ntd1+2+shft) = R0(tn); % cov(X(t1),X (tn)) 0698 BIG(Ntd1+2+shft,Ntd1+2+shft) = R0(1); % cov(X(tn),X (tn)) 0699 %cov(Xd,Xc) 0700 BIG(Ntd1 ,Ntd1+1+shft) = R1(tn); %cov(X'(tn),X(t1)) 0701 BIG(Ntd1 ,Ntd1+2+shft) = 0.d0; %cov(X'(tn),X(tn)) 0702 BIG(Ntd1-1,Ntd1+1+shft) = 0.d0 ; %cov(X'(t1),X(t1)) 0703 BIG(Ntd1-1,Ntd1+2+shft) =-R1(tn); %cov(X'(t1),X(tn)) 0704 0705 0706 if (ts>1) % then 0707 0708 0709 %cov(Xc) 0710 BIG(Ntd1+1,Ntd1+1) = R0(1); % cov(X(ts),X (ts) 0711 BIG(Ntd1+1,Ntd1+2) = R0(ts); % cov(X(ts),X (t1)) 0712 BIG(Ntd1+1,Ntd1+3) = R0(tn+1-ts); % cov(X(ts),X (tn)) 0713 BIG(Ntd1+1,Ntd1+4) = 0.d0; % cov(X(ts),X'(ts)) 0714 0715 BIG(Ntd1+2,Ntd1+4) = R1(ts); % cov(X(t1),X'(ts)) 0716 BIG(Ntd1+3,Ntd1+4) = -R1(tn+1-ts); %cov(X(tn),X'(ts)) 0717 BIG(Ntd1+4,Ntd1+4) = -R2(1); % cov(X'(ts),X'(ts)) 0718 0719 %cov(Xd) 0720 BIG(Ntd1-2,Ntd1-1) = -R3(ts); %cov(X''(ts),X'(t1)) 0721 BIG(Ntd1-2,Ntd1-2) = R4(1); 0722 BIG(Ntd1-2,Ntd1 ) = R3(tn+1-ts); %cov(X''(ts),X'(tn)) 0723 %cov(Xd,Xc) 0724 BIG(Ntd1 ,Ntd1+4) =-R2(tn+1-ts); %cov(X'(tn),X'(ts)) 0725 BIG(Ntd1 ,Ntd1+1) = R1(tn+1-ts); %cov(X'(tn),X (ts)) 0726 0727 BIG(Ntd1-1,Ntd1+4) =-R2(ts); %cov(X'(t1),X'(ts)) 0728 BIG(Ntd1-1,Ntd1+1) =-R1(ts); %cov(X'(t1),X (ts)) 0729 0730 BIG(Ntd1-2,Ntd1+1) = R2(1); %cov(X''(ts),X (ts) 0731 BIG(Ntd1-2,Ntd1+2) = R2(ts); %cov(X''(ts),X (t1)) 0732 BIG(Ntd1-2,Ntd1+3) = R2(tn+1-ts); %cov(X''(ts),X (tn)) 0733 BIG(Ntd1-2,Ntd1+4) = 0.d0; %cov(X''(ts),X'(ts)) 0734 %cov(Xt,Xc) 0735 if tn>2 0736 %do 0737 i = 1:tn-2; 0738 j = (i+1-ts).'; 0739 tau = abs(j)+1; 0740 BIG(i,Ntd1+1) = R0(tau); %cov(X(ti+1),X(ts)) 0741 BIG(i,Ntd1+4) = -abs(R1(tau)).*sign(R1(tau).*j); %cov(X(ti+1),X'(ts)) % check this 0742 0743 %Cov(Xt,Xd)=cov(X(ti+1),X(ts)) 0744 BIG(i,Ntd1-2) = R2(tau); %cov(X(ti+1),X''(ts)) 0745 %enddo 0746 end %tn>2 0747 end %if % ts>1 0748 0749 % make lower triangular part equal to upper 0750 %for j=1:N-1 0751 % for i=j+1:N 0752 % BIG(i,j) = BIG(j,i) 0753 % end %do 0754 %end %do 0755 lp = find(tril(ones(size(BIG)),-1)); % indices to lower triangular part 0756 BIGT = BIG.'; 0757 BIG(lp) = BIGT(lp); 0758 % 0759 return 0760 0761
Comments or corrections to the WAFO group