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
returns the frequency type of a Spectral density struct. | |
returns the constant acceleration of gravity | |
Transforms between different types of spectra | |
Plot a spectral density | |
String representation of date. | |
Check if variables or functions are defined. | |
Current date and time as date number. | |
Convert number to string. (Fast version) | |
Compare strings. | |
Compare strings ignoring case. | |
Trapezoidal numerical integration. |
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
Comments or corrections to the WAFO group