--- /dev/null
+## Copyright (c) 1998-2011 Andrew V. Knyazev <andrew.knyazev@ucdenver.edu>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without",
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+## 3 Neither the name of the author nor the names of its contributors may be
+## used to endorse or promote products derived from this software without
+## specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+% function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
+%
+% Best polynomial approximation in discrete uniform norm
+%
+% INPUT VARIABLES:
+%
+% M : degree of the fitting polynomial
+% N : number of data points
+% X(N) : x-coordinates of data points
+% Y(N) : y-coordinates of data points
+% K : character of the polynomial:
+% K = 0 : mixed parity polynomial
+% K = 1 : odd polynomial ( X(1) must be > 0 )
+% K = 2 : even polynomial ( X(1) must be >= 0 )
+% EPSH : tolerance for leveling. A useful value for 24-bit
+% mantissa is EPSH = 2.0E-7
+% MAXIT : upper limit for number of exchange steps
+% REF0(M2): initial alternating set ( N-vector ). This is an
+% OPTIONAL argument. The length M2 is given by:
+% M2 = M + 2 , if K = 0
+% M2 = integer part of (M+3)/2 , if K = 1
+% M2 = 2 + M/2 (M must be even) , if K = 2
+%
+% OUTPUT VARIABLES:
+%
+% A : polynomial coefficients of the best approximation
+% in order of increasing powers:
+% p*(x) = A(1) + A(2)*x + A(3)*x^2 + ...
+% REF : selected alternating set of points
+% HMAX : maximum deviation ( uniform norm of p* - f )
+% H : pointwise approximation errors
+% R : total number of iterations
+% EQUAL : success of failure of algorithm
+% EQUAL=1 : succesful
+% EQUAL=0 : convergence not acheived
+% EQUAL=-1: input error
+% EQUAL=-2: algorithm failure
+%
+% Relies on function EXCH, provided below.
+%
+% Example:
+% M = 5; N = 10000; K = 0; EPSH = 10^-12; MAXIT = 10;
+% X = linspace(-1,1,N); % uniformly spaced nodes on [-1,1]
+% k=1; Y = abs(X).^k; % the function Y to approximate
+% [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
+% p = polyval(A,X); plot(X,Y,X,p) % p is the best approximation
+%
+% Note: using an even value of M, e.g., M=2, in the example above, makes
+% the algorithm to fail with EQUAL=-2, because of collocation, which
+% appears because both the appriximating function and the polynomial are
+% even functions. The way aroung it is to approximate only the right half
+% of the function, setting K = 2 : even polynomial. For example:
+%
+% N = 10000; K = 2; EPSH = 10^-12; MAXIT = 10; X = linspace(0,1,N);
+% for i = 1:2
+% k = 2*i-1; Y = abs(X).^k;
+% for j = 1:4
+% M = 2^j;
+% [~,~,HMAX] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
+% approxerror(i,j) = HMAX;
+% end
+% end
+% disp('Table 3.1 from Approximation theory and methods, M.J.D.POWELL, p. 27');
+% disp(' ');
+% disp(' n K=1 K=3');
+% disp(' '); format short g;
+% disp([(2.^(1:4))' approxerror']);
+%
+% ALGORITHM:
+%
+% Computation of the polynomial that best approximates the data (X,Y)
+% in the discrete uniform norm, i.e. the polynomial with the minimum
+% value of max{ | p(x_i) - y_i | , x_i in X } . That polynomial, also
+% known as minimax polynomial, is obtained by the exchange algorithm,
+% a finite iterative process requiring, at most,
+% n
+% ( ) iterations ( usually p = M + 2. See also function EXCH ).
+% p
+% since this number can be very large , the routine may not converge
+% within MAXIT iterations . The other possibility of failure occurs
+% when there is insufficient floating point precision for the input
+% data chosen.
+%
+% CREDITS: This routine was developed and modified as
+% computer assignments in Approximation Theory courses by
+% Prof. Andrew Knyazev, University of Colorado Denver, USA.
+%
+% Team Fall 98 (Revision 1.0):
+% Chanchai Aniwathananon
+% Crhistopher Mehl
+% David A. Duran
+% Saulo P. Oliveira
+%
+% Team Spring 11 (Revision 1.1): Manuchehr Aminian
+%
+% The algorithm and the comments are based on a FORTRAN code written
+% by Joseph C. Simpson. The code is available on Netlib repository:
+% http://www.netlib.org/toms/501
+% See also: Communications of the ACM, V14, pp.355-356(1971)
+%
+% NOTES:
+%
+% 1) A may contain the collocation polynomial
+% 2) If MAXIT is exceeded, REF contains a new reference set
+% 3) M, EPSH and REF can be altered during the execution
+% 4) To keep consistency to the original code , EPSH can be
+% negative. However, the use of REF0 is *NOT* determined by
+% EPSH< 0, but only by its inclusion as an input parameter.
+%
+% Some parts of the code can still take advantage of vectorization.
+%
+% Revision 1.0 from 1998 is a direct human translation of
+% the FORTRAN code http://www.netlib.org/toms/501
+% Revision 1.1 is a clean-up and technical update.
+% Tested on MATLAB Version 7.11.0.584 (R2010b) and
+% GNU Octave Version 3.2.4
+
+% $Revision: 1.1 $ $Date: 2011/08/3 $
+
+% ************************************ beginning of POLYFITINF
+function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
+
+ % Preassign output variables A,REF,HMAX,H,R,EQUAL in case of error return
+ A = []; REF = []; HMAX = []; H = []; R = 0; EQUAL = -2;
+ %%%% end preassignment
+
+ % Setting M with respect to K
+ MOLD = M;
+
+ switch K
+ case 1
+ K0 = 0;
+ K1 = 1;
+ Q1 = 1;
+ Q2 = 2;
+ M = (M-Q1)/2;
+ case 2
+ K0 = 0;
+ K1 = 0;
+ Q1 = 0;
+ Q2 = 2;
+
+ % If the user has input odd M, but wants an even polynomial,
+ % subtract 1 from M to prevent errors later. The outputs should be
+ % mathematically equivalent.
+ if mod(M,2) == 1
+ M = M-1;
+ end
+
+ M = (M-Q1)/2;
+ otherwise
+ if (K ~= 0)
+ warning('polyfitinf:MixedParity','Using mixed parity polynomial...');
+ end
+ K0 = 1;
+ K1 = 0;
+ Q1 = 0;
+ Q2 = 1;
+ end
+
+ P = M + 2;
+
+ % Check input data consistency
+
+ if ( (length(X) ~= N) || (length(Y) ~= N) )
+ error('Input Error: check data lengths');
+ end
+
+ if (P > N)
+ error('Input Error: insufficient data points');
+ end
+
+ if (M < 0)
+ error('Input Error: insufficient degree');
+ end
+
+ if ( (K == 2) && (X(1) < 0) ) || ( (K == 1) && (X(1) <= 0) )
+ error('Input Error: X(1) inconsistent with parity');
+ end
+
+ if any(diff(X)<0)
+ error('Input Error: Abscissae out of order');
+ end
+
+ ITEMP = MOLD + 1;
+
+ A = zeros(1,ITEMP);
+
+
+ ITEMP = P + 2;
+ Z = zeros(1,ITEMP);
+ Z(1) = 0;
+ Z(ITEMP) = N + 1;
+
+ EPSH = abs(EPSH);
+
+ % Read initial reference set into Z, if available.
+
+ if (nargin == 8)
+ J = 0;
+
+ Z(2:(P+1))= REF0(1:P);
+
+ % Check if REF is monotonically increasing
+ if ( any(diff(REF0) < 0) || any(REF0 > J) )
+ error('Input Error : Bad initial reference set');
+ end
+
+ else
+
+ % Loads Z with the points closest to the Chebychev abscissas
+
+ X1 = X(1);
+ XE = X(N);
+
+ % Setting parity-dependent parameters
+
+ if (K0 == 1)
+ XA = XE + X1;
+ XE = XE - X1;
+ Q = pi/(M + 1.0);
+ else
+ XA = 0.;
+ XE = XE + XE;
+ ITEMP = 2*(M+1) + Q1;
+ Q = pi/(ITEMP);
+ end
+
+ % Calculate the J-th Chebyshev abcissa and load Z(J+1)
+ % with the appropriate index from the data abcissas
+
+ for JJ = 1:P
+ J = P + 1 - JJ;
+ X1 = XA + XE*( cos(Q*(P-J)) );
+ ITEMP = J + 2;
+ R = Z(ITEMP);
+ HIGH = R - 1;
+ FLAG = 1;
+ if (HIGH >= 2)
+ II = 2;
+ while ( (II <= HIGH) && (FLAG == 1) )
+ I = HIGH + 2 - II;
+ ITEMP = I - 1;
+
+ % If the Chebyschev abscissa is bracketed by
+ % two input abcissas, get out of the while loop
+
+ if (X(I)+X(ITEMP) <= X1)
+ FLAG = 0;
+ end
+ II = II + 1;
+ end
+ end
+
+ if (FLAG == 1)
+ I = 1;
+ end
+ ITEMP = J + 1;
+
+ if (I < R)
+ Z(ITEMP) = I;
+ else
+ Z(ITEMP) = R - 1;
+ end
+ end
+
+ % If the lower Chebyshev abcissas are less than X(1),
+ % load the lower elements of Z with the lowest points
+
+ IND = find(Z(2:end) >= (1:(length(Z)-1)));
+
+ try TEMP = IND(1); % If IND is empty, do nothing.
+ catch exception % The catch will be that IND is an empty array.
+
+ if strcmpi(exception.identifier,'MATLAB:badsubscript')
+ % This will be the exception. Do nothing.
+ end
+ end
+
+ if TEMP~=1
+ Z(2:TEMP) = (1:(TEMP-1))';
+ end
+
+
+ end
+
+ % M1 entry. Initialize variables to prepare for exchange iteration
+
+ ITEMP = M + 1;
+
+
+ % Zero the AA array
+
+
+ AA = zeros(1,ITEMP);
+
+
+ % Load H with the ordinates and XX(I) with the abscissas if the
+ % polynomial is mixed . If it is even or odd , load XX with the
+ % squares of the abscissas.
+
+ H(1:N) = Y(1:N);
+ if (K0 <=0)
+ XX(1:N) = X(1:N).^2;
+ else
+ XX(1:N) = X(1:N);
+ end
+
+ B1 = 0;
+ B2 = 0;
+ B3 = 0;
+ R = -1;
+ T = 0.;
+
+ % Iteration entry. R is the iteration index
+
+ C = zeros(1,P);
+ D = zeros(1,P);
+ DAA = zeros(1,M+1);
+
+ FLAG = 1;
+ while ( (R < MAXIT) && (FLAG == 1) )
+
+ R = R + 1; % LABEL 350
+ %S = 1.;
+
+ % Computation of div. differences schemes
+
+ if (K1 > 0)
+
+ % If the polynomial is mixed or even:
+ %for I = 1:P
+ % S = -S;
+ % ITEMP = I + 1;
+ % J = Z(ITEMP);
+ % Q = X(J);
+ % C(I) = (H(J) + S*T)/Q;
+ % D(I) = S/Q;
+ %end
+
+ I = (1:P);
+ S = (-1).^I;
+ ITEMP = I+1;
+ J = Z(ITEMP);
+ C(I) = (H(J) + S*T)./X(J);
+ D(I) = S./Q;
+ clear I ITEMP S J
+
+ else
+
+ % If the polynomial is odd:
+ %for I = 1:P
+ % S = -S;
+ % ITEMP = I + 1;
+ % ITEMP = Z(ITEMP);
+ % C(I) = H(ITEMP) + S*T;
+ % D(I) = S;
+ %end
+
+ I = (1:P);
+ S = (-1).^I;
+ ITEMP = I+1;
+ C(I) = H( Z(ITEMP) ) + S.*T;
+ D(I) = S;
+ clear I ITEMP S
+
+ end
+
+ for I = 2:P
+ for JJ = I:P
+ J = P + I - JJ;
+ ITEMP = J + 1;
+ ITEMP = Z(ITEMP);
+ QD = XX(ITEMP);
+ ITEMP = 2 + J - I;
+ ITEMP = Z(ITEMP);
+ QD = QD - XX(ITEMP);
+ ITEMP = J - 1;
+ C(J) = (C(J)-C(ITEMP))/QD;
+ D(J) = (D(J)-D(ITEMP))/QD;
+ end
+ end
+
+ DT = -C(P)/D(P);
+ T = T + DT;
+
+ % Computation of polynomial coefficients
+
+ HIGH = M + 1;
+ for II = 1:HIGH
+ I = HIGH - II;
+ ITEMP = I + 1;
+ DAA(ITEMP) = C(ITEMP) + DT*D(ITEMP);
+ ITEMP = I + 2;
+ ITEMP = Z(ITEMP);
+ QD = XX(ITEMP);
+ LOW = I + 1;
+ if (M >= LOW)
+ DAA(LOW:M) = DAA(LOW:M) - QD*DAA(((LOW:M)+1));
+ end
+ end
+
+ AA(1:HIGH) = AA(1:HIGH) + DAA(1:HIGH);
+
+ % Evaluation of the polynomial to get the approximation errors
+
+ MAXX = 0.;
+ H = zeros(1,N);
+ for I = 1:N
+ SD = AA(HIGH);
+ QD = XX(I);
+ if (M > 0)
+ for J = M:-1:1
+ SD = SD*QD + AA(J);
+ end
+ end
+ if (K1 > 0)
+ % If the polynomial is odd, multiply SD by X(I)
+ SD = SD*X(I);
+ end
+
+ QD = Y(I) - SD;
+ H(I) = Y(I) - SD;
+
+ if (abs(QD) > MAXX)
+ % Load MAXX with the largest magnitude
+ % of the approximation array
+ MAXX = abs(QD);
+ end
+ end
+
+ % Test for alternating signs
+
+ ITEMP = Z(2);
+
+ if (H(ITEMP) == 0.)
+
+ % This represents a case where the polynomial
+ % exactly predicts a data point
+
+ warning('polyfitinf:Collocation','Collocation has occured.');
+ if (B3 > 0)
+ B3 = -1;
+ FLAG = 0;
+ else
+ B3 = 1;
+ if (EPSH < MAXX)
+ warning('polyfitinf:AnotherTry','1 more attempt with middle points');
+ LOW = (N+1)/2 - (P+1)/2 + 1;
+ HIGH = LOW + P;
+ Z(LOW:HIGH) = ( (LOW:HIGH) -1);
+
+ else
+ disp('Normal Exit.');
+ FLAG = 0;
+ end
+ end
+ else
+
+ if (H(ITEMP) > 0.)
+ J = -1;
+ else
+ J = 1;
+ end
+
+ I = 2;
+ FLAG2 = 1;
+ while ( (I <= P) && (FLAG2 == 1) )
+ ITEMP = I + 1;
+ ITEMP = Z(ITEMP);
+ if (H(ITEMP) == 0.)
+ J = 0;
+ warning('polyfitinf:Collocation','Collocation has occured.');
+ if (B3 > 0)
+ B3 = -1;
+ FLAG = 0;
+ else
+ B3 = 1;
+ if (EPSH < MAXX)
+ warning('polyfitinf:AnotherTry','1 more attempt with middle points');
+ LOW = (N+1)/2 - (P+1)/2 + 1;
+ HIGH = LOW + P;
+ Z(LOW:HIGH) = ( (LOW:HIGH) -1);
+ else
+ disp('Normal Exit.');
+ FLAG = 0;
+ end
+ end
+ FLAG2 = 0;
+ else
+ if (H(ITEMP) < 0)
+ JJ = -1;
+ else
+ JJ = 1;
+ end
+ if (J~=JJ)
+
+ % Error entry: bad accuracy for calculation
+
+ B1 = 1;
+ FLAG2 = 0;
+ FLAG = 0;
+ else
+ J = -J;
+ end
+ end
+
+ I = I + 1;
+
+ end % end of while
+
+ % Search for another reference
+
+ if (FLAG2*FLAG == 1)
+
+ [H,Z,EQUAL] = exch(N, P, EPSH, H, Z);
+ if (EQUAL > 0)
+ FLAG = 0;
+ else
+ if (R >= MAXIT)
+ B2 = 1;
+ FLAG = 0;
+ end
+ end
+
+ end
+
+ end % end of if over H(ITEMP)
+
+ end; % end of iteration loop
+
+ % M2 entry; load output variables and return
+
+ HIGH = M + 1;
+
+ % Load the coefficients into A array
+
+ A(Q1 + Q2*(((1:HIGH)-1)) + 1) = AA(1:HIGH);
+
+ % Load REF with the final reference points
+
+ REF(1:P) = Z((1:P) + 1);
+
+ HMAX = MAXX;
+
+ if (B3 < 0)
+ EQUAL = -2;
+ warning('polyfitinf:Collocation','polyfitinf terminates');
+ end
+ if (B1 > 0)
+ EQUAL = -2;
+ warning('polyfitinf:NoAlternatingSigns','Alternating signs not observed');
+ end
+ if (B2 > 0)
+ EQUAL = 0;
+ warning('polyfitinf:MaxIterationsReached','MAXIT was reached, current ref. set saved in REF.');
+ end
+
+ % Reverse the order of A to make it compatible with MATLAB'S polyval() function.
+ A = A(end:-1:1);
+
+
+endfunction
+% ****************************************** end of POLYFITINF
+
+function [H,Z,EQUAL] = exch(N, P, EPSH, H, Z)
+% function [H,Z,EQUAL] = exch(N, P, EPSH, H, Z)
+%
+% EXCH: exchange algorithm
+%
+% INPUT VARIABLES:
+% N : number of data points
+% P : number of reference points
+% EPSH : tolerance for leveling.
+% Z : old reference indices
+%
+% OUTPUT VARIABLES:
+% H : pointwise approximation errors
+% Z : new reference indices
+% EQUAL : EQUAL=1 : normal exchange
+% EQUAL=0 : old and new references are equal
+%
+% CREDITS: This routine was developed and modified as
+% computer assignments in Approximation Theory courses by
+% Prof. Andrew Knyazev, University of Colorado Denver, USA.
+%
+% Team Fall 98 (Revision 1.0):
+% Chanchai Aniwathananon
+% Crhistopher Mehl
+% David A. Duran
+% Saulo P. Oliveira
+%
+% Team Spring 11 (Revision 1.1): Manuchehr Aminian
+%
+% The algorithm and the comments are based on a FORTRAN code written
+% by Joseph C. Simpson. The code is available on Netlib repository:
+% http://www.netlib.org/toms/501
+% See also: Communications of the ACM, V14, pp.355-356(1971)
+%
+% Revision 1.0 from 1998 is a direct human translation of
+% the FORTRAN code http://www.netlib.org/toms/501
+% Revision 1.1 is a clean-up and technical update.
+% Tested on MATLAB Version 7.11.0.584 (R2010b) and
+% GNU Octave Version 3.2.4
+
+% License: BSD
+% Copyright 1998-2011 Andrew V. Knyazev
+% $Revision: 1.1 $ $Date: 2011/05/17 $
+
+% ************************************ beginning of exch
+
+ EQUAL = 0;
+ L = 0;
+ ITEMP = Z(2);
+
+ % SIG is arbitrarily chosen equal to the sign of the input
+ % point. This will be adjusted later if necessary.
+
+ if (H(ITEMP) <= 0)
+ SIG = 1.;
+ else
+ SIG = -1.;
+ end
+
+ % The next loop prescans Z to insure it is a proper choice, i.e
+ % resets Z if necessary so that maximum error points are chosen,
+ % given the sign convention mentioned above. In order to work
+ % properly, this section requires Z(1) = 0 and Z(P+2) = N + 1 .
+
+ for I = 1:P
+ MAXX = 0.;
+ SIG = -SIG;
+ ITEMP = I + 2;
+ ZE = Z(ITEMP) - 1;
+ LOW = Z(I) + 1;
+
+ % Scan the open point interval containing only the 1th initial
+ % reference point. In the interval pick the point with largest
+ % magnitude and correct sign. Most of the sorting occurs in
+ % this section. SIG contains the sign assumed for H(I).
+
+ for J = LOW:ZE
+ if (SIG*(H(J)-MAXX) > 0)
+ MAXX = H(J);
+ INDEX = J;
+ end
+ end
+ ITEMP = I + 1;
+ ITEMP = Z(ITEMP);
+ MAXL = abs(MAXX);
+
+ % If the MAX error is significantly greater than the
+ % input point, switch to this point.
+
+ if (abs( MAXX - H(ITEMP) )/MAXL > EPSH)
+ ITEMP = I + 1;
+ Z(ITEMP) = INDEX;
+ L = 1;
+ end
+ end
+ %
+ MAXL = 0.;
+ MAXR = 0.;
+ ITEMP = P + 1;
+ LOW = Z(ITEMP) + 1;
+ %
+ if (LOW <= N)
+
+ % Find the error with largest abs value and proper sign
+ % from among the points above the last reference point.
+ % This section is necessary because the set of points
+ % chosen may begin with the wrong sign alternation.
+
+ for J = LOW:N
+ if (SIG*(MAXR-H(J)) > 0)
+ MAXR = H(J);
+ INDR = J;
+ end
+ end
+ end
+
+ % Find the error with largest abs value and proper sign
+ % from among the points below the 1st reference point.
+ % This section is necessary by the same reason as above.
+
+ ITEMP = Z(2);
+ HZ1 = H(ITEMP);
+ HIGH = ITEMP -1;
+ if (HIGH > 0)
+ if (HZ1 < 0)
+ SIG = -1.;
+ elseif (HZ1 == 0)
+ SIG = 0.;
+ else
+ SIG = 1.;
+ end
+
+ for J = 1:HIGH
+ if (SIG*(MAXL-H(J)) > 0)
+ MAXL = H(J);
+ INDL = J;
+ end
+ end
+
+ end
+
+ % MAXL and MAXR contain the magnitude of the significant
+ % errors outside the reference point set. If either is
+ % zero, the reference point set extends to the end point
+ % on that side of the interval.
+
+ MAXL = abs(MAXL);
+ MAXR = abs(MAXR);
+ HZ1 = abs(HZ1);
+ ITEMP = P + 1;
+ ITEMP = Z(ITEMP);
+ HZP = abs(H(ITEMP));
+
+ % L = 0 implies that the previous prescan did not change
+ % any points. If L = 0 and MAXL, MAXR are not significant
+ % if compared with upper and lower reference point errors,
+ % respectively, use the EQUAL exit.
+
+ FLAG1 = 1;
+ if (L == 0)
+ if ( (MAXL == 0) || (EPSH >= (MAXL-HZP)/MAXL) )
+ if ( (MAXR == 0) || (EPSH >= (MAXR-HZ1)/MAXR) )
+ FLAG1 = 0;
+ EQUAL = 1;
+ end
+ end
+ end
+
+ if ( (MAXL == 0) && (MAXR == 0) )
+ FLAG1 = 0;
+ end
+
+ if ( (MAXL > MAXR) && (MAXL <= HZP) && (MAXR < HZ1) )
+ FLAG1 = 0;
+ end
+
+ if ( (MAXL <= MAXR) && (MAXR <= HZ1) && (MAXL < HZP) )
+ FLAG1 = 0;
+ end
+
+ % If a point outside the present reference set must be
+ % included, (i.e. the sign of the 1st point previously
+ % assumed is wrong) shift to the side of largest
+ % relative error first.
+
+ if (FLAG1 == 1)
+
+ FLAG2 = 1;
+
+ if ( (MAXL > MAXR) && (MAXL > HZP) )
+ FLAG2 = 0;
+ end
+
+ if ( (MAXL <= MAXR) && (MAXR <= HZ1) )
+ FLAG2 = 0;
+ end
+
+ if (FLAG2 == 1)
+
+ % SHR entry. This section inserts a point from
+ % above the prescan point set
+
+ INDEX = Z(2);
+
+ % shift point set down, dropping the lowest point
+
+ Z(2:P) = Z((2:P)+1);
+
+ ITEMP = P + 1;
+
+ % add the next high point
+ Z(ITEMP)=INDR;
+
+ % if MAXL > 0 replace reference points from the left,
+ % stopping the 1st time the candidate for replacement
+ % is greater than in magnitude than the prospective
+ % replacee. Alternation of signs is preserved because
+ % non-replacement immediately terminates the process.
+
+ if (MAXL > 0)
+ I = 2;
+ FLAG3 = 1;
+ while ( (I <= P) && (FLAG3 == 1) )
+ ITEMP = Z(I);
+ if ( abs(H(INDL)) >= abs(H(ITEMP)) )
+ J = Z(I);
+ Z(I) = INDL;
+ INDL = INDEX;
+ INDEX = J;
+ else
+ FLAG3 = 0;
+ end
+ I = I + 1;
+ end
+ end
+
+ else
+
+ % SHL entry. This section inserts a point from below the
+ % prescan point set.
+
+ ITEMP = P + 1 ;
+ INDEX = Z(ITEMP);
+
+ Z((2:P)+1) = Z(2:P);
+
+ % store lowest point in Z(2)
+ Z(2) = INDL;
+
+ % if MAXR > 0 replace reference points from the right,
+ % stopping the 1st time the candidate for replacement
+ % is greater than in magnitude than the prospective
+ % replacee.
+
+ if (MAXR > 0)
+ II = 2;
+ FLAG3 = 1;
+ while ( (II <= P) && (FLAG3 == 1) )
+ I = P + 2 - II;
+ ITEMP = I + 1;
+ HIGH = Z(ITEMP);
+ if ( abs(H(INDR)) >= abs(H(HIGH)) )
+ J = Z(ITEMP);
+ Z(ITEMP) = INDR;
+ INDR = INDEX;
+ INDEX = J;
+ else
+ FLAG3 = 0;
+ end
+ II = II + 1;
+ end
+ end
+
+ end
+
+ end
+
+endfunction
+% ****************************************** end of exch
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff