Home > wafo > spec > spreading.m

spreading

PURPOSE ^

Directional spreading functions

SYNOPSIS ^

D = spreading(th,type,th0,data,w,def)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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