Home > wafo > trgauss > spec2thpdf.m

spec2thpdf

PURPOSE

Joint density of amplitude and period/wave-length characteristics

SYNOPSIS

f = spec2thpdf(spec,utc,def,paramt,h,options,plotflag)

DESCRIPTION

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

CROSS-REFERENCE INFORMATION

This function calls:
 createpdf PDF class constructor dat2gaus Transforms x using the transformation g. fwaitbar Fast display of wait bar. gaus2dat Transforms xx using the inverse of transformation g. initoptions Initializes RIND options according to speed. levels Calculates discrete levels given the parameter matrix. pdfplot Plot contents of pdf structures qlevels Calculates quantile levels which encloses P% of PDF 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. erf Error function. 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. 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. warning Display warning message; disable or enable warning messages.
This function is called by:
 Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics wafofig3 Probability density distributions (pdf) of wave period, Tt, wafofig5 Joint distribution (pdf) of crest front velocity and wave height: wafofig6 Joint distribution (pdf) of crest front period, Tcf, and crest amplitude, Ac wafofig7 Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac wafofig8 Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac for extremal waves

SOURCE CODE

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
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 %
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.
0067 % revised pab 28.03.2000
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
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
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 %
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

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