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:
 createpdf PDF class constructor dat2gaus Transforms x using the transformation g. gaus2dat Transforms xx using the inverse of transformation g. levels Calculates discrete levels given the parameter matrix. loaddata Loads a matrix from a text file. pdfplot Plot contents of pdf structures 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. 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. hold Hold current graph. int2str Convert integer to string (Fast version). isfield True if field is in structure array. lower Convert string to lowercase. num2str Convert number to string. (Fast version) plot Linear plot. strcmp Compare strings.
This function is called by:
 spec2AcAt Evaluates survival function R(h1,h2)=P(Ac>h1,At>h2).

## 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 %
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]';
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
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