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'})

## CROSS-REFERENCE INFORMATION

This function calls:
 freqtype returns the frequency type of a Spectral density struct. gravity returns the constant acceleration of gravity spec2spec Transforms between different types of spectra wspecplot Plot a spectral density datestr String representation of date. exist Check if variables or functions are defined. now Current date and time as date number. num2str Convert number to string. (Fast version) strcmp Compare strings. strcmpi Compare strings ignoring case. trapz Trapezoidal numerical integration.
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 %
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