001 function [f,Hrms,Vrms,fA,fB] = thsspdf(Hd,Scf,Hm0,Tp,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
042
043
044
045
046
047
048
049
050
051
052
053
054 error(nargchk(3,6,nargin))
055 if (nargin < 6|isempty(condon)), condon = 0;end
056 if (nargin < 5|isempty(normalizedInput)), normalizedInput = 0;end
057 if (nargin < 4|isempty(Tp)), Tp = 8;end
058 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end
059
060
061 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1);
062 if multipleSeaStates
063 [errorcode, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp);
064 else
065 [errorcode, Scf,Hd] = comnsize(Scf,Hd);
066 end
067 if errorcode > 0,
068 error('Requires non-scalar arguments to match in size.');
069 end
070 displayWarning = 0;
071 if displayWarning,
072 if any(Hm0>11| Hm0>(Tp-2)*12/11)
073 disp('Warning: Hm0 is outside the valid range')
074 disp('The validity of the Joint (Hd,Scf) distribution in space is questionable')
075 end
076 if any(Tp>20|Tp<3 )
077 disp('Warning: Tp is outside the valid range')
078 disp('The validity of the Joint (Hd,Scf) distribution in space is questionable')
079 end
080 end
081 useWeibull = 1;
082 if useWeibull
083 global THSSPARW
084 if isempty(THSSPARW)
085 THSSPARW = load('thsspar27-Jul-2004.mat');
086 end
087
088
089 A11 = THSSPARW.A11s;
090 B11 = THSSPARW.B11s;
091
092
093 A00 = THSSPARW.A00s;
094 B00 = THSSPARW.B00s;
095 C00 = THSSPARW.C00s;
096
097 Tpp = THSSPARW.Tp;
098 Hm00 = THSSPARW.Hm0;
099 h2 = THSSPARW.h2(:);
100 Tm020 = THSSPARW.Tm02;
101 else
102 global THSSPARG
103 if isempty(THSSPARG)
104 THSSPARG = load('thsspar.mat');
105 end
106
107
108 A11 = THSSPARG.A11s;
109 B11 = THSSPARG.B11s;
110
111
112 A00 = THSSPARG.A00s;
113 B00 = THSSPARG.B00s;
114 C00 = THSSPARG.C00s;
115
116 Tpp = THSSPARG.Tp;
117 Hm00 = THSSPARG.Hm0;
118 h2 = THSSPARG.h2;
119 Tm020 = THSSPARG.Tm02;
120 end
121 if normalizedInput,
122 Hrms = 1;
123 Vrms = 1;
124 else
125 method = 'cubic';
126 [Tp1,Hs1] = meshgrid(Tpp,Hm00);
127 Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method);
128
129
130
131
132 Hrms = Hm0/sqrt(2);
133 Vrms = 2*Hm0./Tm02;
134 end
135
136
137
138
139
140
141
142
143
144
145
146 h = Hd./Hrms;
147 v = Scf./Vrms;
148 cSize = size(h);
149
150 method = '*cubic';
151
152
153 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2);
154 Nh2 = length(h2);
155
156 if multipleSeaStates
157 h = h(:);
158 v = v(:);
159 Tp = Tp(:);
160 Hm0 = Hm0(:);
161 A1 = zeros(length(h),1);
162 B1 = A1;
163 [TpHm0,ix,jx] = unique([Tp,Hm0],'rows');
164 numSeaStates = length(ix);
165 Tpi = zeros(Nh2,1);
166 Hm0i = zeros(Nh2,1);
167 for iz=1:numSeaStates
168 k = find(jx==iz);
169 Tpi(:) = TpHm0(iz,1);
170 Hm0i(:) = TpHm0(iz,2);
171 A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),...
172 1,h(k),1));
173 B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),...
174 1,h(k),1));
175 end
176 else
177 Tpi = repmat(Tp,[Nh2,1]);
178 Hm0i = repmat(Hm0,[Nh2,1]);
179 A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1));
180 B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1));
181 end
182
183 [E1, H1] = meshgrid(Tpp,Hm00);
184 A0 = interp2(E1,H1,A00,Tp,Hm0,method);
185 B0 = interp2(E1,H1,B00,Tp,Hm0,method);
186 C0 = interp2(E1,H1,C00,Tp,Hm0,method);
187
188 if useWeibull
189 switch condon,
190 case 0,
191 f = wtweibpdf(h,A0,B0,C0).*wgampdf(v,A1,B1);
192 case 1,
193 f = wgampdf(v,A1,B1);
194 case 3,
195 f = v.*wgampdf(v,A1,B1);
196 case 4,
197 f = v.^2.*wgampdf(v,A1,B1);
198 case 5,
199 f = wtweibpdf(h,A0,B0,C0).*wgamcdf(v,A1,B1);
200 case 6,
201 f = wgamcdf(v,A1,B1);
202 case 7,
203 f = wtweibpdf(h,A0,B0,C0).*(1-wgamcdf(v,A1,B1));
204 otherwise error('unknown option')
205 end
206 else
207
208 switch condon,
209 case 0,
210 f = wggampdf(h,A0,B0,C0).*wgampdf(v,A1,B1);
211 case 1,
212 f = wgampdf(v,A1,B1);
213 case 3,
214 f = v.*wgampdf(v,A1,B1);
215 case 4,
216 f = v.^2.*wgampdf(v,A1,B1);
217 case 5,
218 f = wggampdf(h,A0,B0,C0).*wgamcdf(v,A1,B1);
219 case 6,
220 f = wgamcdf(v,A1,B1);
221 case 7,
222 f = wggampdf(h,A0,B0,C0).*(1-wgamcdf(v,A1,B1));
223 otherwise error('unknown option')
224 end
225 end
226 if multipleSeaStates
227 f = reshape(f,cSize);
228 end
229
230 if condon~=6
231 f = f./Hrms./Vrms;
232 end
233 f(find(isnan(f)|isinf(f) ))=0;
234 if any(size(f)~=cSize)
235 disp('Wrong size')
236 end
237
238 if nargout>3,
239 fA = createpdf(2);
240 fA.x = {Tpp,Hm00};
241 fA.labx = {'Tp', 'Hm0'};
242 fA(3) = fA(1);
243 fA(2) = fA(1);
244
245 fA(1).f = A00;
246 fA(2).f = B00;
247 fA(3).f = C00;
248
249 fA(1).title = 'wggampdf parameter A';
250 fA(2).title = 'wggampdf parameter B';
251 fA(3).title = 'wggampdf parameter C';
252
253 txt1 = ['The Wggampdf distribution Parameter '];
254 txt2=[' for wave heigth in space as a function of Tp and Hm0 for' ...
255 ' the Torsethaugen spectrum'];
256 fA(1).note =[txt1 'A' txt2];
257 fA(2).note =[txt1 'B' txt2];
258 fA(3).note =[txt1 'C' txt2];
259
260 tmp= [A00(:) B00(:) C00(:)];
261 ra = range(tmp);
262 st = round(min(tmp)*100)/100;
263 en = max(tmp);
264 for ix=1:3
265 fA(ix).cl = st(ix):ra(ix)/20:en(ix);
266 end
267 end
268 if nargout>4,
269 fB = createpdf(3);
270 fB.x = {Tpp,Hm00,h2};
271 fB.labx = {'Tp','Hm0', 'h'};
272 fB(2) = fB(1);
273
274 fB(1).f = A11;
275 fB(2).f = B11;
276
277 txt11 = 'The conditonal Wgampdf distribution Parameter ';
278 txt22 = [' for Scf given h=Hd/Hrms in space as function of Tp' ...
279 ' and Hm0 for the Torsethaugen spectrum'];
280 fB(1).title = 'wgampdf parameter A';
281 fB(2).title = 'wgampdf parameter B';
282
283 fB(1).note = [txt11,'A',txt22];
284 fB(2).note = [txt11,'B',txt22];
285
286
287
288 tmp= [A11(:) B11(:)];
289 ra = range(tmp);
290 st = round(min(tmp)*100)/100;
291 en = max(tmp);
292 for ix=1:2
293 fB(ix).cl = st(ix):ra(ix)/20:en(ix);
294 end
295 end
296 return
297
298
299
300
301
302
303