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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

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

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