Home > wafo > trgauss > spec2Acdf.m

spec2Acdf

PURPOSE ^

Evaluates cdf of crests P(Ac<=h) or troughs P(At<=h).

SYNOPSIS ^

[f] = spec2Acdf(spec,utc,def,paramt,h,nit,speed,bound)

DESCRIPTION ^

 SPEC2ACDF  Evaluates cdf of crests P(Ac<=h) or troughs P(At<=h). 
           
   CALL:  F   = spec2Acdf(S,u,def,param,h,nit,speed,bound); 
  
         F    = cumulative distribution function of crests or trough heights 
         S    = spectral density structure 
         u    = reference level (default the most frequently crossed level). 
        def   = 'Tc',    gives crest height in time, Ac (default). 
                'Tt',    gives trough height in time, At. 
                'Lc' and 'Lt' crest trough in space 
      param  =  [0 tn Nt] defines discretization of crest period (or 
                crest length): tn is the longest period considered while 
                Nt is the number of points, i.e. (Nt-1)/tn is the 
                sampling frequnecy. param = [0 10 51] implies that the 
                crest periods are considered at 51 linearly spaced points 
                in the interval [0,10], i.e. sampling frequency is 5 Hz. 
           h  = vector of  amplitudes; note  h >= 0, 
         nit  =  0,...,9. Dimension of numerical integration  (default 2). 
                -1,-2,-3,... different important sampling type integrations. 
        speed = defines accuraccy of calculations, by choosing different  
                parameters, possible values: 1,2...,9 (9 fastest, default 4). 
        bound = 0 the distribution is approximated (default) and if NIT<0 
              = 1 the upper, lower bounds for the distribution are 
                  computed only for NIT>=0.   
        []    = default values are used. 
  
  SPEC2ACDF calculates P(Ac<=h) or P(At<=h) in a stationary Gaussian 
  transform process X(t) where Y(t) = g(X(t)) (Y zero-mean Gaussian with 
  spectrum given in S).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [f] = spec2Acdf(spec,utc,def,paramt,h,nit,speed,bound) 
0002 %SPEC2ACDF  Evaluates cdf of crests P(Ac<=h) or troughs P(At<=h). 
0003 %          
0004 %  CALL:  F   = spec2Acdf(S,u,def,param,h,nit,speed,bound); 
0005 % 
0006 %        F    = cumulative distribution function of crests or trough heights 
0007 %        S    = spectral density structure 
0008 %        u    = reference level (default the most frequently crossed level). 
0009 %       def   = 'Tc',    gives crest height in time, Ac (default). 
0010 %               'Tt',    gives trough height in time, At. 
0011 %               'Lc' and 'Lt' crest trough in space 
0012 %     param  =  [0 tn Nt] defines discretization of crest period (or 
0013 %               crest length): tn is the longest period considered while 
0014 %               Nt is the number of points, i.e. (Nt-1)/tn is the 
0015 %               sampling frequnecy. param = [0 10 51] implies that the 
0016 %               crest periods are considered at 51 linearly spaced points 
0017 %               in the interval [0,10], i.e. sampling frequency is 5 Hz. 
0018 %          h  = vector of  amplitudes; note  h >= 0, 
0019 %        nit  =  0,...,9. Dimension of numerical integration  (default 2). 
0020 %               -1,-2,-3,... different important sampling type integrations. 
0021 %       speed = defines accuraccy of calculations, by choosing different  
0022 %               parameters, possible values: 1,2...,9 (9 fastest, default 4). 
0023 %       bound = 0 the distribution is approximated (default) and if NIT<0 
0024 %             = 1 the upper, lower bounds for the distribution are 
0025 %                 computed only for NIT>=0.   
0026 %       []    = default values are used. 
0027 % 
0028 % SPEC2ACDF calculates P(Ac<=h) or P(At<=h) in a stationary Gaussian 
0029 % transform process X(t) where Y(t) = g(X(t)) (Y zero-mean Gaussian with 
0030 % spectrum given in S). 
0031  
0032 % Example: % The upper and lower bounds for F=P(Ac<hi; hi in h ),with 
0033 %     S = jonswap; 
0034 %     h = 0:0.1:5; paramt=[0 12 51];% is computed by: 
0035 %             
0036 %     F = spec2Acdf(S,[],'tc',paramt,h,[],[],1); 
0037 % 
0038 % See also  spec2cov, wnormspec, dat2tr, dat2gaus, specdef, wavedef 
0039  
0040 % Tested on : matlab 5.3 
0041 % History: by I. Rychlik 01.10.1998 with name wave_th1.m 
0042 % revised by Per A. Brodtkorb 19.09.1999 
0043 % revised by I.R. 30.09.1999, bugs removing. 
0044 % continued by s.v.i 10.11.1999 by adding crests level in the period densities 
0045 % an then calculation of crests distribution. 
0046 % changed name and introduced possibility of computation of upper and lower 
0047 % bounds by I.R. 17.12.1999. 
0048 % Computation of distribution on denser grid I.R. 26 jan. 2000 
0049 % Revised by es 000322. Made call with directional spectrum possible. 
0050 % revised by IR 6 II 2001 adapted to MATLAB6 
0051 % Revised pab Dec 2003 
0052 %  replaced code with call to spec2cov2 
0053 %  removed tic and toc   
0054    
0055 % TODO % Needs further testing   
0056    
0057 startTime = clock; 
0058 if nargin<3|isempty(def) 
0059   def='tc'; 
0060 end 
0061 if strcmp('l',lower(def(1))) 
0062   spec=spec2spec(spec,'k1d'); 
0063 elseif strcmp('t',lower(def(1))) 
0064   spec=spec2spec(spec,'freq'); 
0065 else 
0066   error('Unknown def') 
0067 end 
0068 switch lower(def) 
0069    case  {'tc','lc'},    defnr = 1; % 'tc' or 'lc' 
0070    case  {'tt','lt'},    defnr =-1; % 'tt' or 'lt' 
0071    otherwise, ,error('Unknown def') 
0072 end                 
0073  
0074 [S, xl4, L0, L2, L4, L1]=wnormspec(spec); 
0075  
0076 A = sqrt(L0/L2); 
0077 SCIS=0; 
0078 if nargin<6|isempty(nit) 
0079   nit=2; 
0080 elseif nit<0 
0081  SCIS=min(abs(nit),2); 
0082  nit=1; 
0083 end 
0084  
0085 if isfield(spec,'tr') 
0086    g=spec.tr; 
0087 else 
0088    g=[]; 
0089 end 
0090 if isempty(g) 
0091   g=[sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)']; 
0092 end 
0093  
0094  
0095  
0096 if nargin<2|isempty(utc) 
0097   utc_d = gaus2dat([0, 0],g); % most frequent crossed level  
0098   utc   = utc_d(1,2); 
0099 end 
0100  
0101 % transform reference level into Gaussian level 
0102 uu = dat2gaus([0., utc],g); 
0103 u  = uu(2); 
0104 disp(['The level u for Gaussian process = ', num2str(u)]) 
0105  
0106  
0107  
0108 if nargin<7|isempty(speed) 
0109   speed=4; 
0110 end                 
0111    
0112 if nargin<8|isempty(bound) 
0113   bound=0; 
0114 end                 
0115    
0116 if SCIS>0 
0117   bound=0; 
0118 end 
0119 ttn    = 1.2*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2)*(0.5+erf(-sign(defnr)*u/sqrt(2))/2)); 
0120 if nargin<4|isempty(paramt) 
0121   Ntime = 51; 
0122   t0    = ttn; 
0123   tn    = 2*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2)*(0.5+erf(-sign(defnr)*u/sqrt(2))/2)); 
0124   paramt=[t0,tn,Ntime]; 
0125 else 
0126   t0     = paramt(1); 
0127   tn     = paramt(2); 
0128   Ntime  = paramt(3); 
0129 end 
0130  
0131 t0=ttn;  
0132 t0=ttn/A; 
0133 tt = tn/A; 
0134 t  = levels([0, tt, Ntime]); % normalized times 
0135 N0 = 2+ceil(ttn/tn*(Ntime-1)); % the starting point to evaluate 
0136 NN0=N0; 
0137 if Ntime>101 
0138   disp('nr. of wavelengths >  101., slow') 
0139 end 
0140   
0141  
0142 nr = 2; 
0143  
0144 if nargin<5|isempty(h)  
0145     hg=[0:0.25:5.]; 
0146     h=tranproc(u+sign(defnr)*hg',fliplr(g))-utc; 
0147 end 
0148 % make sure values are positive and in increasing order 
0149 %Transform h to Gaussian levels:    
0150 h=reshape(h,length(h),1); 
0151 Nx=length(h); 
0152 if Nx>101 
0153   error('nr. of amplitudes limited to 101.') 
0154 end 
0155 hg=tranproc([utc+sign(defnr)*h],g); 
0156  
0157 if exist('h.in'), delete h.in, end 
0158 fid=fopen('h.in','wt'); 
0159 fprintf(fid,'%12.10f\n',hg); 
0160 fclose(fid); 
0161  
0162  
0163 dt = t(2)-t(1); 
0164 % Calculating covariances 
0165 R = spec2cov2(S,nr,Ntime-1,dt); 
0166  
0167 %NB!!! the spec2thpdf.exe programme is very sensitive to how you interpolate  
0168 %      the covariances, especially where the process is very dependent 
0169 %      and the covariance matrix is nearly singular. (i.e. for small t 
0170 %      and high levels of u if Tc and low levels of u if Tt) 
0171 %     The best is to first interpolate through FFT with a 
0172 %     high rate and then interpolate with a spline to obtain the 
0173 %     timepoints one want. However for large t 
0174 %     it often suffices to interpolate linearly provided that 
0175 %     FFT interpolation is good eneough. 
0176  
0177 for k=0:nr 
0178    filename=['Cd' int2str(k) '.in']; 
0179    if exist(filename) 
0180       delete(filename) 
0181    end 
0182    fid=fopen(filename,'wt'); 
0183    fprintf(fid,'%12.10E \n',R(:,k+1)); 
0184    fclose(fid); 
0185 end 
0186 %SCIS=0; 
0187  
0188 if exist('reflev.in'), delete reflev.in, end 
0189 disp('writing data') 
0190 fid=fopen('reflev.in','wt'); 
0191 fprintf(fid,'%12.10E \n',u); 
0192 defnr; 
0193  
0194 fprintf(fid,'%2.0f \n',defnr); 
0195 fprintf(fid,'%2.0f \n',Ntime); 
0196 fprintf(fid,'%2.0f \n',N0); 
0197 fprintf(fid,'%2.0f \n',nit); 
0198 fprintf(fid,'%2.0f \n',speed); 
0199 fprintf(fid,'%2.0f \n',SCIS); 
0200 fprintf(fid,'%2.0f \n',10^9*rand);  % select a random seed for rind  
0201 fprintf(fid,'%2.0f \n',Nx); 
0202 fprintf(fid,'%12.10E \n',dt); 
0203 fclose(fid);  
0204 disp('   Starting Fortran executable.') 
0205  
0206 if bound<0.5 
0207    dos([ wafoexepath, 'sp2Acdf70.exe']);  
0208  else 
0209    dos([ wafoexepath, 'sp2Acdf51.exe']);  
0210 end 
0211  
0212  
0213 switch defnr 
0214   case   1, Htxt = ['Distribution of crest height']; 
0215   case  -1, Htxt = ['Distribution of trough height']; 
0216 end  
0217  
0218 f = createpdf; 
0219 f.title=Htxt; 
0220 f.nit=nit; 
0221 f.speed=speed; 
0222 f.SCIS=SCIS; 
0223 f.u=utc; 
0224 f.labx1='amplitude [m]'; 
0225 ftmp=loaddata('dens.out')/A; 
0226 f.x{1}=h; 
0227  
0228  
0229  
0230 if bound<0.5 
0231   ftmp=reshape(ftmp,Nx,length(t)); 
0232   if (length(t)>2) 
0233     ftmp(:,2)=0.5*(ftmp(:,3)+ftmp(:,1)); 
0234    end 
0235 else 
0236  ftmp_up=reshape(ftmp(:,1),Nx,length(t)); 
0237  ftmp_lo=reshape(ftmp(:,2),Nx,length(t)); 
0238    if (length(t)>2) 
0239      ftmp_up(:,2)=0.5*(ftmp_up(:,3)+ftmp_up(:,1)); 
0240      ftmp_lo(:,2)=0.5*(ftmp_lo(:,3)+ftmp_lo(:,1)); 
0241    end 
0242 end 
0243  
0244  
0245  
0246 % 
0247 %% Here we will compute the beginning of the distribution 
0248 %% om the denser grid. 
0249 % 
0250  
0251 Ntime = 61; 
0252 t0    = 0.; 
0253 tn    = ttn; 
0254 paramt1=[t0, tn, Ntime]; 
0255 tt = tn/A; 
0256 t  = levels([0, tt, Ntime]); % normalized times 
0257 N0 = 1+ceil(t0/tn*(Ntime-1)); % the starting point to evaluate 
0258 dt=t(2)-t(1); 
0259 R = spec2cov2(S,nr,Ntime-1,dt); 
0260 for k=0:nr 
0261    filename=['Cd', int2str(k), '.in']; 
0262    if exist(filename) 
0263       delete(filename) 
0264    end 
0265    fid = fopen(filename,'wt'); 
0266    fprintf(fid,'%12.10E \n',R(:,k+1)); 
0267    fclose(fid); 
0268 end 
0269 if exist('reflev.in'), delete reflev.in, end 
0270 disp('writing data') 
0271 fid=fopen('reflev.in','wt'); 
0272 fprintf(fid,'%12.10E \n',u); 
0273 defnr; 
0274  
0275 fprintf(fid,'%2.0f \n',defnr); 
0276 fprintf(fid,'%2.0f \n',Ntime); 
0277 fprintf(fid,'%2.0f \n',N0); 
0278 fprintf(fid,'%2.0f \n',nit); 
0279 fprintf(fid,'%2.0f \n',speed); 
0280 fprintf(fid,'%2.0f \n',SCIS); 
0281 fprintf(fid,'%2.0f \n',10^9*rand);  % select a random seed for rind  
0282 fprintf(fid,'%2.0f \n',Nx); 
0283 fprintf(fid,'%12.10E \n',dt); 
0284 fclose(fid);  
0285 disp('   Starting Fortran executable.') 
0286 if bound<0.5 
0287    dos([ wafoexepath, 'sp2Acdf70.exe']);  
0288  else 
0289    dos([ wafoexepath, 'sp2Acdf51.exe']);  
0290 end 
0291  
0292  
0293 ftmp1=loaddata('dens.out')/A; 
0294  
0295 if bound<0.5  
0296    ftmp1=reshape(ftmp1,Nx,length(t)); 
0297    if (length(t)>2) 
0298      ftmp1(:,2)=0.5*(ftmp1(:,3)+ftmp1(:,1)); 
0299    end 
0300 else 
0301  ftmp_up1=reshape(ftmp1(:,1),Nx,length(t)); 
0302  ftmp_lo1=reshape(ftmp1(:,2),Nx,length(t)); 
0303    if (length(t)>2) 
0304      ftmp_up1(:,2)=0.5*(ftmp_up1(:,3)+ftmp_up1(:,1)); 
0305      ftmp_lo1(:,2)=0.5*(ftmp_lo1(:,3)+ftmp_lo1(:,1)); 
0306    end 
0307 end 
0308  
0309  
0310 X=[levels(paramt1)]; 
0311 Nt=paramt(3); 
0312 if (NN0<Nt+1) 
0313     t=levels([0, paramt(2), paramt(3)]); 
0314     X=[t(NN0:Nt), X]; 
0315     N=length(X); 
0316     [X, indX]=sort(X); 
0317     dX=X; 
0318     dx(1)=0.; 
0319     dx(N)=X(N)-X(N-1); 
0320     dx(2:N-1)=0.5*(X(3:N)-X(1:N-2)); 
0321     if bound<0.5 
0322        Y=[ftmp(:,NN0:Nt), ftmp1];  
0323        Y=Y(:,indX); 
0324        f.f=min(Y*dx',1); 
0325 %       f.f=Y; 
0326 %       f.x{1}=X; 
0327     else 
0328        Y_up=[ftmp_up(:,NN0:Nt), ftmp_up1]; 
0329        Y_up=Y_up(:,indX); 
0330        Y_lo=[ftmp_lo(:,NN0:Nt), ftmp_lo1]; 
0331        Y_lo=Y_lo(:,indX); 
0332        f.f=min([Y_up*dx', Y_lo*dx'],1); 
0333 %        f.f=[Y_up; Y_lo];        
0334 %        f.x{1}=X; 
0335     end 
0336 else 
0337     N=length(X); 
0338     dX=X; 
0339     dx(1)=0.; 
0340     dx(N)=X(N)-X(N-1); 
0341     dx(2:N-1)=0.5*(X(3:N)-X(1:N-2)); 
0342     if bound<0.5 
0343        Y=[ftmp1];  
0344        f.f=min(Y*dx',1); 
0345 %       f.f=Y; 
0346 %       f.x{1}=X; 
0347     else 
0348        Y_up=[ftmp_up1]; 
0349        Y_lo=[ftmp_lo1]; 
0350        f.f=min([Y_up*dx', Y_lo*dx'],1); 
0351     end 
0352 end 
0353 f.elapsedTime = etime(clock,startTime); 
0354  
0355 pdfplot(f) 
0356 hold on 
0357 if abs(u)<0.1 
0358  xu=tranproc(u+f.x{1},g);  
0359  plot(f.x{1},1-exp(-xu.*xu/2),'b.') 
0360  plot(f.x{1},1-exp(-xu.*xu/2)) 
0361 end  
0362

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