001 function [f,Hrms,Vrms,fA,fB] = thvpdf(Hd,Vcf,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 error(nargchk(3,6,nargin))
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, Vcf,Hd,Hm0,Tp] = comnsize(Vcf,Hd,Hm0,Tp);
055 else
056 [errorcode, Vcf,Hd] = comnsize(Vcf,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,Vcf) 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,Vcf) distribution is questionable')
070 end
071 end
072
073 global THVPAR
074 if isempty(THVPAR)
075 THVPAR = load('thvpar.mat');
076 end
077
078 A11 = THVPAR.A11s;
079 B11 = THVPAR.B11s;
080 e2 = THVPAR.e2;
081 h2 = THVPAR.h2;
082 Tpp = THVPAR.Tp;
083 Hm00 = THVPAR.Hm0;
084 Tm020 = THVPAR.Tm02;
085 eps20 = THVPAR.eps2;
086
087 [Tp1,Hs1] = meshgrid(Tpp,Hm00);
088 method = '*cubic';
089
090 eps2 = interp2(Tp1,Hs1,eps20,Tp,Hm0,method);
091
092 if normalizedInput,
093 Hrms = 1;
094 Vrms = 1;
095 else
096 Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method);
097
098
099
100
101
102 Hrms = Hm0/sqrt(2);
103 Vrms = 2*Hm0./Tm02;
104 end
105
106
107 if displayWarning
108 if eps2<0.4
109 disp('Warning: eps2 is less than 0.4. The pdf returned is questionable')
110 elseif eps2>1.3
111 disp('Warning: eps2 is larger than 1.3. The pdf returned is questionable')
112 end
113 end
114 h = Hd/Hrms;
115 v = Vcf/Vrms;
116 cSize = size(h);
117
118 [E2 H2] = meshgrid(e2,h2);
119 Nh2 = length(h2);
120 if multipleSeaStates
121 h = h(:);
122 v = v(:);
123 Tp = Tp(:);
124 Hm0 = Hm0(:);
125 eps2 = eps2(:);
126 A1 = zeros(length(h),1);
127 B1 = A1;
128 [eps2u,ix,jx] = unique(eps2);
129 numSeaStates = length(ix);
130 eps2i = zeros(Nh2,1);
131 for iz=1:numSeaStates
132 k = find(jx==iz);
133 eps2i(:) = eps2u(iz);
134 A1(k) = smooth(h2,interp2(E2,H2,A11,eps2i,h2,method),1,h(k),1);
135 B1(k) = smooth(h2,interp2(E2,H2,B11,eps2i,h2,method),1,h(k),1);
136 end
137 else
138 eps2i = repmat(eps2,[Nh2,1]);
139 A1 = smooth(h2,interp2(E2,H2,A11,eps2i,h2,method),1,h,1);
140 B1 = smooth(h2,interp2(E2,H2,B11,eps2i,h2,method),1,h,1);
141 end
142
143
144
145
146 [A0, B0, C0] = thwparfun(Hm0,Tp,'time');
147
148 switch condon,
149 case 0,
150 f = wtweibpdf(h,A0,B0,C0).*wweibpdf(v,A1,B1);
151 case 1,
152 f = wweibpdf(v,A1,B1);
153 case 3,
154 f = v.*wweibpdf(v,A1,B1);
155 case 4,
156 f = v.^2.*wweibpdf(v,A1,B1);
157 case 5,
158 f = wtweibpdf(h,A0,B0,C0).*wweibcdf(v,A1,B1);
159 case 6,
160 f = wweibcdf(v,A1,B1);
161 case 7,
162 f = wtweibpdf(h,A0,B0,C0).*(1-wweibcdf(v,A1,B1));
163 otherwise error('unknown option')
164 end
165 if multipleSeaStates
166 f = reshape(f,cSize);
167 end
168
169 if condon~=6
170 f = f./Hrms./Vrms;
171 end
172 f(find(isnan(f)|isinf(f) ))=0;
173 if any(size(f)~=cSize)
174 disp('Wrong size')
175 end
176
177 if nargout>3,
178 fA = createpdf(2);
179 fA.f = A11;
180 fA.x = {e2,h2};
181 fA.labx = {'eps2', 'h'};
182 fA.title = 'wweibpdf parameter A';
183 fA.note = ['The conditonal Weibull distribution Parameter A(h,eps2)/Hrms' ...
184 'for Vcf given h=Hd/Hrms and eps2 for' ...
185 'the Torsethaugen spectrum'];
186
187 ra = range(A11(:));
188 st = round(min(A11(:))*100)/100;
189 en = max(A11(:));
190 fA.cl = st:ra/20:en;
191 end
192 if nargout>4,
193 fB = createpdf(2);
194 fB.f = B11;
195 fB.x = {e2,h2};
196 fB.labx = {'eps2', 'h'};
197 fB.title = 'wweibpdf parameter B';
198 fB.note = ['The conditonal Weibull distribution Parameter B(h,eps2)/Hrms' ...
199 'for Vcf given h=Hd/Hrms and eps2 for' ...
200 'the Torsethaugen spectrum'];
201 ra = range(B11(:));
202 st = round(min(B11(:))*100)/100;
203 en = max(B11(:));
204 fB.cl = st:ra/20:en;
205 end
206 return
207
208
209
210
211
212
213