From sources-request at octave dot org Sun Jun  5 23:13:00 2005
Subject: function to find level crossings
From: "Phillip M. Feldman" <pfeldman at verizon dot net>
To: sources at octave dot org
Date: Sun, 5 Jun 2005 15:29:59 -0500

This is a multi-part message in MIME format.
--------------080500020900030006050004
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit



--------------080500020900030006050004
Content-Type: text/plain;
 name="find_cross.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="find_cross.m"

%                             find_cross.m
%
%                       Phillip M. Feldman, NGST
%
%                       Version 1.1  June 2, 2005
%
% Overview
%
% find_cross() uses a combination of search and interpolation to
% determine the point or points at which a given level is crossed.  If
% there are multiple level crossings, find_cross finds all of the level
% crossings.  If the input sequence "kisses" the specified level but
% does not actually cross it, this is not counted as a crossing.
%
% Calling syntax:
%
% X_cross= find_cross(X, Y, Y_target)
%
% where:
%
% X is a column vector containing a sequence of monotone increasing
% values of the independent variable.
%
% Y is a column vector containing values of the dependent variable.
%
% Y_target is the desired level crossing.
%
% X_cross (the returned value) is the estimated value of X at which the
% dependent variable crosses the level Y_target.


% Revision History

% 1.1 Fixed code to handle situation where Y_target exactly matches
% Y(i), with Y(i-1) and Y(i+1) bracketing Y_target.  Previously, this
% was not being recognized as a crossing.

function X_cross= find_cross(X, Y, Y_target)


% Section 1: Check array dimensions and augment Y array with an extra
% element.

[Xrows Xcols]= size(X);
[Yrows Ycols]= size(Y);

if (Xcols ~= 1 | Ycols ~= 1)
   error('X and Y must be column vectors.');
end

if (Xrows ~= Yrows)
   error('Lengths of X and Y vectors must be the same.');
end

m= Xrows;

% To prevent a subscript-out-of-bounds error in Section 2, augment Y by
% duplicating the last element of the array.  (We are operating on a
% copy of the Y array, so the Y array in the calling program is not
% affected).

Y(m+1)= Y(m);


% Section 2: Search for table entries that bracket the target level
% Y_target.

crossings= 0;
X_cross= []; % to guarantee that array is defined

for i= 1 : m-1

   % Check for up-crossing of level:

   if ( (Y(i)<Y_target & Y(i+1)>Y_target) | ...                % bracket
        (Y(i)<Y_target & Y(i+1)==Y_target & Y(i+2)>Y_target) ) % exact

      crossings= crossings + 1;

      % Estimate crossing point via linear interpolation:

      X_cross(crossings)= ...
        X(i) + (X(i+1) - X(i)) * (Y_target - Y(i)) / (Y(i+1) - Y(i));

   end

   % Check for down-crossing of level:

   if ( (Y(i)>Y_target & Y(i+1)<Y_target) | ...                % bracket
        (Y(i)>Y_target & Y(i+1)==Y_target & Y(i+2)<Y_target) ) % exact

      crossings= crossings + 1;

      % Estimate crossing point via linear interpolation:

      X_cross(crossings)= ...
        X(i) + (X(i+1) - X(i)) * (Y_target - Y(i)) / (Y(i+1) - Y(i));

   end

end

--------------080500020900030006050004--