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: 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; 
264 tmp=loaddata('dens.out')/A;   
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