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
PDF class constructor | |
Transforms x using the transformation g. | |
Transforms xx using the inverse of transformation g. | |
Calculates discrete levels given the parameter matrix. | |
Loads a matrix from a text file. | |
Plot contents of pdf structures | |
Evaluates cdf of crests P(Ac<=h) or troughs P(At<=h). | |
Computes covariance function and its derivatives, alternative version | |
Transforms between different types of spectra | |
Transforms process X and up to four derivatives | |
Returns the path to executables for the WAFO Toolbox | |
Normalize a spectral density such that m0=m2=1 | |
Control axis scaling and appearance. | |
Clear current figure. | |
Current date and time as date vector. | |
Contour plot. | |
Delete file or graphics object. | |
Execute DOS command and return result. | |
Error function. | |
Display message and abort function. | |
Elapsed time. | |
Check if variables or functions are defined. | |
Close file. | |
Create figure window. | |
Open file. | |
Hold current graph. | |
Convert integer to string (Fast version). | |
True if field is in structure array. | |
Convert string to lowercase. | |
Convert number to string. (Fast version) | |
Compare strings. | |
Create axes in tiled positions. |
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
Comments or corrections to the WAFO group