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