| [x1,x2,f1,f2,d1,d2,ne,xe,fe,next,ierr]=dchfev(x1,x2,f1,f2,d1,d2,ne,xe,fe,next,ierr); |
function [x1,x2,f1,f2,d1,d2,ne,xe,fe,next,ierr]=dchfev(x1,x2,f1,f2,d1,d2,ne,xe,fe,next,ierr);
%***BEGIN PROLOGUE DCHFEV
%***PURPOSE Evaluate a cubic polynomial given in Hermite form at an
% array of points. While designed for use by DPCHFE, it may
% be useful directly as an evaluator for a piecewise cubic
% Hermite function in applications, such as graphing, where
% the interval is known in advance.
%***LIBRARY SLATEC (PCHIP)
%***CATEGORY E3
%***TYPE doubleprecision (CHFEV-S, DCHFEV-D)
%***KEYWORDS CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,
% PCHIP
%***AUTHOR Fritsch, F. N., (LLNL)
% Lawrence Livermore National Laboratory
% P.O. Box 808 (L-316)
% Livermore, CA 94550
% FTS 532-4275, (510) 422-4275
%***DESCRIPTION
%
% DCHFEV: Cubic Hermite Function EValuator
%
% Evaluates the cubic polynomial determined by function values
% F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
% XE(J), J=1(1)NE.
%
% ----------------------------------------------------------------------
%
% Calling sequence:
%
% INTEGER NE, NEXT(2), IERR
% doubleprecision X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
%
% CALL DCHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
%
% Parameters:
%
% X1,X2 -- (input) endpoints of interval of definition of cubic.
% (Error return if X1.EQ.X2 .)
%
% F1,F2 -- (input) values of function at X1 and X2, respectively.
%
% D1,D2 -- (input) values of derivative at X1 and X2, respectively.
%
% NE -- (input) number of evaluation points. (Error return if
% NE.LT.1 .)
%
% XE -- (input) real*8 array of points at which the function is to
% be evaluated. If any of the XE are outside the interval
% [X1,X2], a warning error is returned in NEXT.
%
% FE -- (output) real*8 array of values of the cubic function
% defined by X1,X2, F1,F2, D1,D2 at the points XE.
%
% NEXT -- (output) integer array indicating number of extrapolation
% points:
% NEXT(1) = number of evaluation points to left of interval.
% NEXT(2) = number of evaluation points to right of interval.
%
% IERR -- (output) error flag.
% Normal return:
% IERR = 0 (no errors).
% 'Recoverable' errors:
% IERR = -1 if NE.LT.1 .
% IERR = -2 if X1.EQ.X2 .
% (The FE-array has not been changed in either case.)
%
%***REFERENCES (NONE)
%***ROUTINES CALLED XERMSG
%***REVISION HISTORY (YYMMDD)
% 811019 DATE WRITTEN
% 820803 Minor cosmetic changes for release 1.
% 870813 Corrected XERROR calls for d.p. names(s).
% 890206 Corrected XERROR calls.
% 890411 Added SAVE statements (Vers. 3.2).
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890703 Corrected category record. (WRB)
% 890831 Modified array declarations. (WRB)
% 891006 Cosmetic changes to prologue. (WRB)
% 891006 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
%***end PROLOGUE DCHFEV
% Programming notes:
%
% To produce a single precision version, simply:
% a. Change DCHFEV to CHFEV wherever it occurs,
% b. Change the doubleprecision declaration to real, and
% c. Change the constant ZERO to single precision.
%
% DECLARE ARGUMENTS.
%
persistent c2 c3 del1 del2 delta firstCall h i x xma xmi zero ; if isempty(firstCall),firstCall=1;end;
xe_shape=size(xe);xe=reshape(xe,1,[]);
fe_shape=size(fe);fe=reshape(fe,1,[]);
%
% DECLARE LOCAL VARIABLES.
%
if isempty(i), i=0; end;
if isempty(c2), c2=0; end;
if isempty(c3), c3=0; end;
if isempty(del1), del1=0; end;
if isempty(del2), del2=0; end;
if isempty(delta), delta=0; end;
if isempty(h), h=0; end;
if isempty(x), x=0; end;
if isempty(xmi), xmi=0; end;
if isempty(xma), xma=0; end;
if isempty(zero), zero=0; end;
if firstCall, zero=[0.0d0]; end;
firstCall=0;
%
% VALIDITY-CHECK ARGUMENTS.
%
%***FIRST EXECUTABLE STATEMENT DCHFEV
if( ne<1 )
%
% ERROR RETURNS.
%
% NE.LT.1 RETURN.
ierr = -1;
[dumvar1,dumvar2,dumvar3,ierr]=xermsg('SLATEC','DCHFEV','NUMBER OF EVALUATION POINTS LESS THAN ONE',ierr,1);
xe_shape=zeros(xe_shape);xe_shape(:)=xe(1:numel(xe_shape));xe=xe_shape;
fe_shape=zeros(fe_shape);fe_shape(:)=fe(1:numel(fe_shape));fe=fe_shape;
return;
else;
h = x2 - x1;
if( h==zero )
%
% X1.EQ.X2 RETURN.
ierr = -2;
[dumvar1,dumvar2,dumvar3,ierr]=xermsg('SLATEC','DCHFEV','INTERVAL ENDPOINTS EQUAL',ierr,1);
else;
%
% INITIALIZE.
%
ierr = 0;
next(1) = 0;
next(2) = 0;
xmi = min(zero,h);
xma = max(zero,h);
%
% COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
%
delta =(f2-f1)./h;
del1 =(d1-delta)./h;
del2 =(d2-delta)./h;
% (DELTA IS NO LONGER NEEDED.)
c2 = -(del1+del1+del2);
c3 =(del1+del2)./h;
% (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
%
% EVALUATION LOOP.
%
for i = 1 : ne;
x = xe(i) - x1;
fe(i) = f1 + x.*(d1+x.*(c2+x.*c3));
% COUNT EXTRAPOLATION POINTS.
if( x<xmi )
next(1) = fix(next(1) + 1);
end;
if( x>xma )
next(2) = fix(next(2) + 1);
end;
% (NOTE REDUNDANCY--IF EITHER CONDITION IS truemlv, OTHER IS falsemlv.)
end; i = fix(ne+1);
%
% NORMAL RETURN.
%
xe_shape=zeros(xe_shape);xe_shape(:)=xe(1:numel(xe_shape));xe=xe_shape;
fe_shape=zeros(fe_shape);fe_shape(:)=fe(1:numel(fe_shape));fe=fe_shape;
return;
end;
end;
%------------- LAST LINE OF DCHFEV FOLLOWS -----------------------------
xe_shape=zeros(xe_shape);xe_shape(:)=xe(1:numel(xe_shape));xe=xe_shape;
fe_shape=zeros(fe_shape);fe_shape(:)=fe(1:numel(fe_shape));fe=fe_shape;
end
%DECK DCHFIE
|
|
