Home > wafo > trgauss > spec2AcAt.m

# spec2AcAt

## PURPOSE

Evaluates survival function R(h1,h2)=P(Ac>h1,At>h2).

## SYNOPSIS

[f] = spec2AcAt(spec,utc,def,paramtc,paramtt,paramt,h1,h2,nit,speed,bound)

## DESCRIPTION

 SPEC2ACAT   Evaluates survival function R(h1,h2)=P(Ac>h1,At>h2).

CALL: R = spec2AcAt(S,u,def,paramtc,paramtt,paramtcc,h1,h2,nit,speed,bound);

R    =  survival function R=P(Ac>h1,At>h2)
S    = spectral density structure
u    = reference level (default the most frequently crossed level).
def   =  'Tcc', P(Ac>h1,At>h2) for waves in time, (default)
'Lcc', P(Ac>h1,At>h2) for waves in space.
paramtc  = [0 tc Ntc] defines discretization of Tc (or Lc): tc is
the longest Tc considered while Ntc is the number of
points, i.e. (Ntc-1)/tc is the sampling
frequency. paramtc= [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.
paramtt   = [0 tt Ntt] defines discretization of Tt (or Lt): tt is
the longest Tt considered while Ntt is the number of points.
paramtcc   = [0 tcc Ntcc] defines discretization of a period Tcc (or
Lcc): tcc is the longest period considered while Ntcc is
the number of points.
h1 = vector of heights of Ac crest; note  h1 > 0,
h2 = vector of heights of At trough; note  h2 > 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 and lower bounds for the distribution are computed.
[]    = default values are used.

SPEC2ACAT evaluates P(Ac>h1,At>h2), i.e. probability that crest is
higher than h1 and trough lower than h2, for waves in time or in space
in a stationary Gaussian transform process X(t) where Y(t) = g(X(t) (Y
zero-mean Gaussian with spectrum given in S).

Example:%  Very accurate approx. of R=P(Ac>h1,At>h2) (waves in
%  time with unidirectional JONSWAP spectrum) is computed by:

Hm0=7;  Tp=11;
S = jonswap(4*pi/Tp,[Hm0 Tp]);
paramtt = [0 12 51];
paramtc = paramtt;
paramtcc=[0 18 81];
h1=0.5:0.5:6; h2=h1;
F = spec2AcAt(S,[],'Tcc',paramtc,paramtc,paramtcc,h1,h2,-2);

See also  spec2cov, wnormspec, dat2tr, dat2gaus, datastructures, wavedef

## 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 spec2Acdf Evaluates cdf of crests P(Ac<=h) or troughs P(At<=h). 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 axis Control axis scaling and appearance. clf Clear current figure. clock Current date and time as date vector. contour Contour plot. 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. figure Create figure window. 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) strcmp Compare strings. subplot Create axes in tiled positions.
This function is called by:

## SOURCE CODE

001 function [f] = spec2AcAt(spec,utc,def,paramtc,paramtt,paramt,h1,h2,nit,speed,bound)
002 %SPEC2ACAT   Evaluates survival function R(h1,h2)=P(Ac>h1,At>h2).
003 %
004 %  CALL: R = spec2AcAt(S,u,def,paramtc,paramtt,paramtcc,h1,h2,nit,speed,bound);
005 %
006 %        R    =  survival function R=P(Ac>h1,At>h2)
007 %        S    = spectral density structure
008 %        u    = reference level (default the most frequently crossed level).
009 %       def   =  'Tcc', P(Ac>h1,At>h2) for waves in time, (default)
010 %                'Lcc', P(Ac>h1,At>h2) for waves in space.
011 %    paramtc  = [0 tc Ntc] defines discretization of Tc (or Lc): tc is
012 %               the longest Tc considered while Ntc is the number of
013 %               points, i.e. (Ntc-1)/tc is the sampling
014 %               frequency. paramtc= [0 10 51] implies that the crest
015 %               periods are considered at 51 linearly spaced points in
016 %               the interval [0,10], i.e. sampling frequency is 5 Hz.
017 %   paramtt   = [0 tt Ntt] defines discretization of Tt (or Lt): tt is
018 %               the longest Tt considered while Ntt is the number of points.
019 %  paramtcc   = [0 tcc Ntcc] defines discretization of a period Tcc (or
020 %               Lcc): tcc is the longest period considered while Ntcc is
021 %               the number of points.
022 %          h1 = vector of heights of Ac crest; note  h1 > 0,
023 %          h2 = vector of heights of At trough; note  h2 > 0.
024 %        nit  =  0,...,9. Dimension of numerical integration  (default 2).
025 %               -1,-2,-3,... different important sampling type integrations.
026 %       speed = defines accuraccy of calculations, by choosing different
027 %               parameters, possible values: 1,2...,9 (9 fastest, default 4).
028 %       bound = 0 the distribution is approximated (default), and if  nit<0
029 %             = 1 the upper and lower bounds for the distribution are computed.
030 %       []    = default values are used.
031 %
032 % SPEC2ACAT evaluates P(Ac>h1,At>h2), i.e. probability that crest is
033 % higher than h1 and trough lower than h2, for waves in time or in space
034 % in a stationary Gaussian transform process X(t) where Y(t) = g(X(t) (Y
035 % zero-mean Gaussian with spectrum given in S).
036 %
037 % Example:%  Very accurate approx. of R=P(Ac>h1,At>h2) (waves in
038 %         %  time with unidirectional JONSWAP spectrum) is computed by:
039 %
040 %   Hm0=7;  Tp=11;
041 %   S = jonswap(4*pi/Tp,[Hm0 Tp]);
042 %   paramtt = [0 12 51];
043 %   paramtc = paramtt;
044 %   paramtcc=[0 18 81];
045 %   h1=0.5:0.5:6; h2=h1;
046 %   F = spec2AcAt(S,[],'Tcc',paramtc,paramtc,paramtcc,h1,h2,-2);
047 %
048 % See also  spec2cov, wnormspec, dat2tr, dat2gaus, datastructures, wavedef
049
050 % Tested on : matlab 5.3
051 % History: by I. Rychlik 01.10.1998 with name wave_t.m
052 % Revised by I. R.   13.10.1999
053 % Revised by Sylvie van Iseghem 10.11.1999- adding levels of crests and troughs
054 % Revised by I. R.   07.12.1999 The upper and lower bounds to  the
055 %  density included.
056 % Revised by I. R.   17.12.1999 The case def=-1,-2 are included.
057 % Revised by es 000322. Made call with directional spectrum possible.
058 % Revised by IR 29 VI 2000, introducing 'def' to consider also waves in space.
059 %                           and removing bug from definition of transformation.
060 % Revised IR 6 II 2001 adapted for MATLAB6
061 % REvised pab Dec 2003
062 % replace tic-toc statements with clock and etime
063 startTime = clock;
064 if nargin<3|isempty(def)
065   def='tcc';
066 end
067 if strcmp('l',lower(def(1)))
068   spec=spec2spec(spec,'k1d');
069 elseif strcmp('t',lower(def(1)))
070   spec=spec2spec(spec,'freq');
071 else
072   error('Unknown def')
073 end
074
075 if nargin<11|isempty(bound)
076    bound=0;
077 end
078 [S, xl4, L0, L2, L4, L1]=wnormspec(spec);
079 A = sqrt(L0/L2);
080 SCIS=0;
081 if nargin<9|isempty(nit)
082   nit=2;
083 elseif nit<0
084  SCIS=min(abs(nit),2);
085  bound=0;
086 end
087
088 if isfield(spec,'tr')
089    g=spec.tr;
090 else
091    g=[];
092 end
093 if isempty(g)
094   g = [sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)'];
095 end
096
097 if nargin<2|isempty(utc)
098   utc_d = gaus2dat([0, 0],g); % most frequent crossed level
099   utc   = utc_d(1,2);
100 end
101
102 % transform reference level into Gaussian level
103 uu = dat2gaus([0., utc],g);
104 u  = uu(2);
105 disp(['The level u for Gaussian process = ', num2str(u)])
106
107
108 if nargin<10|isempty(speed)
109   speed=4;
110 end
111
112 if nargin<4|isempty(paramtc)
113   defnr=1;
114   paramtc=[0., 2*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2)*(0.5+erf(-sign(defnr)*u/sqrt(2))/2)), 51];
115 end
116 if nargin<5|isempty(paramtt)
117   defnr=-1;
118   paramtt=[0., 2*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2)*(0.5+erf(-sign(defnr)*u/sqrt(2))/2)), 51];
119 end
120 if nargin<6|isempty(paramt)
121    Ntime = 51;
122    t0    = 0.;
123    tn    = 1.5*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2));
124   paramt=[t0,tn,Ntime];
125   else
126    t0     = paramt(1);
127    tn     = paramt(2);
128    Ntime  = paramt(3);
129 end
130 t  = levels([0, tn/A, Ntime]); % normalized times
131 N0 = 1+ceil(t0/tn*(Ntime-1)); % the starting point to evaluate
132 if Ntime>101
133   disp('Note: nr. of wavelengths (periods) exceeds 101 (slow).')
134 end
135
136 nr = 2;
137
138 px1=gaus2dat([0., u;1, 4],g);
139 px1=abs(px1(2,2)-px1(1,2));
140 px2=gaus2dat([0., u;1, 4],g);
141 px2=abs(px2(2,2)-px2(1,2));
142
143 if nargin<7|isempty(h1)
144   Nx1 = 6;
145   h1=levels([0, px1, 6]);
146 else
147   h1=sort(abs(h1(:))); % make sure values are positive and in increasing order
148   Nx1=length(h1);
149 end
150 figure(1)
151 clf
152 subplot(221)
153 nit0=nit;
154 if nit0>=0
155   nit0=nit0+3;
156 end
157 if strcmp('l',lower(def(1)))
158   FAc=spec2Acdf(spec,utc,'Lc',paramtc,h1,nit0,speed,bound);
159 elseif strcmp('t',lower(def(1)))
160   FAc=spec2Acdf(spec,utc,'Tc',paramtc,h1,nit0,speed,bound);
161 else
162   error('Unknown def')
163 end
164 pdfplot(FAc);
165 hold on
166 if nargin<7|isempty(h2)
167   Nx2 = 6;
168   h2=levels([0, px2, 6]);
169 else
170   h2=sort(abs(h2(:))); % make sure values are positive and  increasing
171   Nx2=length(h2);
172 end
173 subplot(222)
174 if strcmp('l',lower(def(1)))
175   FAt=spec2Acdf(spec,utc,'Lt',paramtc,h1,nit0,speed,bound);
176 elseif strcmp('t',lower(def(1)))
177   FAt=spec2Acdf(spec,utc,'Tt',paramtc,h1,nit0,speed,bound);
178 else
179   error('Unknown def')
180 end
181
182 pdfplot(FAt);
183
184
185
186 Nx=Nx1*Nx2;
187
188 if Nx>101
189    disp('Note: nr. of amplitudes exceeds 101, (slow).')
190 end
191
192 %Transform amplitudes to Gaussian levels:
193 der1=ones(Nx1,1); % dh/dh=1
194 h1=reshape(h1,Nx1,1);
195 hg1=tranproc([utc+h1, der1],g);
196 der1=abs(hg1(:,2));
197 hg1=hg1(:,1); % Gaussian level
198 der2=ones(Nx2,1); % dh/dh=1
199 h2=reshape(h2,Nx2,1);
200 hg2=tranproc([utc-h2, der2],g);
201 der2=abs(hg2(:,2));
202 hg2=hg2(:,1); % Gaussian level
203
204 if exist('h.in'), delete h.in, end
205 fid=fopen('h.in','wt');
206 fprintf(fid,'%12.10f\n',hg1);
207 fprintf(fid,'%12.10f\n',hg2);
208 fclose(fid);
209
210 dt=t(2)-t(1);
211  % Calculating covariances
212 R = spec2cov2(S,nr,Ntime-1,dt);
213 %NB!!! the spec2thpdf.exe programme is very sensitive to how you interpolate
214 %      the covariances, especially where the process is very dependent
215 %      and the covariance matrix is nearly singular. (i.e. for small t
216 %      and high levels of u if Tc and low levels of u if Tt)
217 %     The best is to first interpolate through FFT with a
218 %     high rate and then interpolate with a spline to obtain the
219 %     timepoints one want. However for large t
220 %     it often suffices to interpolate linearly provided that
221 %     FFT interpolation is good eneough.
222
223 for k=0:nr
224   filename=['Cd', int2str(k), '.in'];
225   if exist(filename)
226     delete(filename)
227   end
228   fid=fopen(filename,'wt');
229   fprintf(fid,'%12.10E \n',R(:,k+1));
230   fclose(fid);
231 end
232 %SCIS=0;
233 if exist('reflev.in'), delete reflev.in, end
234 disp('writing data')
235 fid=fopen('reflev.in','wt');
236 fprintf(fid,'%12.10E \n',u);
237 fprintf(fid,'%2.0f \n',Ntime);
238 fprintf(fid,'%2.0f \n',N0);
239 fprintf(fid,'%2.0f \n',nit);
240 fprintf(fid,'%2.0f \n',speed);
241 fprintf(fid,'%2.0f \n',SCIS);
242 fprintf(fid,'%2.0f \n',10^9*rand);  % select a random seed for rind
243 fprintf(fid,'%2.0f %2.0f\n',[Nx1, Nx2]);
244 fprintf(fid,'%12.10E \n',dt);
245 fclose(fid);
246 disp('   Starting Fortran executable.')
247 if bound==0
248   dos([ wafoexepath, 'sp2tccpdf70.exe']);
249   %compiled spec2tccpdf.f with rind60.f and intmodule.f
250   %   dos([ wafoexepath, 'sp2tccpdf50.exe']); %compiled spec2tccpdf.f with rind50.f
251 else
252    dos([ wafoexepath, 'sp2tccpdf51.exe']); %compiled spec2tccpdf1.f with rind51.f
253 end
254
255
256 f=createpdf;
257 f.labx{1}='amplitude [m]';
258 f.x{1}=h1';
259 f.x{2}=h2';
260 f.nit=nit;
261 f.speed=speed;
262 f.SCIS=SCIS;
263 f.u=utc;
265 t=t'*A;
266 dt=t(2)-t(1);
267 ff1=reshape(tmp(:,1),Nx1*Nx2,length(t));
268 if bound==1
269    gg1=reshape(tmp(:,2),Nx1*Nx2,length(t));
270    fff=zeros(Nx1,Nx2);
271    ggg=fff;
272    for j=1:Nx1
273     for i=1:Nx2
274      ii1=i+(j-1)*Nx2;
275      fff(j,i)=max(dt*sum(ff1(ii1,:))+1.-FAc.f(2,j)-FAt.f(2,i),0.);
276      fff(j,i)=min(fff(j,i),1.);
277     end
278    end
279    for j=1:Nx1
280     for i=1:Nx2
281      ii1=i+(j-1)*Nx2;
282      ggg(j,i)=max(dt*sum(gg1(ii1,:))+1.-FAc.f(1,j)-FAt.f(1,i),0.);
283      ggg(j,i)=min(ggg(j,i),1.);
284     end
285    end
286    subplot(223)
287    contour(h1,h2,fff',[1, 0.9, 0.75, 0.5, 0.25, 0.1, 0.05])
288    subplot(224)
289    contour(h1,h2,ggg',[1, 0.9, 0.75, 0.5, 0.25, 0.1, 0.05])
290    f.f{1}=fff';
291    f.f{2}=ggg';
292 else
293    subplot(212)
294    fff=zeros(Nx1,Nx2);
295    for j=1:Nx1
296     for i=1:Nx2
297      ii1=i+(j-1)*Nx2;
298      fff(j,i)=max(dt*sum(ff1(ii1,:))+1.-FAc.f(j,1)-FAt.f(i,1),0.);
299      fff(j,i)=min(fff(j,i),1.);
300     end
301    end
302    f.f=fff';
303    contour(h1,h2,fff',[1, 0.9, 0.75, 0.5, 0.25, 0.1, 0.05])
304    axis('square')
305 %   pdfplot(f);
306 end
307 f.elapsedTime = etime(clock,startTime);
308
309
310
311

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