001 function [f,Hrms,Vrms] = jhsnlpdf(Hd,Scf,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
042
043
044 error(nargchk(3,7,nargin))
045
046
047 if (nargin < 7|isempty(condon)), condon = 0; end
048 if (nargin < 6|isempty(normalizedInput)), normalizedInput = 0;end
049
050 if (nargin < 4|isempty(Tp)), Tp = 8;end
051 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end
052 if (nargin < 5|isempty(gam))
053
054 end
055
056 gam = getjonswappeakedness(Hm0,Tp);
057 multipleSeaStates = any(prod(size(Hm0))>1|...
058 prod(size(Tp)) >1|...
059 prod(size(gam))>1);
060 if multipleSeaStates
061 [errorcode, Scf,Hd,Hm0,Tp,gam] = comnsize(Scf,Hd,Hm0,Tp,gam);
062 else
063 [errorcode, Scf,Hd] = comnsize(Scf,Hd);
064 end
065 if errorcode > 0
066 error('Requires non-scalar arguments to match in size.');
067 end
068 displayWarning = 0;
069 if displayWarning
070 if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0)))
071 disp('Warning: Hm0,Tp is outside the JONSWAP range')
072 disp('The validity of the parameters returned are questionable')
073 end
074 end
075 k = find(gam<1);
076 if any(k)
077 if displayWarning,
078 warning('Peakedness parameter less than 1. Must be larger than 1.')
079 end
080 gam(k)=1;
081 end
082 k1 = find(7<gam);
083 if any(k1)
084 if displayWarning
085 warning('Peakedness parameter larger than 7. The pdf returned is questionable')
086 end
087 gam(k1) = 7;
088 end
089
090 global JHSNLPAR
091 if isempty(JHSNLPAR)
092
093 JHSNLPAR = load('jhsnlpar21-Jul-2004.mat');
094 end
095
096 A11 = JHSNLPAR.A11s;
097 B11 = JHSNLPAR.B11s;
098 C1 = 1;
099
100
101
102 if normalizedInput,
103 Hrms = 1;
104 Vrms = 1;
105 else
106
107
108 c1 =[ 0.16183666835624 1.53691936441548 1.55852759524555];
109 c2 =[ 0.15659478203944 1.15736959972513 1];
110 Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam));
111 Hrms = Hm0/sqrt(2);
112 Vrms = 1.25*Hm0./(Tm02.^2);
113 end
114
115 h = Hd./Hrms;
116 v = Scf./Vrms;
117 cSize = size(h);
118
119 h2 = JHSNLPAR.h2(:);
120 Nh2 = length(h2);
121 method ='*cubic';
122 if 1,
123 Tpp = JHSNLPAR.Tp;
124 Hm00 = JHSNLPAR.Hm0;
125 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2);
126 Nh2 = length(h2);
127 if multipleSeaStates
128 h = h(:);
129 v = v(:);
130 Tp = Tp(:);
131 Hm0 = Hm0(:);
132 gam = gam(:);
133 A1 = zeros(length(h),1);
134 B1 = A1;
135 [TpHm0,ix,jx] = unique([Tp,Hm0],'rows');
136 numSeaStates = length(ix);
137 Tpi = zeros(Nh2,1);
138 Hm0i = zeros(Nh2,1);
139 for iz=1:numSeaStates
140 k = find(jx==iz);
141 Tpi(:) = TpHm0(iz,1);
142 Hm0i(:) = TpHm0(iz,2);
143 A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),...
144 1,h(k),1));
145 B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),...
146 1,h(k),1));
147 end
148 else
149 Tpi = repmat(Tp,[Nh2,1]);
150 Hm0i = repmat(Hm0,[Nh2,1]);
151 A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1));
152 B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1));
153 end
154 else
155 e2 = JHSNLPAR.gam(:);
156 h2 = JHSNLPAR.h2(:);
157 [E2 H2] = meshgrid(e2,h2);
158 if multipleSeaStates
159 h = h(:);
160 v = v(:);
161 Tp = Tp(:);
162 Hm0 = Hm0(:);
163 gam = gam(:);
164
165 A1 = zeros(length(h),1);
166 B1 = A1;
167 [gamu,ix,jx] = unique(gam);
168 numSeaStates = length(gamu);
169 gami = zeros(Nh2,1);
170 for iz=1:numSeaStates
171 k = find(jx==iz);
172 gami(:) = gamu(iz);
173 A1(k) = exp(smooth(h2,interp2(E2,H2,log(A11),gami,h2,method),1,h(k),1));
174 B1(k) = exp(smooth(h2,interp2(E2,H2,log(B11),gami,h2,method),1,h(k),1));
175 end
176 else
177 A1 = exp(smooth(h2,interp2(E2,H2,log(A11),...
178 gam(ones(size(h2))),h2,method),1,h,1));
179 B1 = exp(smooth(h2,interp2(E2,H2,log(B11),...
180 gam(ones(size(h2))),h2,method),1,h,1));
181 end
182 end
183
184
185 [A0, B0, C0] = jhnlwparfun(Hm0,Tp,gam,'time');
186 switch condon,
187 case 0,
188 f = wtweibpdf(h,A0,B0,C0).*wggampdf(v,A1,B1,C1);
189 case 1,
190 f = wggampdf(v,A1,B1,C1);
191 case 3,
192 f = v.*wggampdf(v,A1,B1,C1);
193 case 4,
194 f = v.^2.*wggampdf(v,A1,B1,C1);
195 case 5,
196 f = wtweibpdf(h,A0,B0,C0).*wggamcdf(v,A1,B1,C1);
197 case 6,
198 f = wggamcdf(v,A1,B1,C1);
199 case 7,
200 f = wtweibpdf(h,A0,B0,C0).*(1-wggamcdf(v,A1,B1,C1));
201 otherwise error('unknown option')
202 end
203 if multipleSeaStates
204 f = reshape(f,cSize);
205 end
206
207 if condon~=6
208 f = f./Hrms./Vrms;
209 end
210 f(find(isnan(f)|isinf(f) ))=0;
211 if any(size(f)~=cSize)
212 disp('Wrong size')
213 end
214
215 return
216
217
218
219
220