Home > wafo > trgauss > spec2tpdf.m

# spec2tpdf

## PURPOSE

Evaluates densities for crest-,trough-period, length.

## SYNOPSIS

[f] = spec2tpdf(spec,utc,def,paramt,h,nit,speed,bound,plotflag)

## DESCRIPTION

 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]);
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

## 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 delete Delete file or graphics object. dos Execute DOS command and return result. erf Error function. error Display message and abort function. 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. tic Start a stopwatch timer. toc Read the stopwatch timer.
This function is called by:
 Chapter1 % CHAPTER1 demonstrates some applications of WAFO Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics

## SOURCE CODE

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

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