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)

## CROSS-REFERENCE INFORMATION

This function calls:
 createspec Spectrum structure constructor wspecplot Plot a spectral density fmin fminbnd Scalar bounded nonlinear function minimization. linspace Linearly spaced vector. num2str Convert number to string. (Fast version) str2num Convert string matrix to numeric array. version MATLAB version number.
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 %
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