Home > wafo > spec > scalespec.m

scalespec

PURPOSE ^

Scale spectral density so that the moments equals m0,m2.

SYNOPSIS ^

[Sn,mn,mom]=scalespec(So,m0n,m2n,plotflag)

DESCRIPTION ^

 SCALESPEC  Scale spectral density so that the moments equals m0,m2.  
 
   CALL: [Sn,mn,mo]=scalespec(So,m0n,m2n,plotflag);
 
         Sn  = the scaled spectrum struct
         mn  = [m0n m1n m2n m4n] the spectral moments of Sn.
         mo  = [m0o m1o m2o m4o] the spectral moments of So.
         So  = the original spectrum struct
    plotflag = 0, do not plot the normalized spectrum (default).
               1, plot the normalized spectrum.
            
   The Scaling is performed so that
 
     INT Sn(freq) dfreq = m0n  INT freq^2  Sn(freq) dfreq = m2n
 
   where  integration limits is given by freq and Sn(freq) is the 
   spectral density. freq can be frequency or wave number.
   Default values for m0n=m2n=1. The normalization is defined by:
 
       freq'=freq*A; S(freq')=S(freq)*B;
   where
       A=sqrt(m0o/m2o)/sqrt(m0n/m2n); B=m0n/(A*m0o); 
 
   If S is a directional spectrum then a normalized gravity (.g) is added
   to Sn, such that mxx normalizes to 1, as well as m0 and mtt.
   (See spec2mom for notation of moments)
 
  Example: Transform spectra from a model scale
 
    Hm0 = 0.133; Tp = 1.36;
    Sj = jonswap(linspace(0,125,1025),[Hm0,Tp,3]); 
    ch = spec2char(Sj,{'Hm0','Tm02','Ss'});
        % to the corresponding spectrum with Hm0=12 and Ss=ch(3)
    Ss=ch(3);Tm02=ch(2);Hm0b=12;
    m0n = (Hm0b/4)^2; 
    m2n = 4*pi^2*m0n*Hm0/(Tm02^2*Hm0b); 
    Sn = scalespec(Sj,m0n,m2n,1);
    ch2 = spec2char(Sn,{'Hm0','Tm02','Ss'})
 
  See also  wnormspec, spec2mom

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [Sn,mn,mom]=scalespec(So,m0n,m2n,plotflag)
002 %SCALESPEC  Scale spectral density so that the moments equals m0,m2.  
003 %
004 %  CALL: [Sn,mn,mo]=scalespec(So,m0n,m2n,plotflag);
005 %
006 %        Sn  = the scaled spectrum struct
007 %        mn  = [m0n m1n m2n m4n] the spectral moments of Sn.
008 %        mo  = [m0o m1o m2o m4o] the spectral moments of So.
009 %        So  = the original spectrum struct
010 %   plotflag = 0, do not plot the normalized spectrum (default).
011 %              1, plot the normalized spectrum.
012 %           
013 %  The Scaling is performed so that
014 %
015 %    INT Sn(freq) dfreq = m0n  INT freq^2  Sn(freq) dfreq = m2n
016 %
017 %  where  integration limits is given by freq and Sn(freq) is the 
018 %  spectral density. freq can be frequency or wave number.
019 %  Default values for m0n=m2n=1. The normalization is defined by:
020 %
021 %      freq'=freq*A; S(freq')=S(freq)*B;
022 %  where
023 %      A=sqrt(m0o/m2o)/sqrt(m0n/m2n); B=m0n/(A*m0o); 
024 %
025 %  If S is a directional spectrum then a normalized gravity (.g) is added
026 %  to Sn, such that mxx normalizes to 1, as well as m0 and mtt.
027 %  (See spec2mom for notation of moments)
028 %
029 % Example: Transform spectra from a model scale
030 %
031 %   Hm0 = 0.133; Tp = 1.36;
032 %   Sj = jonswap(linspace(0,125,1025),[Hm0,Tp,3]); 
033 %   ch = spec2char(Sj,{'Hm0','Tm02','Ss'});
034 %       % to the corresponding spectrum with Hm0=12 and Ss=ch(3)
035 %   Ss=ch(3);Tm02=ch(2);Hm0b=12;
036 %   m0n = (Hm0b/4)^2; 
037 %   m2n = 4*pi^2*m0n*Hm0/(Tm02^2*Hm0b); 
038 %   Sn = scalespec(Sj,m0n,m2n,1);
039 %   ch2 = spec2char(Sn,{'Hm0','Tm02','Ss'})
040 %
041 % See also  wnormspec, spec2mom
042  
043 % Tested on: Matlab 5.3
044 % History:
045 % revised pab jan2004  
046 % by pab 20.09.2000
047   
048 % TODO % Scaling is not correct for directional spectra.
049 % TODO % Needs testing for directional spectra.
050 
051 
052 if nargin<3|isempty(m2n),m2n=1;end
053 if nargin<2|isempty(m0n),m0n=1;end   
054 if (nargin<4)|isempty(plotflag)
055   if nargout==0,
056     plotflag=1;
057   else
058     plotflag=0;
059   end
060 end
061 
062 if strcmpi(So.type(end-2:end),'k2d')
063   intype=So.type;
064   So=spec2spec(So,'dir');
065 end
066   
067 S = So.S;   % size np x nf
068 ftype = freqtype(So);
069 switch ftype 
070 case 'w',  f=So.w(:);  
071 case 'k',  f=So.k(:);
072 otherwise
073   f=2*pi*So.f(:);
074   S=S/2/pi;
075 end
076 S1=abs(S);
077 if strcmp(So.type(end-2:end),'dir')
078    S2=trapz(So.theta(:),S1.*(cos(So.theta(:)*ones(size(f'))).^2),1).'; % integrate out theta
079    S1=trapz(So.theta(:),S1,1).'; % integrate out theta
080    m20=trapz(f,f.^4/gravity^2.*S2,1);
081 else
082    S1=S1(:);
083 end
084 mom = trapz(f,[S1  f.*S1 f.^2.*S1 f.^4.*S1] ,1);
085 m0 = mom(1); 
086 m1 = mom(2); 
087 m2 = mom(3);  
088 m4 = mom(4);
089 
090 SM0 = sqrt(m0); 
091 SM2 = sqrt(m2); 
092 SM0n=sqrt(m0n); 
093 SM2n=sqrt(m2n); 
094 A = SM0*SM2n/(SM2*SM0n);
095 B = SM2*SM0n*m0n/(SM0*m0*SM2n);
096 
097 f  = f*A;
098 S1 = S*B;
099 
100 if (nargout > 1)
101    mn=trapz(f,[S1  f.*S1 f.^2.*S1 f.^4.*S1] );
102 end
103 
104 Sn=So;
105 Sn.S=S1;
106 
107 switch ftype 
108 case 'w',  Sn.w=f; 
109 case 'k',  Sn.k=f; 
110 otherwise
111    Sn.f=f/2/pi;
112    Sn.S=S.S*2*pi;
113 end
114 
115 if strcmp(So.type,'dir')
116   Sn.g=gravity*sqrt(m0*m20)/m2;
117 end
118 
119 Sn.norm=(m0n==1 & m2n==1);
120 if Sn.norm,
121    Sn.note=[Sn.note ' Scaled'];
122 else    
123    Hm0 = 4*sqrt(m0n);Tm02=2*pi*SM0n/SM2n;
124    Sn.note=[Sn.note,' Scaled to Hm0 =' num2str(Hm0),' Tm02 =' num2str(Tm02)]; 
125 end
126    
127 %Sn.date=strvcat(datestr(now),Sn.date);
128 Sn.date=datestr(now);
129    
130 if exist('intype')
131   Sn=spec2spec(Sn,intype);
132 end  
133 if plotflag
134   wspecplot(Sn)
135 end
136

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