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
Spectrum structure constructor | |
Plot a spectral density | |
Gamma function. | |
Linearly spaced vector. | |
Convert number to string. (Fast version) |
Calculates and plots a bimodal Ochi-Hubble spectral density | |
Calculates Bimodal Ochi-Hubble spectral densities |
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 % 027 % See also ohspec2, jonswap, torsethaugen 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
Comments or corrections to the WAFO group