Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
144 lines (131 sloc) 4.81 KB
function [newval,req,dispos,avma]=discinter(pos,val,req)
% [newval,req,dispos,avma]=DISCINTER(pos,val,req)
%
% For non-monotonous vectors (such as with discontinuities, specified
% by two or three times the same value), interpolates, i.e. finds the
% interpolated value of input values.
% Returns different values than requested!
%
% Returns NaN outside the range, as in regular interpolation.
% Also for regular interpolation, i.e. without discontinuities.
% 'pos' must not start with a discontinuity
%
% The 'newval' belong to an upgoing sorted 'req', possibly with
% fewer entries than the input 'req', if the latter was discontinuous
% either as 2p or 3p.
%
% Discontinuities at zero are treated as normal ones, so they might
% end up having negative values. So when doing inner products
% and stuff, make sure to specify values a little bit out of the range.
%
% Everything discontinuity related will get differences (intervals)
% of 0.5 and 1.5 times smint10, all the rest are real intervals.
% Changed smint/10 to smallest of 1 m or smint/10.
%
% Also returns the vector 'dispos' from DISCIDENT(POS)
% Also returns the averaging matrix 'avma' from DISCIDENT(POS)
%
% See DISCINTEREX, and NODINTER
%
% For example, an request vector [15 30 80 140 200 300 400 400] has
% a discontinuity at 400. But the position/value pair has discons at
% 15 80 and 400, specified as a repeated value (2p).
% The one at 400 is common, but 15 and 80 aren't
% The result is: you get values at
% [14 14.5 16 30 79 80.5 81 140 200 300 399 401]
% If you ask what the value at 15 is, you can't get a straight answer, but
% get three values instead.
% If you ask what the value is at 400 and 400, you get the values at 399 and 401
%
% Last updated by fjsimons-at-alum.mit.edu, October 18th, 2002
[pos,val,req]=deal(pos(:),val(:),req(:));
if pos(1)>pos(end)
pos=flipud(pos);
val=flipud(val);
fli=1;
else
fli=0;
end
% If requested values are 2p-3p-discontinuities in the request array, then
% destroy, not even leaving the middle value intact for 3p.
difreq=diff(req);
smint=min(difreq(~~difreq));
smint10=min(1,smint/10);
if sum(~difreq)
[a,b]=discident(req);
% Leave quasi-3p for quadrature; work with NaN's
% Smallest interval of req
req(~isnan(a))=req(~isnan(a))-(smint10)*a(~isnan(a));
% For 3p, adjust the middle value, for 2p, leave intact
req(~b)=req(~b)-(smint10)/2;
end
% This makes [400 400] into [399 401]
req=sort(req);
outsid=req<min(pos) | req>max(pos);
reqout=req(outsid);
req=req(~outsid);
outvol=repmat(NaN,sum(outsid),1);
% What are the discontinuities of 'pos'?
[c,d,posdisc,avma]=discident(pos);
dispos=c;
% If it was flipped, unflip for output
if fli==1
dispos=flipud(dispos);
avma=flipud(avma);
end
% So now the requested vector 'req' is free of discontinuities.
% Now need to make sure that the values of 'req'
% that would be equal to an element of 'pos' are treated
% properly. But 'req' can still have values equal to a discontinuity in pos
% this was the original point of this program.
flag1=0;
if ~isempty(intersect(pos,req))
[pdi,inp,inr]=intersect(posdisc,req);
if ~isempty(pdi)
% Remediate in the request vector, make it want a "3p-discontinuity"
% 3p for integration purposes
req(inr)=pdi(:)-repmat(smint10,length(pdi),1);
req=[req ; pdi(:)+repmat(smint10,length(pdi),1)];
req=[req ; pdi(:)+repmat(smint10/2,length(pdi),1)];
req=sort(req);
end
flag1=1;
% Else you treat those values
% I can take them out, knowing they are not discs in either
% vector, lookup their unambiguous value and put them in at
% the end, knowing the output was sorted anyway.
[reqvol,posinpos,posinreq]=intersect(pos,req);
newvol=val(posinpos);
req=req(~ismember(1:length(req),posinreq));
end
% If there are 2p or 3p discontinuities in the positions
if sum(~diff(pos))
posorg=pos;
% This removes the discontinuities from 'pos';
% it is assumed that 'pos' and 'req' have no common members.
% This is where 'req' gets sorted!
pos=unique([pos ; req]);
ins=[1:length(pos)]';
% reqpos is the index of req in [pos-disc+req]
reqpos=ins(ismember(pos,req));
% reqpos is index of one above req in [pos-disc]
reqpos=reqpos-[0:length(reqpos)-1]';
% This need to be the position shift needed to translate
% the original pos index into pos-disc
[a,ind,c]=degamini(posorg);
ind=cumsum([0 ind-1]); ind=ind(:);
% What is the index, in the original value matrix, of the
% parameter just down the one we are requesting?
% And then ups-1 is the one right below
ups=reqpos+ind(reqpos);
newval=val(ups-1)+(req-posorg(ups-1))./...
(posorg(ups)-posorg(ups-1)).*(val(ups)-val(ups-1));
else
newval=interp1(pos,val,req,'linear');
end
if flag1==1
outrows=sortrows([req newval ; reqvol newvol ; reqout outvol]);
else
outrows=sortrows([req newval ; reqout outvol]);
end
[req,newval]=deal(outrows(:,1),outrows(:,2));