Home > wafo > spec > mccormick.m

mccormick

PURPOSE ^

Calculates (and plots) a McCormick spectral density.

SYNOPSIS ^

S1=mccormick(w1,sdata,plotflag)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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