Home > wafo > spec > ohspec.m

ohspec

PURPOSE

Calculates (and plots) a Ochi-Hubble spectral density.

SYNOPSIS

S1=ohspec(w1,sdata,plotflag)

DESCRIPTION

 OHSPEC Calculates (and plots) a Ochi-Hubble spectral density.

CALL:  S = ohspec(w,data,plotflag);
S = ohspec(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 L]
Hm0 = significant wave height (default 7 (m))
Tp  = peak period (default 11 (sec))
L   = spectral shape parameter (default 3)
plotflag = 0, do not plot the spectrum (default).
1, plot the spectrum.

The OH spectrum parameterization used is

S(w) = B^L*Hm0^2/(4*gamma(L)*w^(4*L+1))*exp(-B/w^4)
where
B = (L+1/4)*(2*pi/Tp)^4

Example:
S = ohspec(0.95);

See also  ohspec2, jonswap, torsethaugen

CROSS-REFERENCE INFORMATION

This function calls:
 createspec Spectrum structure constructor wspecplot Plot a spectral density gamma Gamma function. linspace Linearly spaced vector. num2str Convert number to string. (Fast version)
This function is called by:
 ohspec2 Calculates and plots a bimodal Ochi-Hubble spectral density ohspec3 Calculates Bimodal Ochi-Hubble spectral densities

SOURCE CODE

001 function S1=ohspec(w1,sdata,plotflag)
002 %OHSPEC Calculates (and plots) a Ochi-Hubble spectral density.
003 %
004 % CALL:  S = ohspec(w,data,plotflag);
005 %        S = ohspec(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 L]
011 %               Hm0 = significant wave height (default 7 (m))
012 %               Tp  = peak period (default 11 (sec))
013 %               L   = spectral shape parameter (default 3)
014 %    plotflag = 0, do not plot the spectrum (default).
015 %               1, plot the spectrum.
016 %
017 %  The OH spectrum parameterization used is
018 %
019 %     S(w) = B^L*Hm0^2/(4*gamma(L)*w^(4*L+1))*exp(-B/w^4)
020 %  where
021 %       B = (L+1/4)*(2*pi/Tp)^4
022 %
023 %
024 % Example:
025 %  S = ohspec(0.95);
026 %
028
029 % References:
030 % Ochi, M.K. and Hubble, E.N. (1976)
031 % 'On six-parameter wave spectra.'
032 % In Proc. 15th Conf. Coastal Engng., Vol.1, pp301-328
033
034 % Tested on: matlab 6.0, 5.3
035 % History:
036 % Revised pab Apr2005
037 % Fixed bug:
038 % Revised pab Dec2004
039 %  Fixed bug: w was previously not initiliazed.
040 % Revised jr 03.04.2001
041 % - added wc to input
042 % - updated information
043 % Revised pab 24.11.2000
044 % - fixed bug: wrong sign L-0.25 is replaced with L+0.25
045 % revised pab 16.02.2000
046 % by pab 01.12.99
047
048 monitor=0;
049 w = [];
050 if nargin<3|isempty(plotflag)
051   plotflag=0;
052 end
053 % Old call
054 %if nargin<1|isempty(w)
055 %  w=linspace(0,3,257).';
056 %end
057 sdata2=[7 11 3]; % default values
058 if nargin<2|isempty(sdata)
059 else
060  ns1=length(sdata);
061  k=find(~isnan(sdata(1:min(ns1,3))));
062  if any(k)
063    sdata2(k)=sdata(k);
064  end
065 end %
066 if nargin<1|isempty(w1),
067    wc = 33/sdata2(2);
068 elseif length(w1)==1,
069    wc = w1;
070 else
071    w = w1 ;
072 end
073 nw = 257;
074 if isempty(w),
075    w = linspace(0,wc,nw).';
076 end
077 n       = length(w);
078
079 S1      = createspec;
080 S1.S    = zeros(n,1);
081 S1.w    = w(:);
082 S1.norm = 0; % The spectrum is not normalized
083
084
085 Hm0     = sdata2(1);
086 Tp      = sdata2(2);
087 L       = sdata2(3);
088 S1.note = ['Ochi-Hubble, Hm0 = ' num2str(Hm0)  ', Tp = ' num2str(Tp), ...
089     ', L = ' num2str(L)];
090 if monitor
091   disp(['Hm0, Tp      = ' num2str([Hm0 Tp])])
092 end
093 wp = 2*pi/Tp;
094 k  = find(w>0);  % avoid division by zero
095 %Old call
096 %S1.S(k)=0.25*((L-0.25)*wp^4)^L/gamma(L)*Hm0^2./(w(k).^(4*L+1)).*exp(-(L-0.25)*(wp./w(k)).^4);
097 % New call pab 24.11.2000
098 B = (L+0.25);
099 S1.S(k)=0.25*(B*wp^4)^L/gamma(L)*Hm0^2./(w(k).^(4*L+1)).*exp(-B*(wp./w(k)).^4);
100
101 % for w<=0
102
103 %S1.S(~k)=0;
104
105 if plotflag
106   wspecplot(S1,plotflag)
107 end
108

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