OHHSPDF Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. CALL: f = ohhspdf(Hd,Scf,Hm0,def) f = pdf Hd = zero down crossing wave height Scf = crest front steepness Hm0 = significant wave height [m]. def = defines the parametrization of the spectral density (default 1) 1 : The most probable spectrum (default) 2,3,...11 : gives 95% Confidence spectra OHHSPDF approximates the joint distribution of (Scf, Hd) in time, i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process with a bimodal Ochi-Hubble spectral density (ohspec2). The empirical parameters of the model is fitted by least squares to simulated (Scf,Hd) data for 24 classes of Hm0. Between 50000 and 150000 zero-downcrossing waves were simulated for each class of Hm0. OHHSPDF is restricted to the following range for Hm0: 0.5 < Hm0 [m] < 12 The size of f is the common size of the input arguments, Hd, Scf and Hm0. Example: Hm0 = 6;Tp = 8;def= 2; h = linspace(0,4*Hm0/sqrt(2))'; v = linspace(0,6*1.25*Hm0/Tp^2)'; [V,H] = meshgrid(v,h); f = ohhspdf(H,V,Hm0,def); contourf(v,h,f) See also ohhpdf, thspdf
PDF class constructor | |
Wave height, Hd, distribution parameters for Ochi-Hubble spectra. | |
Calculates quantile levels which encloses P% of PDF | |
Calculates the difference between the maximum and minimum values. | |
Calculates a smoothing spline. | |
Generalized Gamma probability density function | |
Weibull cumulative distribution function | |
Weibull probability density function | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
1-D interpolation (table lookup) | |
2-D interpolation (table lookup). | |
X and Y arrays for 3-D plots. | |
Convert number to string. (Fast version) | |
Write formatted data to string. | |
Remove singleton dimensions. | |
Set unique. | |
Display warning message; disable or enable warning messages. |
Joint (Scf,Hd) CDF for linear waves with Ochi-Hubble spectra. | |
Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. |
001 function [f,Hrms,Vrms,fA,fB] = ohhspdf(Hd,Scf,Hm0,def,normalizedInput, ... 002 condon) 003 %OHHSPDF Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. 004 % 005 % CALL: f = ohhspdf(Hd,Scf,Hm0,def) 006 % 007 % f = pdf 008 % Hd = zero down crossing wave height 009 % Scf = crest front steepness 010 % Hm0 = significant wave height [m]. 011 % def = defines the parametrization of the spectral density (default 1) 012 % 1 : The most probable spectrum (default) 013 % 2,3,...11 : gives 95% Confidence spectra 014 % 015 % OHHSPDF approximates the joint distribution of (Scf, Hd) in time, 016 % i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, 017 % for a Gaussian process with a bimodal Ochi-Hubble spectral density 018 % (ohspec2). The empirical 019 % parameters of the model is fitted by least squares to simulated 020 % (Scf,Hd) data for 24 classes of Hm0. Between 50000 and 150000 021 % zero-downcrossing waves were simulated for each class of Hm0. 022 % OHHSPDF is restricted to the following range for Hm0: 023 % 0.5 < Hm0 [m] < 12 024 % The size of f is the common size of the input arguments, Hd, Scf and 025 % Hm0. 026 % 027 % Example: 028 % Hm0 = 6;Tp = 8;def= 2; 029 % h = linspace(0,4*Hm0/sqrt(2))'; 030 % v = linspace(0,6*1.25*Hm0/Tp^2)'; 031 % [V,H] = meshgrid(v,h); 032 % f = ohhspdf(H,V,Hm0,def); 033 % contourf(v,h,f) 034 % 035 % See also ohhpdf, thspdf 036 037 % Reference 038 % P. A. Brodtkorb (2004), 039 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 040 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 041 % Trondheim, Norway. 042 043 % History 044 % revised pab 09.09.2003 045 % By pab 06.02.2001 046 047 error(nargchk(3,6,nargin)) 048 if (nargin < 6|isempty(condon)), condon = 0;end 049 if (nargin < 5|isempty(normalizedInput)), normalizedInput = 0;end 050 if (nargin < 4|isempty(def)), def=1;end 051 052 multipleSeaStates = any(prod(size(Hm0))>1); 053 if multipleSeaStates 054 [errorcode, Scf,Hd,Hm0] = comnsize(Scf,Hd,Hm0); 055 else 056 [errorcode, Scf,Hd] = comnsize(Scf,Hd); 057 end 058 if errorcode > 0 059 error('Requires non-scalar arguments to match in size.'); 060 end 061 062 if any(Hm0>12| Hm0<=0.5) 063 disp('Warning: Hm0 is outside the valid range') 064 disp('The validity of the Hd distribution is questionable') 065 end 066 067 if def>11|def<1 068 Warning('DEF is outside the valid range') 069 def = mod(def-1,11)+1; 070 end 071 072 global OHHSPAR 073 if isempty(OHHSPAR) 074 OHHSPAR = load('ohhspar.mat'); 075 end 076 077 078 %w = linspace(0,100,16*1024+1).'; %original spacing 079 %ch = spec2char(ohspec2(w,[Hm0,def]),{'Tm02','eps2'}); 080 %Tm02 = ch(1); 081 %Vrms = 2*Hm0/Tm02; 082 Hm00 = OHHSPAR.Hm0; 083 084 if normalizedInput, 085 Hrms = 1; 086 Vrms = 1; 087 else 088 method = 'cubic'; %'spline' 089 Tm020 = OHHSPAR.Tm02; 090 Hrms = Hm0/sqrt(2); 091 Tm02 = interp1(Hm00,Tm020(:,def),Hm0,method); 092 Vrms = 1.25*Hm0./(Tm02.^2); % Erms 093 end 094 095 096 h = Hd./Hrms; 097 v = Scf./Vrms; 098 cSize = size(h); % common size 099 100 % Logarithm of Weibull distribution parameters as a function of Hm0 and h2 101 A11 = squeeze(OHHSPAR.A11s(:,def,:)); 102 B11 = squeeze(OHHSPAR.B11s(:,def,:)); 103 104 h2 = OHHSPAR.h2; %Hd/Hrms 105 [E2 H2] = meshgrid(Hm00,h2); 106 method = '*cubic'; %'spline' 107 if multipleSeaStates 108 h = h(:); 109 v = v(:); 110 Hm0 = Hm0(:); 111 A1 = zeros(length(h),1); 112 B1 = A1; 113 [Hm0u,ix,jx] = unique(Hm0); 114 numSeaStates = length(ix); 115 Nh2 = length(h2); 116 Hm0i = zeros(Nh2,1); 117 for iz=1:numSeaStates 118 k = find(jx==iz); 119 Hm0i(:) = Hm0u(iz); 120 A1(k) = exp(smooth(h2,interp2(E2,H2,A11.',Hm0i,h2,method),... 121 1,h(k),1)); 122 B1(k) = (smooth(h2,interp2(E2,H2,exp(B11).',Hm0i,h2,method),... 123 1,h(k),1)); 124 end 125 else 126 Hm0i = Hm0(ones(size(h2))); 127 A1 = exp(smooth(h2,interp2(E2,H2,A11.', ... 128 Hm0i,h2,method),1,h,1)); 129 B1 = (smooth(h2,interp2(E2,H2,exp(B11).', ... 130 Hm0i,h2,method),1,h,1)); 131 end 132 %fh = ohhpdf(h(:)/Hrms,Hm0,def,'time',1); 133 % Fh = fh.f; 134 % Waveheight distribution in time 135 % Generalized gamma distribution parameters as a function of Hm0 136 [A0 B0 C0] = ohhgparfun(Hm0,def,'time'); 137 138 139 switch condon, 140 case 0, % regular pdf is returned 141 f = wggampdf(h,A0,B0,C0).*wweibpdf(v,A1,B1); 142 case 1, %pdf conditioned on x1 ie. p(x2|x1) 143 f = wweibpdf(v,A1,B1); 144 case 3, % secret option used by XXstat: returns x2*p(x2|x1) 145 f = v.*wweibpdf(v,A1,B1); 146 case 4, % secret option used by XXstat: returns x2.^2*p(x2|x1) 147 f = v.^2.*wweibpdf(v,A1,B1); 148 case 5, % p(h)*P(V|h) is returned special case used by thscdf 149 f = wggampdf(h,A0,B0,C0).*wweibcdf(v,A1,B1); 150 case 6, % P(V|h) is returned special case used by thscdf 151 f = wweibcdf(v,A1,B1); 152 case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf 153 f = wggampdf(h,A0,B0,C0).*(1-wweibcdf(v,A1,B1)); 154 otherwise error('unknown option') 155 end 156 if multipleSeaStates 157 f = reshape(f,cSize); 158 end 159 160 if condon~=6 161 f = f./Hrms./Vrms; 162 end 163 f(find(isnan(f)|isinf(f) ))=0; 164 if any(size(f)~=cSize) 165 disp('Wrong size') 166 end 167 if nargout>3, 168 fA = createpdf(2); 169 fA.f = A11; 170 fA.x = {Hm00, h2}; 171 fA.labx = {'Hm0', 'h'}; 172 fA.note = sprintf('%s \n',... 173 'The conditonal Weibull distribution Parameter',... 174 'A(h)/Hrms for wave heigth as a function of',... 175 'h=Hd/Hrms and Hm0 for the ohspec2 spectrum, ',... 176 sprintf('given def = %d.',def)); 177 178 ra = range(A11(:)); 179 st = round(min(A11(:))*100)/100; 180 en = max(A11(:)); 181 fA.cl = st:ra/20:en; 182 end 183 if nargout>4, 184 fB = createpdf(2); 185 fB.f = B11; 186 fB.x = {Hm00,h2}; 187 fB.labx = {'Hm0','h'}; 188 fB.note = sprintf('%s \n',... 189 'The conditonal Weibull distribution Parameter',... 190 'B(h)/Hrms for wave heigth as a function of',... 191 'h=Hd/Hrms and Hm0 for the ohspec2 spectrum, ',... 192 sprintf('given def = %d.',def)); 193 ra = range(B11(:)); 194 st = round(min(B11(:))*100)/100; 195 en = max(B11(:)); 196 fB.cl = st:ra/20:en; 197 end 198 return 199 200 201 202 203 204 205 206 return % old call 207 Fh = wggampdf(h(:)/Hrms,A0,B0,C0); 208 209 210 211 %C1 = exp(smooth(h2,interp2(E2,H2,log(squeeze(C11(:,def,:))).',Hm0(ones(size(h2))),h2,method),1,h/Hrms,1)); 212 %figure(1), plot(h/Hrms, A1,h2,squeeze(exp(A11(10,def,:))),'r') 213 %figure(2), plot(h/Hrms,B1,h2,squeeze(exp(B11(10,def,:))),'r') 214 [V1 A1] = meshgrid(v/Vrms,A1); 215 [V1 B1] = meshgrid(v/Vrms,B1); 216 %[V1 C1] = meshgrid(v/Vrms,C1); 217 218 f = createpdf(2); 219 f.title = 'Joint distribution of (Hd,Scf) in time'; 220 if norm==0 221 f.f = repmat(Fh/Hrms,[1 length(v)]).*wweibpdf(V1,A1,B1)/Vrms; 222 f.x = {v,h}; 223 f.norm=0; 224 f.labx={'Scf', 'Hd [m]'}; 225 else 226 f.f = repmat(Fh,[1 length(v)]).*wweibpdf(V1,A1,B1); 227 f.x = {v/Vrms,h/Hrms}; 228 f.labx={'Scf', 'Hd'}; 229 f.norm = 1; 230 end 231 f.note = ['Ochi-Hubble Hm0=' num2str(Hm0) ' def = ' num2str(def)]; 232 if 1, 233 f.f(isinf(f.f)|isnan(f.f))=0; 234 % Some places numerical problems occur 235 %f.f(f.f>30)=0; % Needed to calculate the contour levels 236 [f.cl,f.pl] = qlevels(f.f);%,[10:20:90 95 99 99.9 99.99 99.999 99.9999]); 237 end 238
Comments or corrections to the WAFO group