001 function [f,Hrms,Vrms,fA,fB] = jhvnlpdf(Hd,Vcf,Hm0,Tp,gam,normalizedInput,condon)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041 error(nargchk(4,7,nargin))
042 if (nargin < 7|isempty(condon)), condon = 0; end
043 if (nargin < 6|isempty(normalizedInput)), normalizedInput = 0;end
044 if (nargin < 5|isempty(gam))
045 gam = getjonswappeakedness(Hm0,Tp);
046 end
047
048 multipleSeaStates = any(prod(size(Hm0))>1|...
049 prod(size(Tp)) >1|...
050 prod(size(gam))>1);
051 if multipleSeaStates
052 [errorcode, Vcf,Hd,Hm0,Tp,gam] = comnsize(Vcf,Hd,Hm0,Tp,gam);
053 else
054 [errorcode, Vcf,Hd] = comnsize(Vcf,Hd);
055 end
056 if errorcode > 0
057 error('Requires non-scalar arguments to match in size.');
058 end
059
060 displayWarning = 0;
061 if displayWarning
062 if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0)))
063 disp('Warning: Hm0,Tp is outside the JONSWAP range')
064 disp('The validity of the parameters returned are questionable')
065 end
066 end
067 k = find(gam<1);
068 if any(k)
069 if displayWarning,
070 warning('Peakedness parameter less than 1. Must be larger than 1.')
071 end
072 gam(k)=1;
073 end
074 k1 = find(7<gam);
075 if any(k1)
076 if displayWarning
077 warning('Peakedness parameter larger than 7. The pdf returned is questionable')
078 end
079 gam(k1) = 7;
080 end
081
082 global JHVNLPAR
083 if isempty(JHVNLPAR)
084 JHVNLPAR = load('jhvnlpar.mat');
085 end
086
087 A11 = JHVNLPAR.A11s;
088 B11 = JHVNLPAR.B11s;
089
090
091
092
093
094 c1 =[ 0.16183666835624 1.53691936441548 1.55852759524555];
095 c2 =[ 0.15659478203944 1.15736959972513 1];
096 Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam));
097
098
099
100
101
102
103
104 if normalizedInput
105 Hrms = 1;
106 Vrms = 1;
107 else
108 Hrms = Hm0/sqrt(2);
109 Vrms = 2*Hm0./Tm02;
110 end
111
112 v = Vcf./Vrms;
113 h = Hd./Hrms;
114 cSize = size(h);
115
116 if gam<1
117 disp('Warning: Peakedness parameter less than 1. Must be larger than 1.')
118 gam=1;
119 elseif gam>7
120 disp('Warning: Peakedness parameter larger than 7. The pdf returned is questionable')
121 end
122
123
124 h2 = JHVNLPAR.h2(:);
125 method ='*cubic';
126 Nh2 = length(h2);
127 if 1,
128
129 Tpp = JHVNLPAR.Tp;
130 Hm00 = JHVNLPAR.Hm0;
131 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2);
132 Nh2 = length(h2);
133 if multipleSeaStates
134 h = h(:);
135 v = v(:);
136 Tp = Tp(:);
137 Hm0 = Hm0(:);
138 gam = gam(:);
139 A1 = zeros(length(h),1);
140 B1 = A1;
141 [TpHm0,ix,jx] = unique([Tp,Hm0],'rows');
142 numSeaStates = length(ix);
143 Tpi = zeros(Nh2,1);
144 Hm0i = zeros(Nh2,1);
145 for iz=1:numSeaStates
146 k = find(jx==iz);
147 Tpi(:) = TpHm0(iz,1);
148 Hm0i(:) = TpHm0(iz,2);
149 A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),...
150 1,h(k),1));
151 B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),...
152 1,h(k),1));
153 end
154 else
155 Tpi = repmat(Tp,[Nh2,1]);
156 Hm0i = repmat(Hm0,[Nh2,1]);
157 A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1));
158 B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1));
159 end
160 else
161 e2 = JHVNLPAR.gam;
162
163 [E2 H2] = meshgrid(e2,h2);
164 if multipleSeaStates
165 h = h(:);
166 v = v(:);
167 Tp = Tp(:);
168 Hm0 = Hm0(:);
169 gam = gam(:);
170
171 A1 = zeros(length(h),1);
172 B1 = A1;
173 [gamu,ix,jx] = unique(gam);
174 numSeaStates = length(gamu);
175 gami = zeros(Nh2,1);
176 for iz=1:numSeaStates
177 k = find(jx==iz);
178 gami(:) = gamu(iz);
179 A1(k) = smooth(h2,interp2(E2,H2,A11,gami,h2,method),1,h(k),1);
180 B1(k) = smooth(h2,interp2(E2,H2,B11,gami,h2,method),1,h(k),1);
181 end
182 else
183 A1 = smooth(h2,interp2(E2,H2,A11,gam(ones(size(h2))),h2,method),1,h,1);
184 B1 = smooth(h2,interp2(E2,H2,B11,gam(ones(size(h2))),h2,method),1,h,1);
185 end
186 end
187 [A0,B0,C0] = jhnlwparfun(Hm0,Tp,gam);
188
189 switch condon,
190 case 0,
191 f = wtweibpdf(h,A0,B0,C0).*wweibpdf(v,A1,B1);
192 case 1,
193 f = wweibpdf(v,A1,B1);
194 case 3,
195 f = v.*wweibpdf(v,A1,B1);
196 case 4,
197 f = v.^2.*wweibpdf(v,A1,B1);
198 case 5,
199 f = wtweibpdf(h,A0,B0,C0).*wweibcdf(v,A1,B1);
200 case 6,
201 f = wweibcdf(v,A1,B1);
202 case 7,
203 f = wtweibpdf(h,A0,B0,C0).*(1-wweibcdf(v,A1,B1));
204 otherwise error('unknown option')
205 end
206 if multipleSeaStates
207 f = reshape(f,cSize);
208 end
209
210 if condon~=6
211 f = f./Hrms./Vrms;
212 end
213 f(find(isnan(f)|isinf(f) ))=0;
214 if any(size(f)~=cSize)
215 disp('Wrong size')
216 end
217
218 if nargout>3,
219 fA = createpdf(2);
220 fA.f = A11;
221 fA.x = {e2,h2};
222 fA.labx = {'Gamma', 'h'};
223 fA.note = ['The conditonal Weibull distribution Parameter A(h,gamma)/Hrms' ...
224 'for wave heigth as a function of h=Hd/Hrms and gamma for' ...
225 'the Jonswap spectrum'];
226
227 ra = range(A11(:));
228 st = round(min(A11(:))*100)/100;
229 en = max(A11(:));
230 fA.cl = st:ra/20:en;
231 end
232 if nargout>4,
233 fB = createpdf(2);
234 fB.f = B11;
235 fB.x = {e2,h2};
236 fB.labx = {'Gamma', 'h'};
237 fB.note = ['The conditonal Weibull distribution Parameter B(h,gamma)/Hrms' ...
238 'for wave heigth as a function of h=Hd/Hrms and gamma for' ...
239 'the Jonswap spectrum'];
240 ra = range(B11(:));
241 st = round(min(B11(:))*100)/100;
242 en = max(B11(:));
243 fB.cl = st:ra/20:en;
244 end
245 return
246
247
248
249
250
251
252