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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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

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