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: This function is called by:

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

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