SPEC2TCCPDF Evaluates densities of wave period Tcc, wave lenght Lcc. CALL: f = spec2tccpdf(S,u,def,paramt,h1,h2,nit,speed,bound); f = density structures of wave periods of waves. S = spectral density structure. u = reference level (default the most frequently crossed level). def = 't<', period for waves with Ac<h1, At<h2, (default) 't>', period for waves with Ac>h1, At>h2. 'L<' and 'L>' ditto for wave length paramt = [t0 tn Nt] where t0, tn and Nt are the first value, last value and the number of points, respectively, for which the density will be computed. paramt= [5 5 51] implies that the density is computed only for T=5 and using 51 linearly spaced points in the interval [0,5]. h1 = height of Ac crest; note only one positive value h1 > 0, h2 = height of At trough; note only one positive value h2 > 0. nit = 0,...,9. Dimension of numerical integration (default 2), -1,-2,... 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) = 1 the upper, lower bounds for the distribution are computed only for NIT>=0. OBSERVE that options '>' and bound=1 the bounds converges very slowly, i.e. it requires very high NIT (options '<' and bound=1 is faster) plotflag= if 0 then do not plot, else plot (default 0) [] = default values are used. Calculates pdf of period Tcc (crest to crest period) (or wave lenght Lcc) such that Ac>h1 and At>h2, or Ac<h1 and At<h2 in a stationary Gaussian transform process X(t), where Y(t) = g(X(t)).(Y zero-mean Gaussian with spectrum given in S). The transformation, g, can be estimated using lc2tr or dat2tr. Example: %Compute the pdf of Tcc with crest and trough higher then 0.5*Hs: S=jonswap; L=spec2mom(S); NIT=0; f = spec2tccpdf(S,[],'t>',[],[2*sqrt(L(1))],[2*sqrt(L(1))],NIT); See also spec2cov2, wnormspec, dat2tr, dat2gaus, datastructures, perioddef
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 | |
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 | |
Current date and time as date vector. | |
Delete file or graphics object. | |
Execute DOS command and return result. | |
Display message and abort function. | |
Elapsed time. | |
Check if variables or functions are defined. | |
Close file. | |
Open file. | |
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. |
% CHAPTER3 Demonstrates distributions of wave characteristics |
001 function [f] = spec2tccpdf(spec,utc,def,paramt,h1,h2,nit,speed,bound,plotflag) 002 %SPEC2TCCPDF Evaluates densities of wave period Tcc, wave lenght Lcc. 003 % 004 % CALL: f = spec2tccpdf(S,u,def,paramt,h1,h2,nit,speed,bound); 005 % 006 % f = density structures of wave periods of waves. 007 % S = spectral density structure. 008 % u = reference level (default the most frequently crossed level). 009 % def = 't<', period for waves with Ac<h1, At<h2, (default) 010 % 't>', period for waves with Ac>h1, At>h2. 011 % 'L<' and 'L>' ditto for wave length 012 % paramt = [t0 tn Nt] where t0, tn and Nt are the first value, last value 013 % and the number of points, respectively, for which 014 % the density will be computed. paramt= [5 5 51] implies 015 % that the density is computed only for T=5 and 016 % using 51 linearly spaced points in the interval [0,5]. 017 % h1 = height of Ac crest; note only one positive value h1 > 0, 018 % h2 = height of At trough; note only one positive value h2 > 0. 019 % nit = 0,...,9. Dimension of numerical integration (default 2), 020 % -1,-2,... different important sampling type integrations. 021 % speed = defines accuraccy of calculations, by choosing different. 022 % parameters, possible values: 1,2,...,9 (9 fastest, default 4). 023 % bound = 0 the distribution is approximated (default) 024 % = 1 the upper, lower bounds for the distribution are 025 % computed only for NIT>=0. 026 % OBSERVE that options '>' and bound=1 the bounds 027 % converges very slowly, i.e. it requires very high NIT 028 % (options '<' and bound=1 is faster) 029 % plotflag= if 0 then do not plot, else plot (default 0) 030 % [] = default values are used. 031 % 032 % Calculates pdf of period Tcc (crest to crest period) (or wave lenght 033 % Lcc) such that Ac>h1 and At>h2, or Ac<h1 and At<h2 in a stationary 034 % Gaussian transform process X(t), where Y(t) = g(X(t)).(Y zero-mean 035 % Gaussian with spectrum given in S). The transformation, g, can be 036 % estimated using lc2tr or dat2tr. 037 % 038 % Example: %Compute the pdf of Tcc with crest and trough higher then 0.5*Hs: 039 % S=jonswap; L=spec2mom(S); NIT=0; 040 % f = spec2tccpdf(S,[],'t>',[],[2*sqrt(L(1))],[2*sqrt(L(1))],NIT); 041 % 042 % See also spec2cov2, wnormspec, dat2tr, dat2gaus, datastructures, perioddef 043 044 % Tested on : matlab 5.3 045 % History: by I. Rychlik 01.10.1998 with name wave_t.m 046 % Revised by I. R. 13.10.1999 047 % Revised by Sylvie van Iseghem 10.11.1999- adding levels of crests and troughs 048 % Revised by I. R. 07.12.1999 The upper and lower bounds to the density included. 049 % Revised by I. R. 27.12.1999 and changed name to spec2tccpdf. 050 % Revised by es 000322. Made call with directional spectrum possible. 051 % revised by IR removing error in transformation 29 VI 2000 052 % revised by IR addapting to MATLAB6.0 8 II 2001 053 % Revisde by pab 054 % replaced some code with call to spec2cov2 055 056 startTime = clock; 057 058 [S, xl4, L0, L2, L4, L1]=wnormspec(spec); 059 A=sqrt(L0/L2); 060 SCIS=0; 061 if nargin<7|isempty(nit) 062 nit=2; 063 elseif nit<0 064 SCIS=min(abs(nit),2); 065 end 066 if nargin<10|isempty(plotflag) 067 plotflag=0; 068 end 069 070 if isfield(spec,'tr') 071 g=spec.tr; 072 else 073 g=[]; 074 end 075 if isempty(g) 076 g=[sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)']; 077 end 078 079 080 if nargin<2|isempty(utc) 081 utc_d = gaus2dat([0, 0],g); % most frequent crossed level 082 utc = utc_d(1,2); 083 end 084 085 % transform reference level into Gaussian level 086 uu = dat2gaus([0., utc],g); 087 u = uu(2); 088 disp(['The level u for Gaussian process = ' num2str(u)]) 089 090 if nargin<3|isempty(def) 091 def='t>'; 092 else 093 if strcmp('l',lower(def(1))) 094 spec=spec2spec(spec,'k1d'); 095 elseif strcmp('t',lower(def(1))) 096 spec=spec2spec(spec,'freq'); 097 else 098 error('Unknown def') 099 end 100 101 switch lower(def) 102 case {'l<','t<'}, defnr = 1; 103 case {'l>','t>'}, defnr =-1; 104 otherwise, error('Unknown def') 105 end 106 end 107 108 if nargin<9|isempty(bound) 109 bound=0; 110 end 111 if SCIS>0 112 bound=0; 113 end 114 115 if bound>0 116 defnr = 2*defnr; 117 end 118 if nargin<8|isempty(speed) 119 speed=4; 120 end 121 122 if nargin<4|isempty(paramt) 123 Ntime = 61; 124 t0 = 0.; 125 tn = 2.*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2)); 126 else 127 t0 = paramt(1); 128 tn = paramt(2); 129 Ntime = paramt(3); 130 end 131 132 t = levels([0, tn/A, Ntime]); % normalized times 133 N0 = 1+ceil(t0/tn*(Ntime-1)); % the starting point to evaluate 134 if Ntime>101 135 disp('Note: nr. of wavelengths (periods) exceeds 101 (slow).') 136 end 137 138 nr = 2; 139 140 px1=gaus2dat([0., u;1, 5],g); 141 px1=abs(px1(2,2)-px1(1,2)); 142 px2=gaus2dat([0., u;1, 5],g); 143 px2=abs(px2(2,2)-px2(1,2)); 144 Nx1 = 1; 145 Nx2 = 1; 146 147 if nargin<5|isempty(h1) 148 if(defnr>0) 149 h1=px1; 150 else 151 h1=0.; 152 end 153 else 154 h1=min(h1); 155 end 156 if nargin<6|isempty(h2) 157 if(defnr>0) 158 h2=px2; 159 else 160 h2=0.; 161 end 162 else 163 h2=min(h2); 164 end 165 h01=h1; 166 h02=h2; 167 168 if defnr<0 169 if h1<px1 170 Nx1=2; 171 h1=[h1; px1]; 172 else 173 error('For option ">" the value of h1 is too large') 174 end 175 if h2<px2 176 Nx2=2; 177 h2=[h2; px2]; 178 else 179 error('For option ">" the value of h2 is too large') 180 end 181 end 182 Nx=Nx1*Nx2; 183 %Transform amplitudes to Gaussian levels: 184 der1=ones(Nx1,1); % dh/dh=1 185 hg1=tranproc([utc+h1, der1],g); 186 der1=abs(hg1(:,2)); 187 hg1=hg1(:,1); % Gaussian level 188 der2=ones(Nx2,1); % dh/dh=1 189 hg2=tranproc([utc-h2, der2],g); 190 der2=abs(hg2(:,2)); 191 hg2=hg2(:,1); % Gaussian level 192 193 if exist('h.in'), delete h.in, end 194 fid=fopen('h.in','wt'); 195 fprintf(fid,'%12.10f\n',hg1); 196 fprintf(fid,'%12.10f\n',hg2); 197 fclose(fid) 198 199 dt = t(2)-t(1); 200 R = spec2cov2(S,nr,Ntime-1,dt); 201 for k=0:nr 202 filename=['Cd', int2str(k), '.in']; 203 if exist(filename) 204 delete(filename) 205 end 206 fid=fopen(filename,'wt'); 207 fprintf(fid,'%12.10E \n',R(:,k+1)); 208 fclose(fid); 209 end 210 %SCIS=0; 211 if exist('reflev.in'), delete reflev.in, end 212 disp('writing data') 213 fid=fopen('reflev.in','wt'); 214 fprintf(fid,'%12.10E \n',u); 215 fprintf(fid,'%2.0f \n',Ntime); 216 fprintf(fid,'%2.0f \n',N0); 217 fprintf(fid,'%2.0f \n',nit); 218 fprintf(fid,'%2.0f \n',speed); 219 fprintf(fid,'%2.0f \n',SCIS); 220 fprintf(fid,'%2.0f \n',10^9*rand); % select a random seed for rind 221 fprintf(fid,'%2.0f %2.0f\n',[Nx1 Nx2]); 222 fprintf(fid,'%12.10E \n',dt); 223 fclose(fid); 224 disp(' Starting Fortran executable.') 225 if abs(defnr)<1.5 226 dos([ wafoexepath 'sp2tccpdf70.exe']); 227 else 228 dos([ wafoexepath 'sp2tccpdf51.exe']); 229 end 230 231 232 f=createpdf; 233 234 switch lower(def) 235 case 't>' 236 Htxt = ['Density of Tcc with Ac>', num2str(h01), ' and At>', num2str(h02)]; 237 xtxt = ['T [s]']; 238 case 'l>' 239 Htxt = ['Density of Lcc with Ac>', num2str(h01), ' and At>', num2str(h02)]; 240 xtxt = ['L [m]']; 241 case 't<' 242 Htxt = ['Density of Tcc with Ac<', num2str(h01), ' and At<', num2str(h02)]; 243 xtxt = ['T [s]']; 244 case 'l<' 245 Htxt = ['Density of Lcc with Ac<', num2str(h01), ' and At<', num2str(h02)]; 246 xtxt = ['L [m]']; 247 end 248 249 f.title = Htxt; 250 f.labx{1} = xtxt; 251 t=t'*A; 252 f.x{1}=t; 253 f.nit=nit; 254 f.speed=speed; 255 f.SCIS=SCIS; 256 f.u=utc; 257 tmp=loaddata('dens.out')/A; 258 259 if defnr==1 260 f.f=tmp(:,1); 261 end 262 if defnr==2 263 f.f=tmp(:,1:2); 264 end 265 if defnr==-2 266 fup=reshape(tmp(:,1),4,length(t)); 267 flo=reshape(tmp(:,2),4,length(t)); 268 ff=fup(4,:)+fup(1,:)-flo(2,:)-flo(3,:); 269 gg=max(flo(4,:)+flo(1,:)-fup(2,:)-fup(3,:),0.); 270 f.f=[ff', gg']'; 271 %f.f=[fup;flo]; 272 end 273 if defnr==-1 274 fup=reshape(tmp(:,1),4,length(t)); 275 f.f=fup(4,:)+fup(1,:)-fup(2,:)-fup(3,:); 276 end 277 278 f.elapsedTime = etime(clock,startTime) 279 280 if plotflag 281 pdfplot(f) 282 end 283 284 285 286 287
Comments or corrections to the WAFO group