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