Home > wafo > trgauss > spec2cmat.m

spec2cmat

PURPOSE ^

Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough)

SYNOPSIS ^

[f, fmm] = spec2cmat(spec,utc,def,paramt,paramu,nit)

DESCRIPTION ^

 SPEC2CMAT Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) 
            
  CALL:  f   = spec2cmat(S,u,def,paramt,paramu,nit);  
   
         f    = pdf (density structure) of crests (trough) heights  
         S    = spectral density structure  
         u    = reference level (default the most frequently crossed level). 
        def   = 'Mm'  : gives maximum and the following minimum height. 
                'rfc' : gives maximum and the rainflow minimum height. 
                'AcAt': gives (crest,trough) heights  (this option needs 
                        more work).  
      paramt  = [0 tn Nt] defines discretization of half period: tn is 
                the longest period considered while Nt is the number of 
                points, i.e. (Nt-1)/tn is the sampling frequnecy.  
                paramt=[0 10 51] implies that the halfperiods are 
                considered at 51 linearly spaced points in the interval 
                [0,10], i.e. sampling frequency is 5 Hz.  
      paramu  = [u v N] defines discretization of maxima and minima ranges:  
                u is the lowest minimum considered, v the heighest 
                maximum and N is the number of levles (u,v) included.  
         nit  =  0,...,9. Dimension of numerical integration (only 
                 positive nit are allowed). (default nit=1).  
        []    = default values are used.  
   
        
   The model for loads is a stationary Gaussian transformed process X(t), 
   where  Y(t) = g(X(t)) is a zero-mean Gaussian with spectrum, S. 
    
   Note: algorithm uses Markov Chain approximation to the sequence of 
   turning points in Y.    
  
  Example: % The intensity matrix of rainflow cycles is computed by:  
        S  = jonswap;       
        L0 = spec2mom(S,1);  
        paramu = [sqrt(L0)*[-4 4] 41];  
        frfc   = spec2cmat(S,[],'rfc',[],paramu);  
   
  See also  spec2mmtpdf spec2cov, wnormspec, dat2tr, dat2gaus, wavedef

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [f, fmm] = spec2cmat(spec,utc,def,paramt,paramu,nit)  
002 %SPEC2CMAT Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) 
003 %           
004 % CALL:  f   = spec2cmat(S,u,def,paramt,paramu,nit);  
005 %  
006 %        f    = pdf (density structure) of crests (trough) heights  
007 %        S    = spectral density structure  
008 %        u    = reference level (default the most frequently crossed level). 
009 %       def   = 'Mm'  : gives maximum and the following minimum height. 
010 %               'rfc' : gives maximum and the rainflow minimum height. 
011 %               'AcAt': gives (crest,trough) heights  (this option needs 
012 %                       more work).  
013 %     paramt  = [0 tn Nt] defines discretization of half period: tn is 
014 %               the longest period considered while Nt is the number of 
015 %               points, i.e. (Nt-1)/tn is the sampling frequnecy.  
016 %               paramt=[0 10 51] implies that the halfperiods are 
017 %               considered at 51 linearly spaced points in the interval 
018 %               [0,10], i.e. sampling frequency is 5 Hz.  
019 %     paramu  = [u v N] defines discretization of maxima and minima ranges:  
020 %               u is the lowest minimum considered, v the heighest 
021 %               maximum and N is the number of levles (u,v) included.  
022 %        nit  =  0,...,9. Dimension of numerical integration (only 
023 %                positive nit are allowed). (default nit=1).  
024 %       []    = default values are used.  
025 %  
026 %       
027 %  The model for loads is a stationary Gaussian transformed process X(t), 
028 %  where  Y(t) = g(X(t)) is a zero-mean Gaussian with spectrum, S. 
029 %   
030 %  Note: algorithm uses Markov Chain approximation to the sequence of 
031 %  turning points in Y.    
032 % 
033 % Example: % The intensity matrix of rainflow cycles is computed by:  
034 %       S  = jonswap;       
035 %       L0 = spec2mom(S,1);  
036 %       paramu = [sqrt(L0)*[-4 4] 41];  
037 %       frfc   = spec2cmat(S,[],'rfc',[],paramu);  
038 %  
039 % See also  spec2mmtpdf spec2cov, wnormspec, dat2tr, dat2gaus, wavedef  
040   
041  
042 % Tested on : matlab 5.3  
043 % History: by I. Rychlik 01.10.1998 with name minmax.m  
044 % bounds by I.R. 02.01.2000.  
045 % Revised by es 000322. Made call with directional spectrum possible. 
046 % revised by ir 000612. Help and plots improved.  
047 % revised by IR removing error in transformation 29 VI 2000 
048 % revised by I.R. 01.20.2001 Change name minmax to wminmax 
049 % revised by I.R. 6 II 2001 adapted for MATLAB 6 
050 % revised pab 30nov2003 
051  
052 % TODO % AcAt option needs more work   
053 startTime = clock;  
054 [S, xl4, L0, L2, L4, L1]=wnormspec(spec); 
055  
056   
057 A = sqrt(L0/L2);  
058 SCIS=0;  
059 if nargin<6|isempty(nit)  
060   nit=1;  
061 elseif nit<0  
062   warning('Only postive nit allowed') 
063  nit=1;  
064 end  
065   
066 if isfield(spec,'tr')  
067    g = spec.tr;  
068 else  
069    g = [];  
070 end  
071 if isempty(g)  
072   g = [sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)'];  
073 end 
074 S.tr = g; 
075  
076 if nargin<7|isempty(speed)  
077    speed=5;  
078 end                  
079 if nargin<2|isempty(utc) 
080     utc_d=gaus2dat([0, 0],g); % most frequent crossed level  
081     utc=utc_d(1,2); 
082 end 
083  
084 % transform reference level into Gaussian level 
085 uu = dat2gaus([0., utc],g); 
086 u  = uu(2); 
087 disp(['The level u for Gaussian process = ', num2str(u)]) 
088  
089  
090 if nargin<4|isempty(paramt)  
091   distanceBetweenExtremes = 5*pi*sqrt(L2/L4); %(2.5 * mean distance between extremes) 
092   paramt = [0 distanceBetweenExtremes,41];    
093 end 
094 t0     = paramt(1);  
095 tn     = paramt(2);  
096 Ntime  = paramt(3);  
097 t      = levels([0, tn/A, Ntime]); % normalized times  
098  
099  
100 if nargin<3|isempty(def)  
101     defnr=0;      
102 else  
103  switch lower(def)  
104  case  'mm',    defnr = 1; 
105  case  'rfc',   defnr = 0; 
106  case  'acat',  defnr =-1;  
107  otherwise, error('Unknown def')  
108  end  
109 end                  
110  
111   
112 nr = 4;  
113   
114 if nargin<5|isempty(paramu)   
115    paramu=[-4*sqrt(L0), 4*sqrt(L0), 41];  
116 end  
117  %Transform amplitudes to Gaussian levels:     
118 h   = levels(paramu);  
119 h   = reshape(h,length(h),1);  
120 Nx  = length(h);  
121 der = ones(Nx,1);   
122 hg  = tranproc([h, der],g);  
123 der = abs(hg(:,2));  
124 hg  = hg(:,1); % Gaussian level  
125  
126 %if exist('h.in'), delete h.in, end 
127 %fid=fopen('h.in','wt'); 
128 %fprintf(fid,'%12.10f\n',hg);  
129 %fclose(fid); 
130   
131   
132  
133 %paru=paramu; 
134 % paru(1:2)=paru(1:2)/sqrt(L0); 
135  fmM = wminmax(S,nit,paramu,t); 
136   
137   
138 Htxt=['Joint density of crest and trough'];   
139 labx='crest [m]'; 
140 laby='trough [m]'; 
141 if (defnr==0) 
142    Htxt = ['Joint density of maximum and rainflow minimum']; 
143    labx='max [m]'; 
144    laby='rainflow min [m]'; 
145 end  
146 f=createpdf;  
147 f.title=Htxt;  
148 f.nit=nit;  
149 f.speed=speed;  
150 f.SCIS=SCIS; 
151 f.labx{1}=labx; 
152 f.labx{2}=laby; 
153 f.x{1}=h;  
154 fmm=createpdf;  
155 fmm.title = ['Joint density of maximum and minimum']; 
156 fmm.labx{1}='max [m]'; 
157 fmm.labx{2}='min [m]';  
158 fmm.nit=nit;  
159 fmm.speed=speed;  
160 fmm.SCIS=SCIS; 
161 dh=h(2)-h(1); 
162 f.x{1}=h; 
163 f.x{2}=h; 
164 fmm.x{1}=h;   
165 fmm.x{2}=h;   
166 ftmp=fmM; 
167  
168 for i=1:Nx  
169   ftmp(:,i)=dh*dh*ftmp(:,i);%.*der*der(i);%* sqrt(-R(1,6)/R(1,4))/2/pi; 
170 end  
171 fmm.f=ftmp; 
172        
173 if (defnr==0) 
174   f.f=fliplr(mctp2rfc(fliplr(ftmp)));%* sqrt(-R(1,6)/R(1,4))/2/pi; 
175   fmm.f=ftmp; 
176 end 
177 if (defnr==1) 
178   f.f=ftmp; 
179 end 
180 if (defnr==-1)    
181   f.f    = fliplr(mctp2tc(fliplr(ftmp),utc,paramu)); 
182   index1 = find(f.x{1}>0); 
183   index2 = find(f.x{2}<0); 
184   f.f    = flipud(f.f(index2,index1)); 
185   f.x{1} = f.x{1}(index1); 
186   f.x{2} = abs(flipud(f.x{2}(index2))); 
187 end 
188 [f.cl,f.pl] = qlevels(f.f,[10, 30, 50, 70, 90, 95, 99, 99.9],f.x{1},f.x{2}); 
189  
190 f.elapsedTime = etime(clock,startTime); 
191   
192

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