Home > wafo > spec > spreading.m

## DESCRIPTION

 SPREADING Directional spreading functions

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:
data =  [spa spb wc ma mb wlim],
wc = cut over frequency (usually the peak frequency, wp or fp)
m  = shape parameter defining the frequency dependent
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
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];
% Frequency dependent direction
th0 = linspace(0,pi/2,257)';

See also  mkdspec, wspecplot, spec2spec

## CROSS-REFERENCE INFORMATION

This function calls:
 sinc Sin(pi*x)/(pi*x) function. besseli Modified Bessel function of the first kind. error Display message and abort function. fft Discrete Fourier transform. fzero Scalar nonlinear zero finding. gamma Gamma function. interp1 1-D interpolation (table lookup) linspace Linearly spaced vector. lower Convert string to lowercase. optimset Create/alter OPTIM OPTIONS structure. sprintf Write formatted data to string. str2num Convert string matrix to numeric array. strncmp Compare first N characters of strings. struct Create or convert to structure array. version MATLAB version number. warning Display warning message; disable or enable warning messages.
This function is called by:
 Chapter1 % CHAPTER1 demonstrates some applications of WAFO Chapter2 % CHAPTER2 Modelling random loads and stochastic waves Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics demospec Loads a precreated spectrum of chosen type mkdspec Make a directional spectrum wafofig4 Directional spectra using Thorsethaugen and cos2s spreading

## SOURCE CODE

0001 function D = spreading(th,type,th0,data,w,def)
0003 %
0004 % CALL:  D = spreading(th,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:
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
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];
0092 %        % Frequency dependent direction
0093 %   th0 = linspace(0,pi/2,257)';
0095 %
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
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
0140 % revised pab 6 April 2001
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
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
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);
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);
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;
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));
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);
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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003