SPEC2TPDF Evaluates densities for crest-,trough-period, length. CALL: F = spec2tpdf(S,u,def,paramt,h,nit,speed,bound,plotflag); F = density structure. S = spectral density structure u = reference level (default the most frequently crossed level). def = 'Tc', gives half wave period, Tc (default). 'Tt', gives half wave period, Tt 'Lc' and 'Lt' ditto for wave length. paramt = [t0 tn Nt] where t0, tn and Nt is 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 equidistant points in the interval [0,5]. h = amplitude condition: density for waves with crests above h. (note h >= 0), (default 0.) if abs(h) less 1% of standard deviation of process X then 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 5). bound = 0 the distribution is approximated (default) = 1 the upper and lower bounds for the distribution are computed. plotflag= if 0 then do not plot, else plot (default 0). [] = default values are used. Calculates pdf of halfperiods Tc, Tt, Lc or Lt such that Ac>h or 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). The tr. g, can be estimated using lc2tr, dat2tr, hermitetr or ochitr. Example:% For the directional sea compute density of encountered Tc in % the direction pi/4 from the principal wave direction (0) at % points 0.,0.1,...,10. Hm0=7; Tp=11; S = jonswap(6*pi/Tp,[Hm0 Tp]); D=spreading(linspace(-pi,pi,51),'cos2s',[],[],S.w); Sdir=mkdspec(S,D,1); Senc=spec2spec(Sdir,'enc',pi/8,10); paramt=[0 10 101]; f = spec2tpdf(Senc,[],'Tc',paramt,[],-1,[],[],1); hold on; f1 = spec2tpdf(Sdir,[],'Tc',paramt,[],-1,[],[],1); See also spec2cov, wnormspec, datastructures, wavedef, wafomenu
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 | |
Delete file or graphics object. | |
Execute DOS command and return result. | |
Error function. | |
Display message and abort function. | |
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. | |
Start a stopwatch timer. | |
Read the stopwatch timer. |
% CHAPTER1 demonstrates some applications of WAFO | |
% CHAPTER3 Demonstrates distributions of wave characteristics |
001 function [f] = spec2tpdf(spec,utc,def,paramt,h,nit,speed,bound,plotflag) 002 %SPEC2TPDF Evaluates densities for crest-,trough-period, length. 003 % 004 % CALL: F = spec2tpdf(S,u,def,paramt,h,nit,speed,bound,plotflag); 005 % 006 % F = density structure. 007 % S = spectral density structure 008 % u = reference level (default the most frequently crossed level). 009 % def = 'Tc', gives half wave period, Tc (default). 010 % 'Tt', gives half wave period, Tt 011 % 'Lc' and 'Lt' ditto for wave length. 012 % paramt = [t0 tn Nt] where t0, tn and Nt is 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 equidistant points in the interval [0,5]. 017 % h = amplitude condition: density for waves with crests above h. 018 % (note h >= 0), (default 0.) 019 % if abs(h) less 1% of standard deviation of process X then h=0. 020 % nit = 0,...,9. Dimension of numerical integration (default 2). 021 % -1,-2,-3,... different important sampling type integrations. 022 % speed = defines accuraccy of calculations by choosing different 023 % parameters, possible values: 1,2...,9 (9 fastest, default 5). 024 % bound = 0 the distribution is approximated (default) 025 % = 1 the upper and lower bounds for the distribution are computed. 026 % plotflag= if 0 then do not plot, else plot (default 0). 027 % [] = default values are used. 028 % 029 % Calculates pdf of halfperiods Tc, Tt, Lc or Lt such that Ac>h or At>h, 030 % in a stationary Gaussian transform process X(t), where Y(t) = g(X(t)) 031 % (Y zero-mean Gaussian with spectrum given in S). The tr. g, can be estimated 032 % using lc2tr, dat2tr, hermitetr or ochitr. 033 % 034 % Example:% For the directional sea compute density of encountered Tc in 035 % % the direction pi/4 from the principal wave direction (0) at 036 % % points 0.,0.1,...,10. 037 % 038 % Hm0=7; Tp=11; S = jonswap(6*pi/Tp,[Hm0 Tp]); 039 % D=spreading(linspace(-pi,pi,51),'cos2s',[],[],S.w); 040 % Sdir=mkdspec(S,D,1); 041 % Senc=spec2spec(Sdir,'enc',pi/8,10); paramt=[0 10 101]; 042 % f = spec2tpdf(Senc,[],'Tc',paramt,[],-1,[],[],1); 043 % hold on; f1 = spec2tpdf(Sdir,[],'Tc',paramt,[],-1,[],[],1); 044 % 045 % See also spec2cov, wnormspec, datastructures, wavedef, wafomenu 046 047 % Tested on : matlab 5.3 048 % History: by I. Rychlik 01.10.1998 with name wave_th1.m 049 % revised by Per A. Brodtkorb 19.09.1999 050 % revised by I.R. 30.09.1999, bugs removing. 051 % continued by s.v.i 10.11.1999 by adding crests level in the period densities 052 % an then calculation of crests distribution. 053 % changed name and introduced possibility of computation of upper and lower 054 % bounds by I.R. 17.12.1999. 055 % revised by es 28.01.2000 Adjusting for directional spectrum and wave length. 056 % revised by es 28.01.2000 help text 057 % revised by IR removing error in transformation 29 VI 2000 058 % revised by IR adopting to Matlab 6.0 - 6 II 2001 059 % revised pab 30nov2003 060 tic 061 if nargin<3|isempty(def) 062 def='tc'; 063 end 064 if strcmp('l',lower(def(1))) 065 spec=spec2spec(spec,'k1d'); 066 elseif strcmp('t',lower(def(1))) 067 spec=spec2spec(spec,'freq'); 068 else 069 error('Unknown def') 070 end 071 switch lower(def) 072 case {'tc','lc'}, defnr = 1; % 'tc' or 'lc' 073 case {'tt','lt'}, defnr =-1; % 'tt' or 'lt' 074 otherwise, ,error('Unknown def') 075 end 076 [S, xl4, L0, L2, L4, L1]=wnormspec(spec); 077 078 A=sqrt(L0/L2); 079 SCIS=0; 080 if nargin<6|isempty(nit) 081 nit=2; 082 elseif nit<0 083 SCIS=min(abs(nit),11); 084 nit=1; 085 end 086 087 if isfield(spec,'tr') 088 g=spec.tr; 089 else 090 g=[]; 091 end 092 if isempty(g) 093 g=[(-5:0.02:5)', (-5:0.02:5)']; 094 g(:,1)=sqrt(L0)*g(:,1); 095 end 096 097 098 if nargin<2|isempty(utc) 099 utc_d=gaus2dat([0, 0],g); % most frequent crossed level 100 utc=utc_d(1,2); 101 end 102 103 % transform reference level into Gaussian level 104 uu=dat2gaus([0., utc],g); 105 u=uu(2); 106 disp(['The level u for Gaussian process = ' num2str(u)]) 107 108 109 110 if nargin<7|isempty(speed) 111 speed=5; 112 end 113 114 if nargin<8|isempty(bound) 115 bound=0; 116 end 117 if nargin<9|isempty(plotflag) 118 plotflag=0; 119 end 120 121 if SCIS>0 122 bound=0; 123 end 124 125 if nargin<4|isempty(paramt) 126 Ntime = 51; 127 t0 = 0.; 128 tn = 2*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2)*(0.5+erf(-sign(defnr)*u/sqrt(2))/2)); 129 else 130 t0 = paramt(1); 131 tn = paramt(2); 132 Ntime = paramt(3); 133 end 134 135 t = levels([0, tn/A, Ntime]); % normalized times 136 137 N0 = 1+round(t0/tn*(Ntime-1)); % the starting point to evaluate 138 %if Ntime>101 139 % disp('nr. of wavelengths limited to 101.') 140 %end 141 142 143 nr = 2; 144 px=gaus2dat([0., u;1, 5],g); 145 px=abs(px(2,2)-px(1,2)); 146 Nx = 1; 147 if nargin<5|isempty(h) 148 h=px; 149 h0=0.; 150 else 151 h=abs(min(h)); 152 h0=h; 153 if h0>0.01*sqrt(L0) 154 Nx=2; 155 h=[h; px]; 156 else 157 h=px; 158 h0=0.; 159 end 160 end 161 162 h=reshape(h,length(h),1); 163 hg=tranproc([utc+sign(defnr)*h],g); 164 165 166 167 dt = t(2)-t(1); 168 R = spec2cov2(S,nr,Ntime-1,dt); 169 170 for k=0:nr 171 filename=['Cd', int2str(k), '.in']; 172 if exist(filename) 173 delete(filename) 174 end 175 fid = fopen(filename,'wt'); 176 fprintf(fid,'%12.10E \n',R(:,k+1)); 177 fclose(fid); 178 end 179 %SCIS=0; 180 181 if exist('h.in'), delete('h.in'), end 182 183 fid = fopen('h.in','wt'); 184 fprintf(fid,'%12.10f\n',hg); 185 fclose(fid); 186 187 188 if exist('reflev.in'), delete('reflev.in'), end 189 disp('writing data') 190 fid=fopen('reflev.in','wt'); 191 fprintf(fid,'%12.10E \n',u); 192 defnr; 193 194 fprintf(fid,'%2.0f \n',defnr); 195 fprintf(fid,'%2.0f \n',Ntime); 196 fprintf(fid,'%2.0f \n',N0); 197 fprintf(fid,'%2.0f \n',nit); 198 fprintf(fid,'%2.0f \n',speed); 199 fprintf(fid,'%2.0f \n',SCIS); 200 fprintf(fid,'%2.0f \n',10^9*rand); % select a random seed for rind 201 fprintf(fid,'%2.0f \n',Nx); 202 fprintf(fid,'%12.10E \n',dt); 203 fclose(fid); 204 disp(' Starting Fortran executable.') 205 if bound<0.5 206 dos([ wafoexepath, 'sp2Acdf70.exe']); %compiled spec2Acdf.f with rind60.f and intmodule.f 207 else 208 dos([ wafoexepath, 'sp2Acdf51.exe']); %compiled spec2Acdf1.f with rind51.f 209 end 210 211 f=createpdf; 212 switch lower(def) 213 case 'tc' 214 Htxt = ['Density of Tc with Ac>', num2str(h0), ' u=',num2str(utc)]; 215 xtxt = ['T [s]']; 216 case 'tt' 217 Htxt = ['Density of Tt with At>', num2str(h0), ' u=',num2str(utc)]; 218 xtxt = ['T [s]']; 219 case 'lc' 220 Htxt = ['Density of Lc with Ac>', num2str(h0), ' u=',num2str(utc)]; 221 xtxt = ['L [m]']; 222 case 'lt' 223 Htxt = ['Density of Lt with At>', num2str(h0), ' u=',num2str(utc)]; 224 xtxt = ['L [m]']; 225 end 226 227 f.title=Htxt; 228 f.labx{1}=xtxt; 229 ftmp=loaddata('dens.out')/A; 230 f.x{1}=t*A; 231 if bound<0.5 232 ftmp=reshape(ftmp,Nx,length(t)); 233 if (length(t)>2) 234 ftmp(:,2)=0.5*(ftmp(:,3)+ftmp(:,1)); 235 end 236 if (Nx>1) 237 ft_up=ftmp(2,1:Ntime)-ftmp(1,1:Ntime); 238 else 239 ft_up=ftmp(1,1:Ntime); 240 end 241 f.f=ft_up; 242 else 243 ftmp_up=reshape(ftmp(:,1),Nx,length(t)); 244 ftmp_lo=reshape(ftmp(:,2),Nx,length(t)); 245 if (length(t)>2) 246 ftmp_up(:,2)=0.5*(ftmp_up(:,3)+ftmp_up(:,1)); 247 ftmp_lo(:,2)=0.5*(ftmp_lo(:,3)+ftmp_lo(:,1)); 248 end 249 if (Nx>1) 250 ft_lo=max(ftmp_lo(2,1:Ntime)-ftmp_up(1,1:Ntime),0.); 251 ft_up=ftmp_up(2,1:Ntime)-ftmp_lo(1,1:Ntime); 252 else 253 ft_up=ftmp_up(1,1:Ntime); 254 ft_lo=ftmp_lo(1,1:Ntime); 255 end 256 f.f=[ft_up', ft_lo']; 257 f.x{1}=t*A; 258 end 259 f.nit=nit; 260 f.speed=speed; 261 f.SCIS=SCIS; 262 f.u=utc; 263 toc 264 if plotflag 265 pdfplot(f) 266 end 267 268
Comments or corrections to the WAFO group