Home > wafo > trgauss > spec2tccpdf.m

# spec2tccpdf

## PURPOSE

Evaluates densities of wave period Tcc, wave lenght Lcc.

## SYNOPSIS

[f] = spec2tccpdf(spec,utc,def,paramt,h1,h2,nit,speed,bound,plotflag)

## DESCRIPTION

 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

## 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 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 clock Current date and time as date vector. delete Delete file or graphics object. dos Execute DOS command and return result. error Display message and abort function. etime Elapsed time. exist Check if variables or functions are defined. fclose Close file. fopen Open file. 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.
This function is called by:
 Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics

## SOURCE CODE

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 %
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;
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

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