SPREADING Directional spreading functions CALL: D = spreading(th,type,th0,data,w,def); D = spreading(Nt,type,th0,data,w,def); D = spreading function, struct reminding of spectrum struct th = vector of direction angles, (default linspace(-pi,pi,Nt)) Nt = scalar defining the length of th. (default 101) type = type of spreading function, see options below (default 'cos2s') th0 = vector or a scalar defining average direction at every frequency. (length 1 or length == length(w)) (default 0) data = vector of spreading parameters, see options below. (default [15 15 0.52 5 -2.5 0 1 inf]) w = frequency vector, if frequency dependent spreading (default linspace(0,3,257)) def = 0 No frequency dependence of spreading function 1 Mitsuyasu et al., Hasselman et al (JONSWAP experiment) frequency dependent parametrization (default) The different types of spreading functions implemented are: type = 'cos2s' : cos-2s spreading N(S)*[cos((th-th0)/2)]^(2*S) (0 < S) 'box' : Box-car spreading N(A)*I( -A < th-th0 < A) (0 < A < pi) 'mises' : von Mises spreading N(K)*exp(K*cos(th-th0)) (0 < K) 'poisson': Poisson spreading N(X)/(1-2*X*cos(th-th0)+X^2) (0 < X < 1) 'sech2' : sech-2 spreading N(B)*0.5*B*sech(B*(th-th0))^2 (0 < B) 'wnormal': Wrapped Normal [1 + 2*sum exp(-(n*D)^2/2)*cos(n*(th-th0))]/(2*pi) (0 < D) (N(.) = normalization factor) (the first letter is enough for unique identification) Here the S-parameter of the COS-2S spreading function is used as a measure of spread. All the parameters of the other distributions are related to this S-parameter throug the first Fourier coefficient, R1, of the directional distribution as follows: R1 = S/(S+1) or S = R1/(1-R1). where Box-car spreading : R1 = sin(A)/A Von Mises spreading: R1 = besseli(1,K)/besseli(0,K), Poisson spreading : R1 = X sech-2 spreading : R1 = pi/(2*B*sinh(pi/(2*B)) Wrapped Normal : R1 = exp(-D^2/2) Use of data vector: frequency dependent spreading: data = [spa spb wc ma mb wlim], sp = maximum spread wc = cut over frequency (usually the peak frequency, wp or fp) m = shape parameter defining the frequency dependent spreading parameter, S=S(w), where S(w) = sp *(w/wc)^m, with sp=spa, m=ma for wlim(1) <= w/wc < wlim(2) sp=spb, m=mb for wlim(2) <= w/wc < wlim(3) = 0 otherwise frequency independent spreading: data = S, default value: 15 which corresponds to 'cos2s' : S=15 'box' : A=0.62 'sech2' : B=0.89 'mises' : K=8.3 'poisson': X=0.94 'wnormal': D=0.36 The 'cos2s' is the most frequently used spreading in engineering practice. Apart from the current meter/pressure cell data in WADIC all instruments seem to support the 'cos2s' distribution for heavier sea states, (Krogstad and Barstow, 1999). For medium sea states a spreading function between 'cos2s' and 'poisson' seem appropriate, while 'poisson' seems appropriate for swell. For the 'cos2s' Mitsuyasu et al. parameterized SPa = SPb = 11.5*(U10/Cp) where Cp = g/wp is the deep water phase speed at wp and U10 the wind speed at reference height 10m. Hasselman et al. (1980) parameterized mb = -2.33-1.45*(U10/Cp-1.17). Mitsuyasu et al. (1975) showed that SP for wind waves varies from 5 to 30 being a function of dimensionless wind speed. However, Goda and Suzuki (1975) proposed SP = 10 for wind waves, SP = 25 for swell with short decay distance and SP = 75 for long decay distance. Compared to experiments Krogstad et al. (1998) found that ma = 5 +/- eps and that -1< mb < -3.5. Values given in the litterature: [spa spb wc ma mb wlim(1:3) ] (Mitsuyasu: spa == spb) (cos-2s) [15 15 0.52 5 -2.5 0 1 inf] (Hasselman: spa ~= spb) (cos-2s) [6.97 9.77 0.52 4.06 -2.52 0 1 inf] NOTE: - by specifying NaN's in the data vector default values will be used. - if length(data) is shorter than the parameters needed then the default values are used Example:% Set spa = 10, wc = 0.43 and spb, ma, mb, wlim to their % default values, respectively: data = [10, nan, .43]; D = spreading(51,'cos2s',0,data) % Frequency dependent direction th0 = linspace(0,pi/2,257)'; D = spreading(51,'cos2s',th0,data) See also mkdspec, wspecplot, spec2spec
Sin(pi*x)/(pi*x) function. | |
Modified Bessel function of the first kind. | |
Display message and abort function. | |
Discrete Fourier transform. | |
Scalar nonlinear zero finding. | |
Gamma function. | |
1-D interpolation (table lookup) | |
Linearly spaced vector. | |
Convert string to lowercase. | |
Create/alter OPTIM OPTIONS structure. | |
Write formatted data to string. | |
Convert string matrix to numeric array. | |
Compare first N characters of strings. | |
Create or convert to structure array. | |
MATLAB version number. | |
Display warning message; disable or enable warning messages. |
% CHAPTER1 demonstrates some applications of WAFO | |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
Loads a precreated spectrum of chosen type | |
Make a directional spectrum | |
Directional spectra using Thorsethaugen and cos2s spreading |
0001 function D = spreading(th,type,th0,data,w,def) 0002 %SPREADING Directional spreading functions 0003 % 0004 % CALL: D = spreading(th,type,th0,data,w,def); 0005 % D = spreading(Nt,type,th0,data,w,def); 0006 % 0007 % D = spreading function, struct reminding of spectrum struct 0008 % th = vector of direction angles, (default linspace(-pi,pi,Nt)) 0009 % Nt = scalar defining the length of th. (default 101) 0010 % type = type of spreading function, see options below (default 'cos2s') 0011 % th0 = vector or a scalar defining average direction at every frequency. 0012 % (length 1 or length == length(w)) (default 0) 0013 % data = vector of spreading parameters, see options below. 0014 % (default [15 15 0.52 5 -2.5 0 1 inf]) 0015 % w = frequency vector, if frequency dependent spreading 0016 % (default linspace(0,3,257)) 0017 % def = 0 No frequency dependence of spreading function 0018 % 1 Mitsuyasu et al., Hasselman et al (JONSWAP experiment) 0019 % frequency dependent parametrization (default) 0020 % 0021 % The different types of spreading functions implemented are: 0022 % type = 'cos2s' : cos-2s spreading N(S)*[cos((th-th0)/2)]^(2*S) (0 < S) 0023 % 'box' : Box-car spreading N(A)*I( -A < th-th0 < A) (0 < A < pi) 0024 % 'mises' : von Mises spreading N(K)*exp(K*cos(th-th0)) (0 < K) 0025 % 'poisson': Poisson spreading N(X)/(1-2*X*cos(th-th0)+X^2) (0 < X < 1) 0026 % 'sech2' : sech-2 spreading N(B)*0.5*B*sech(B*(th-th0))^2 (0 < B) 0027 % 'wnormal': Wrapped Normal 0028 % [1 + 2*sum exp(-(n*D)^2/2)*cos(n*(th-th0))]/(2*pi) (0 < D) 0029 % (N(.) = normalization factor) 0030 % (the first letter is enough for unique identification) 0031 % 0032 % Here the S-parameter of the COS-2S spreading function is used as a 0033 % measure of spread. All the parameters of the other distributions are 0034 % related to this S-parameter throug the first Fourier coefficient, R1, of the 0035 % directional distribution as follows: 0036 % R1 = S/(S+1) or S = R1/(1-R1). 0037 % where 0038 % Box-car spreading : R1 = sin(A)/A 0039 % Von Mises spreading: R1 = besseli(1,K)/besseli(0,K), 0040 % Poisson spreading : R1 = X 0041 % sech-2 spreading : R1 = pi/(2*B*sinh(pi/(2*B)) 0042 % Wrapped Normal : R1 = exp(-D^2/2) 0043 % 0044 % Use of data vector: 0045 % frequency dependent spreading: 0046 % data = [spa spb wc ma mb wlim], 0047 % sp = maximum spread 0048 % wc = cut over frequency (usually the peak frequency, wp or fp) 0049 % m = shape parameter defining the frequency dependent 0050 % spreading parameter, S=S(w), where 0051 % S(w) = sp *(w/wc)^m, with sp=spa, m=ma for wlim(1) <= w/wc < wlim(2) 0052 % sp=spb, m=mb for wlim(2) <= w/wc < wlim(3) 0053 % = 0 otherwise 0054 % frequency independent spreading: 0055 % data = S, default value: 15 which corresponds to 0056 % 'cos2s' : S=15 0057 % 'box' : A=0.62 0058 % 'sech2' : B=0.89 0059 % 'mises' : K=8.3 0060 % 'poisson': X=0.94 0061 % 'wnormal': D=0.36 0062 % 0063 % The 'cos2s' is the most frequently used spreading in engineering practice. 0064 % Apart from the current meter/pressure cell data in WADIC all 0065 % instruments seem to support the 'cos2s' distribution for heavier sea 0066 % states, (Krogstad and Barstow, 1999). For medium sea states 0067 % a spreading function between 'cos2s' and 'poisson' seem appropriate, 0068 % while 'poisson' seems appropriate for swell. 0069 % For the 'cos2s' Mitsuyasu et al. parameterized SPa = SPb = 0070 % 11.5*(U10/Cp) where Cp = g/wp is the deep water phase speed at wp and 0071 % U10 the wind speed at reference height 10m. Hasselman et al. (1980) 0072 % parameterized mb = -2.33-1.45*(U10/Cp-1.17). 0073 % Mitsuyasu et al. (1975) showed that SP for wind waves varies from 0074 % 5 to 30 being a function of dimensionless wind speed. 0075 % However, Goda and Suzuki (1975) proposed SP = 10 for wind waves, SP = 25 0076 % for swell with short decay distance and SP = 75 for long decay distance. 0077 % Compared to experiments Krogstad et al. (1998) found that ma = 5 +/- eps and 0078 % that -1< mb < -3.5. 0079 % Values given in the litterature: [spa spb wc ma mb wlim(1:3) ] 0080 % (Mitsuyasu: spa == spb) (cos-2s) [15 15 0.52 5 -2.5 0 1 inf] 0081 % (Hasselman: spa ~= spb) (cos-2s) [6.97 9.77 0.52 4.06 -2.52 0 1 inf] 0082 % 0083 % NOTE: - by specifying NaN's in the data vector default values will be used. 0084 % - if length(data) is shorter than the parameters needed then the 0085 % default values are used 0086 % 0087 % Example:% Set spa = 10, wc = 0.43 and spb, ma, mb, wlim to their 0088 % % default values, respectively: 0089 % 0090 % data = [10, nan, .43]; 0091 % D = spreading(51,'cos2s',0,data) 0092 % % Frequency dependent direction 0093 % th0 = linspace(0,pi/2,257)'; 0094 % D = spreading(51,'cos2s',th0,data) 0095 % 0096 % See also mkdspec, wspecplot, spec2spec 0097 0098 % References 0099 % Krogstad, H.E. and Barstow, S.F. (1999) 0100 % "Directional Distributions in Ocean Wave Spectra" 0101 % Proceedings of the 9th ISOPE Conference, Vol III, pp. 79-86 0102 % 0103 % Goda, Y. (1999) 0104 % "Numerical simulation of ocean waves for statistical analysis" 0105 % Marine Tech. Soc. Journal, Vol. 33, No. 3, pp 5--14 0106 % 0107 % Banner, M.L. (1990) 0108 % "Equilibrium spectra of wind waves." 0109 % J. Phys. Ocean, Vol 20, pp 966--984 0110 % 0111 % Donelan M.A., Hamilton J, Hui W.H. (1985) 0112 % "Directional spectra of wind generated waves." 0113 % Phil. Trans. Royal Soc. London, Vol A315, pp 387--407 0114 % 0115 % Hasselmann D, Dunckel M, Ewing JA (1980) 0116 % "Directional spectra observed during JONSWAP." 0117 % J. Phys. Ocean, Vol.10, pp 1264--1280 0118 % 0119 % Mitsuyasu, H, et al. (1975) 0120 % "Observation of the directional spectrum of ocean waves using a 0121 % coverleaf buoy." 0122 % J. Physical Oceanography, Vol.5, No.4, pp 750--760 0123 0124 % Some of this might be included in help header: 0125 % cos-2s: 0126 % NB! The generally strong frequency dependence in directional spread 0127 % makes it questionable to run load tests of ships and structures with a 0128 % directional spread independent of frequency (Krogstad and Barstow, 1999). 0129 0130 % Parameterization of B 0131 % def = 2 Donelan et al freq. parametrization for 'sech2' 0132 % def = 3 Banner freq. parametrization for 'sech2' 0133 % (spa ~= spb) (sech-2) [2.61 2.28 0.52 1.3 -1.3 0.56 0.95 1.6] 0134 0135 0136 % Tested on: Matlab 5.3 0137 % History: 0138 % revisd pab 17.06.2001 0139 % - added wrapped normal spreading 0140 % revised pab 6 April 2001 0141 % - added fcof2par 0142 % - Fixed the normalization of sech2 spreading 0143 % revised by PAB and IR 1 April 2001: Introducing the azymuth as a 0144 % standard parameter in order to avoid rotations of the directions 0145 % theta. The x-axis is always pointing into the principal direction 0146 % as defined in the spreading function D(omega,theta). The actual 0147 % principal direction is defined by means of field D.phi. 0148 % revised es 06.06.2000, commented away: if ((ma==0) & (mb==0)), ..., 0149 % hoping that the check is unnecessary 0150 % revised pab 13.06.2000 0151 % - fixed a serious bug: made sure -pi<= th-th0 <=pi 0152 % revised pab 16.02.2000 0153 % -fixed default value for Hasselman parametrization 0154 % revised pab 02.02.2000 0155 % - Nt or th may be specified + check on th 0156 % - added frequency dependence for sech-2 0157 % - th0 as separate input 0158 % - updated header info 0159 % - changed check for nargins 0160 % - added the possibility of nan's in data vector 0161 % Revised by jr 2000.01.25 0162 % - changed check of nargins 0163 % - frequency dependence only for cos-2s 0164 % - updated information 0165 % By es, jr 1999.11.25 0166 0167 error(nargchk(0,6,nargin)) 0168 % Default values 0169 %~~~~~~~~~~~~~~~ 0170 Nt=101;Nw = 257; 0171 if nargin<1|length(th)<2 0172 if (nargin>0) & (length(th)==1), Nt=abs(th);end 0173 th = linspace(-pi,pi,Nt).'; 0174 elseif abs(th(end)-th(1))<2*pi-eps | abs(th(end)-th(1))>2*pi+eps 0175 error('theta must cover all angles -pi -> pi') 0176 else 0177 Nt = length(th); 0178 end 0179 if Nt<40, warning('Number of angles is less than 40. Spreading too sparsely sampled!'), end 0180 if nargin<2|isempty(type), type = 'cos2s'; end 0181 if nargin<3|isempty(th0), th0 = 0; end 0182 if nargin<5|isempty(w), w = linspace(0,3,Nw).'; end 0183 if nargin<6|isempty(def), def = 1; end 0184 0185 % Determine the default values 0186 0187 data2=[15 15 0.52 5 -2.5 0 1 inf]; 0188 spb=[];wc =[]; ma=[]; mb=[];wlim=[]; 0189 nd2 = length(data2); 0190 if (nargin>3) & ~isempty(data), 0191 nd = length(data); 0192 ind = find(~isnan(data(1:min(nd,nd2)))); 0193 if any(ind) % replace default values with those from input data 0194 data2(ind)=data(ind); 0195 end 0196 %if ((def==1) & strncmp(type,'cos2s',1) & ((nd<2)|isempty(ind)|~any(ind==2))), 0197 % data2(2)=data2(1); 0198 %end 0199 end 0200 0201 if (nd2>0), spa = data2(1); end 0202 if (nd2>1), spb = data2(2); end 0203 if def~=0 0204 if (nd2>2), wc = data2(3); end 0205 if (nd2>3), ma = data2(4); end 0206 if (nd2>4), mb = data2(5); end 0207 if (nd2>5), wlim = data2(6:5+3); end 0208 end 0209 %th = mod(th(:)+pi,2*pi)-pi; 0210 % Make sure theta is from -pi to pi 0211 th = th(:); 0212 phi0 = th(1)+pi; 0213 th = th-phi0; 0214 th(end) = pi; 0215 phi0 = mod(phi0,2*pi); 0216 0217 th0 = th0(:).'; 0218 Nt0 = length(th0); 0219 w = w(:); 0220 Nw = length(w); 0221 0222 if Nt0==Nw, 0223 % frequency dependent spreading and/or 0224 % frequency dependent direction 0225 TH = mod(th(:,ones(1,Nw))-th0(ones(Nt,1),:)+pi,2*pi)-pi; % make sure -pi<=TH<pi 0226 if def==0, def =1;end 0227 elseif Nt0~=1, 0228 error('The length of th0 must equal to 1 or the length of w') 0229 else 0230 % If ma==0 and mb==0 then frequency independent spreading is wanted 0231 %if ((nd2>4) & (ma==0) & (mb==0)), def=0;spb=[];wc =[]; ma=[]; mb=[];wlim=[];end 0232 TH = mod(th-th0+pi,2*pi)-pi; % make sure -pi<=TH<pi 0233 if def~=0, % frequency dependent spreading 0234 TH = TH(:,ones(1,Nw)); 0235 end 0236 end 0237 0238 % Createspreading 0239 %~~~~~~~~~~~~~~~~ 0240 if def>0 0241 D = struct('S',[],'w',w,'theta',th,'type','dir','phi',phi0,'note',[]); % frequency dependent 0242 if ~isempty(wc),wn = w./wc; end % normalized frequency 0243 else 0244 D = struct('S',[],'theta',th,'type','dir','phi',phi0,'note',[]); % frequency independent 0245 end 0246 0247 0248 D.th0 = th0; 0249 D.data = [spa spb wc ma mb wlim]; % save the spreading parameters used 0250 if def==0|isempty(wc), 0251 % no frequency dependent spreading, but possible frequency dependent 0252 % direction 0253 s = spa; 0254 else 0255 k = find(wlim(3)<wn); 0256 % Mitsuyasu et. al and Hasselman et. al parametrization of 0257 % frequency dependent spreading 0258 s=[ zeros(length(find(wn<=wlim(1))),1) ; ... 0259 spa*(wn((wlim(1)<wn) & (wn<=wlim(2)))).^ma ;... 0260 spb*(wn((wlim(2)<wn) & (wn<=wlim(3)))).^mb ;... 0261 zeros(length(k),1) ].'; 0262 if def==2 & any(k), 0263 % Donelan et. al. parametrization for B in SECH-2 0264 s(k) = s(max(find(s~=0))); 0265 end 0266 if def==3 & any(k), 0267 % Banner parametrization for B in SECH-2 0268 s(k) = 10.^(-0.4+0.8393*exp(-0.567*log(wn(k).^2))); 0269 end 0270 end 0271 0272 if any(s<0), error('Spreading: SP value must be larger than 0'),end 0273 0274 if strncmp(type,'cos2s',1), 0275 if isempty(wc),S=s;else,S = s(ones(Nt,1),:);end 0276 D.S = gamma(S+1)/2/sqrt(pi)./gamma(S+1/2).*cos(TH/2).^(2*S); 0277 D.note='Spreading: cos2s'; 0278 return 0279 end 0280 0281 r1 = abs(s./(s+1)); % First Fourier coefficient of the directional spreading function. 0282 %plot(s,r1) 0283 if strncmp(type,'poisson',1), 0284 if isempty(wc),X = r1; else, X = r1(ones(Nt,1),:);end 0285 if any(X>=1), error('POISSON spreading: X value must be less than 1'),end 0286 D.S = (1-X.^2)./(1-(2*cos(TH)-X).*X)/(2*pi); 0287 D.note = 'Spreading: Poisson'; 0288 return 0289 end 0290 if strncmp(type,'wnormal',1), 0291 D1 = repmat(inf,size(r1)); 0292 ind = find(r1); % avoid log of 0 0293 D1(ind) = sqrt(-2*log(r1(ind))); 0294 %D1 = sqrt(-2*log(r1)); 0295 if any(D1<0|~isreal(D1)), 0296 error('WRAPPED NORMAL spreading: D value must be greater than 0.'), 0297 end 0298 D1 = D1.^2/2; 0299 ix = (1:Nt-1).'; 0300 ix2 = ix.^2; 0301 if ~isempty(wc), 0302 D1 = D1(ones(Nt-1,1),:); 0303 end 0304 Nd2 = size(D1,2); 0305 Fcof = [ ones(1,Nd2)/2;... 0306 exp(-ix2(:,ones(1,Nd2)).*D1) ]/pi; 0307 Pcor = [ ones(1,Nd2 ); exp(sqrt(-1)*ix*TH(1,:))]; % correction term to get 0308 % the correct integration limits 0309 Fcof = Fcof(1:Nt,:).*conj(Pcor); 0310 D.S = real(fft(Fcof)); 0311 D.S(D.S<0)=0; 0312 D.note = 'Spreading: Wrapped Normal'; 0313 return 0314 end 0315 % Find distribution parameter from first Fourier coefficient. 0316 par = fcof2par(r1,type); 0317 0318 0319 if strncmp(type,'sech2',1), 0320 NB = tanh(pi*par); % Normalization factor. 0321 NB(NB==0) = 1; % Avoid division by zero 0322 if isempty(wc),B=par; else,B = par(ones(Nt,1),:); NB = NB(ones(Nt,1),:); end 0323 D.S = 0.5*B.*sech(B.*TH).^2./NB; 0324 D.note='Spreading: sech2'; 0325 return 0326 end 0327 0328 if strncmp(type,'mises',1), 0329 if isempty(wc),K = par; else,K = par(ones(Nt,1),:);end 0330 D.S=1./(2*pi*besseli(0,K,1)).*exp(K.*(cos(TH)-1)); 0331 D.note='Spreading: von Mises'; 0332 return 0333 end 0334 0335 if strncmp(type,'box',1), 0336 if any(par>pi), error('BOX-CAR spreading: The A value must be less than pi'),end 0337 if isempty(wc),A=par; else, A = par(ones(Nt,1),:);end 0338 D.S=1./(2.*A).*(-A<=TH & TH<=A); 0339 D.note='Spreading: Box-car'; 0340 return 0341 end 0342 0343 error('TYPE of spreading function unknown') 0344 0345 return 0346 0347 function x=fcof2par(r1,type) 0348 %FCOF2PAR Fourier coefficients to distribution parameter 0349 0350 mvrs = version;ix=find(mvrs=='.'); 0351 mvnr = str2num(mvrs(1:ix(2)-1)); 0352 0353 x = zeros(size(r1)); 0354 0355 switch type(1), 0356 case 's' % sech2 0357 x0 = [linspace(eps,5,513) linspace(5.0001,100)]'; 0358 x0 = interp1(0.5*pi./(x0.*sinh(.5*pi./x0)),x0,r1(:)); 0359 k1 = find(isnan(x0)); 0360 if any(k1), 0361 x0(k1)= 0; 0362 x0(k1)=max(x0); 0363 end 0364 0365 ix0 = find(r1~=0); 0366 if mvnr>5.2, 0367 for ix=ix0, 0368 x(ix) = abs(fzero('(0.5*pi./(sinh(.5*pi./x)))-x.*(P1)',x0(ix),optimset,r1(ix))); 0369 end 0370 else 0371 for ix=ix0, 0372 x(ix) = abs(fzero('(0.5*pi./(sinh(.5*pi./x)))-x.*(P1)',x0(ix),[],[],r1(ix))); 0373 end 0374 end 0375 case 'm', % mises 0376 x0 = [linspace(0,10,513) linspace(10.00001,100)]'; 0377 x0 = interp1(besseli(1,x0,1)./besseli(0,x0,1),x0,r1(:)); 0378 k1 = find(isnan(x0)); 0379 if any(k1), 0380 x0(k1)= 0; 0381 x0(k1)=max(x0); 0382 end 0383 0384 if mvnr>5.2, 0385 for ix=1:length(r1), 0386 x(ix) = fzero('besseli(1,x,1)./besseli(0,x,1)-P1',x0(ix),optimset,r1(ix)); 0387 end 0388 else 0389 for ix=1:length(r1), 0390 x(ix) = fzero('besseli(1,x,1)./besseli(0,x,1)-P1',x0(ix),[],[],r1(ix)); 0391 end 0392 end 0393 case 'b', % Box-char 0394 % Initial guess 0395 x0 = linspace(0,pi+0.1,1025).'; 0396 x0 = interp1(sinc(x0/pi),x0,r1(:)); 0397 x = x0.'; 0398 0399 % Newton-Raphson 0400 dx = ones(size(r1)); 0401 iy=0; 0402 max_count=100; 0403 ix =find(x); 0404 while (any(ix) & iy<max_count) 0405 xi=x(ix); 0406 dx(ix) = (sin(xi) -xi.*r1(ix)) ./(cos(xi)-r1(ix)); 0407 xi = xi - dx(ix); 0408 % Make sure that the current guess is larger than zero and less than pi 0409 %x(ix) = xi + 0.1*(dx(ix) - 9*xi).*(xi<=0) +... 0410 % 0.38*(dx(ix)-6.2*xi +6.2).*(xi>pi) ; 0411 % Make sure that the current guess is larger than zero. 0412 x(ix) = xi + 0.5*(dx(ix) - xi).* (xi<=0); 0413 0414 iy=iy+1; 0415 ix=find((abs(dx) > sqrt(eps)*abs(x)) & abs(dx) > sqrt(eps)); 0416 %disp(['Iteration ',num2str(iy),... 0417 %' Number of points left: ' num2str(length(ix)) ]), 0418 end 0419 if iy == max_count, 0420 disp('Warning: fcof2par did not converge.'); 0421 str = 'The last step was: '; 0422 outstr = sprintf([str,'%13.8f'],dx(ix)); 0423 fprintf(outstr); 0424 end 0425 x(x==0)=eps; % Avoid division by zero 0426 end 0427 %semilogy(x,abs(x(:)-x0)) 0428 %plot(x,x(:)-x0) 0429 0430 0431 0432 return 0433 0434 0435 0436 if 0, %old call kept just in case 0437 switch lower(type(1)) 0438 case {'s'}, % sech-2 0439 if def==0, 0440 data2=2; 0441 else 0442 def=1;data2=[2.61 2.28 0.52 1.3 -1.3 0.56 0.95 1.6]; 0443 end 0444 case {'b'}, % box-car 0445 def=0;data2=pi/3; 0446 case {'p'}, % poisson 0447 def=0;data2=0.7; 0448 case {'m'}, % von mises 0449 def=0;data2=2; 0450 otherwise, % cos-2s 0451 type='cos2s'; 0452 switch def 0453 case 0, data2=15; 0454 case 2, data2=[6.97 9.77 0.52 4.06 -2.52 0 1 inf]; 0455 otherwise, %def =1 0456 def=1; data2=[15 15 0.52 5 -2.5 0 1 inf]; 0457 end 0458 end 0459 end 0460 0461 % box-char 0462 ix0 = find(r1<1); 0463 if mvnr>5.2, 0464 for ix=ix0, 0465 x(ix) = abs(fzero('sin(x)-x.*P1',x0(ix) ,optimset,r1(ix))); 0466 end 0467 else 0468 for ix=ix0, 0469 x(ix) = abs(fzero('sin(x)-x.*P1',x0(ix) ,[],[],r1(ix))); 0470 end 0471 end 0472 0473 0474 0475
Comments or corrections to the WAFO group