MCCORMICK Calculates (and plots) a McCormick spectral density. CALL: S = mccormick(w,data,plotflag); S = mccormick(wc,data,plotflag); S = a struct containing the spectral density, see datastructures. w = angular frequency (default linspace(0,33/Tp,257)) wc = angular cutoff frequency (default 33/Tp) data = [Hm0 Tp Tz] Hm0 = significant wave height (default 7 (m)) Tp = peak period (default 11 (sec)) Tz = zero-down crossing period (default 0.8143*Tp) plotflag = 0, do not plot the spectrum (default). 1, plot the spectrum. The McCormick spectrum parameterization used is S(w) = (M+1)*(Hm0/4)^2/wp*(wp./w)^(M+1)*exp(-(M+1)/M*(wp/w)^M); where Tp/Tz=(1+1/M)^(1/M)/gamma(1+1/M) Example: S = mccormick(1.1,[6.5 10]), wspecplot(S) See also jonswap, torsethaugen, simpson
Spectrum structure constructor | |
Plot a spectral density | |
Scalar bounded nonlinear function minimization. | |
Linearly spaced vector. | |
Convert number to string. (Fast version) | |
Convert string matrix to numeric array. | |
MATLAB version number. |
001 function S1=mccormick(w1,sdata,plotflag) 002 %MCCORMICK Calculates (and plots) a McCormick spectral density. 003 % 004 % CALL: S = mccormick(w,data,plotflag); 005 % S = mccormick(wc,data,plotflag); 006 % 007 % S = a struct containing the spectral density, see datastructures. 008 % w = angular frequency (default linspace(0,33/Tp,257)) 009 % wc = angular cutoff frequency (default 33/Tp) 010 % data = [Hm0 Tp Tz] 011 % Hm0 = significant wave height (default 7 (m)) 012 % Tp = peak period (default 11 (sec)) 013 % Tz = zero-down crossing period (default 0.8143*Tp) 014 % plotflag = 0, do not plot the spectrum (default). 015 % 1, plot the spectrum. 016 % 017 % The McCormick spectrum parameterization used is 018 % 019 % S(w) = (M+1)*(Hm0/4)^2/wp*(wp./w)^(M+1)*exp(-(M+1)/M*(wp/w)^M); 020 % where 021 % Tp/Tz=(1+1/M)^(1/M)/gamma(1+1/M) 022 % 023 % 024 % Example: 025 % S = mccormick(1.1,[6.5 10]), wspecplot(S) 026 % 027 % See also jonswap, torsethaugen, simpson 028 029 030 % References: 031 % M.E. McCormick (1999) 032 % "Application of the Generic Spectral Formula to Fetch-Limited Seas" 033 % Marine Technology Society, Vol 33, No. 3, pp 27-32 034 035 036 % Tested on: matlab 6.0, 5.3 037 % History: 038 % revised pab nov 2004 039 % -replaced fmin with fminbnd 040 % revised jr 03.04.2001 041 % - added wc to input 042 % - updated information 043 % by pab 01.12.99 044 045 046 monitor=0; 047 048 if nargin<3|isempty(plotflag) 049 plotflag=0; 050 end 051 if nargin<2|isempty(sdata) 052 sdata=[7 11 11/1.228]; 053 else 054 switch length(sdata) 055 case 1, sdata=[sdata, 11, 11/1.228]; 056 case 2, sdata=[sdata, 11/1.228]; 057 end 058 end 059 060 w = []; 061 if nargin<1|isempty(w1), wc = 33/sdata(2); 062 elseif length(w1)==1, wc = w1; 063 else w = w1 ; end 064 nw = 257; 065 if isempty(w), w = linspace(0,wc,nw).'; end 066 067 % Old definition of w 068 %if nargin<1|isempty(w) 069 % w=linspace(0,33/sdata(2),257).'; 070 %end 071 072 073 n=length(w); 074 S1=createspec; 075 S1.S=zeros(n,1); 076 S1.w=w; 077 S1.norm=0; % The spectrum is not normalized 078 079 080 Hm0=sdata(1); 081 if 0, 082 Tz=sdata(2); 083 Tp=sdata(3); 084 S1.note=['McCormick, Hm0, Hm0 = ' num2str(Hm0) ', Tm01 = ' num2str(Tp)]; 085 else 086 Tp=sdata(2); 087 Tz=sdata(3); 088 S1.note=['McCormick, Hm0 = ' num2str(Hm0) ', Tp = ' num2str(Tp) ', Tz = ' num2str(Tz) ]; 089 end 090 wp=2*pi/Tp; 091 092 if monitor 093 disp(['Hm0, Tp, Tz = ' num2str([Hm0 Tp Tz])]) 094 end 095 096 mvrs=version;ix=find(mvrs=='.'); 097 if str2num(mvrs(1:ix(2)-1))>5.2, 098 M=fminbnd(['((1+1./x).^(1./x)/gamma(1+1./x)-' num2str(Tp/Tz) ').^2' ],1,7); 099 else 100 M=fmin(['((1+1./x).^(1./x)/gamma(1+1./x)-' num2str(Tp/Tz) ').^2' ],1,7); 101 end 102 103 104 105 % for w>0 % avoid division by zero 106 k=find(w>0); 107 S1.S(k)=(M+1)*(Hm0/4)^2/wp*(wp./w(k)).^(M+1).*exp(-(M+1)/M*(wp./w(k)).^M); 108 109 110 if plotflag 111 wspecplot(S1,plotflag) 112 end 113
Comments or corrections to the WAFO group