% function findzeros3
% for use on output from flexural models of volcano loading (in cylindrical
% axisymmetric coordinates r,z,theta)
% find all zero intersections
%
% warning: as currently constitued, this function may have problems with
% consecutive exact zeroes.
function [rinter] = findzeros3(nzeros)
% input arguments:
% vectors to be integrated are passed in global vectors below
% nzeros: number of zeros to be found, starting from first
% matrix entry and working in direction of increasing radius
% output quantities
% rinter: a vector containing the zero crossings (intercepts)
% parameters passed in global arrays
% radinp: radius vector (should be from zero to max radius)
% zinp: vector of quantity being integrated
%
global radinp zinp zoffset;
rinter = zeros(1,nzeros);
lenz = length(zinp); % should also be length of radinp
zsign = sign(zinp);
warnfind = find(zsign == 0);
if (~(isempty(warnfind)))
disp(['Warning: ' num2str(length(warnfind)) ' exact zeroes found by findzeros3'])
end
% sign function = 1 for > 0, 0 for = 0, and -1 for < 0
zsndiff = diff(zsign); % should be of length lenz-1
zcfindex = find(zsndiff ~= 0);
zcfnum = length(zcfindex);
% nonzero entries in zsndiff correspond to places where zero has been
% reached or crossed
strzcfnum = num2str(zcfnum);
if nzeros > zcfnum % too many zeros requested, get outa here
error(['user requested ' num2str(nzeros) ' zeros, at most ' strzcfnum 'found'])
else
disp(['findzeros3 found ' strzcfnum ' zeros']);
end
zsncount = 1; % counter for nonzero entries in zsndiff (which we need because
% of possibility of exactly hitting a zero (see below)
nzfound = 0; % count of number of zeros actually found
while (nzfound < nzeros)
% find next occurence of a zero or zero crossing as determined by the
% matricies zsign and zsndiff
if abs(zsndiff(zcfindex(zsncount))) == 2
%difference will be +/- 2 if going from pos. to neg. or vice versa
% use "fzero" to find the zero in this area
% interz is a interpolating function that has access to the
% global variables listed above
rinter(nzfound+1) = fzero('interz',radinp(zcfindex(zsncount)));
zsncount = zsncount + 1;
elseif abs(zsndiff(zcfindex(zsncount))) == 1
% difference will be +/- if going from a zero to pos. to neg.:
% i.e. we have already found the zero
rinter(nzfound+1) = radinp(zcfindex(zsncount+1));
% skip over next entry in zsndiff (which should be another +/- 1) by
% incrementing zsncount by 2
% exception: 1) if zero is found at first r value, increment by 1
if zcfindex(zsncount) == 1
zsncount = zsncount + 1;
else
zsncount = zsncount + 2;
end
if zsncount > zcfnum
error('not enough potential zero crossings in findzeros3')
end
end % if loop
nzfound = nzfound + 1;
end% while loop