source: LMDZ5/trunk/libf/misc/chfev.F @ 3719

Last change on this file since 3719 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 4.9 KB
Line 
1*DECK CHFEV
2      SUBROUTINE CHFEV (X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)
3C***BEGIN PROLOGUE  CHFEV
4C***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an
5C            array of points.  While designed for use by PCHFE, it may
6C            be useful directly as an evaluator for a piecewise cubic
7C            Hermite function in applications, such as graphing, where
8C            the interval is known in advance.
9C***LIBRARY   SLATEC (PCHIP)
10C***CATEGORY  E3
11C***TYPE      SINGLE PRECISION (CHFEV-S, DCHFEV-D)
12C***KEYWORDS  CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,
13C             PCHIP
14C***AUTHOR  Fritsch, F. N., (LLNL)
15C             Lawrence Livermore National Laboratory
16C             P.O. Box 808  (L-316)
17C             Livermore, CA  94550
18C             FTS 532-4275, (510) 422-4275
19C***DESCRIPTION
20C
21C          CHFEV:  Cubic Hermite Function EValuator
22C
23C     Evaluates the cubic polynomial determined by function values
24C     F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
25C     XE(J), J=1(1)NE.
26C
27C ----------------------------------------------------------------------
28C
29C  Calling sequence:
30C
31C        INTEGER  NE, NEXT(2), IERR
32C        REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
33C
34C        CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
35C
36C   Parameters:
37C
38C     X1,X2 -- (input) endpoints of interval of definition of cubic.
39C           (Error return if  X1.EQ.X2 .)
40C
41C     F1,F2 -- (input) values of function at X1 and X2, respectively.
42C
43C     D1,D2 -- (input) values of derivative at X1 and X2, respectively.
44C
45C     NE -- (input) number of evaluation points.  (Error return if
46C           NE.LT.1 .)
47C
48C     XE -- (input) real array of points at which the function is to be
49C           evaluated.  If any of the XE are outside the interval
50C           [X1,X2], a warning error is returned in NEXT.
51C
52C     FE -- (output) real array of values of the cubic function defined
53C           by  X1,X2, F1,F2, D1,D2  at the points  XE.
54C
55C     NEXT -- (output) integer array indicating number of extrapolation
56C           points:
57C            NEXT(1) = number of evaluation points to left of interval.
58C            NEXT(2) = number of evaluation points to right of interval.
59C
60C     IERR -- (output) error flag.
61C           Normal return:
62C              IERR = 0  (no errors).
63C           "Recoverable" errors:
64C              IERR = -1  if NE.LT.1 .
65C              IERR = -2  if X1.EQ.X2 .
66C                (The FE-array has not been changed in either case.)
67C
68C***REFERENCES  (NONE)
69C***ROUTINES CALLED  XERMSG
70C***REVISION HISTORY  (YYMMDD)
71C   811019  DATE WRITTEN
72C   820803  Minor cosmetic changes for release 1.
73C   890411  Added SAVE statements (Vers. 3.2).
74C   890531  Changed all specific intrinsics to generic.  (WRB)
75C   890703  Corrected category record.  (WRB)
76C   890703  REVISION DATE from Version 3.2
77C   891214  Prologue converted to Version 4.0 format.  (BAB)
78C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
79C***END PROLOGUE  CHFEV
80C  Programming notes:
81C
82C     To produce a double precision version, simply:
83C        a. Change CHFEV to DCHFEV wherever it occurs,
84C        b. Change the real declaration to double precision, and
85C        c. Change the constant ZERO to double precision.
86C
87C  DECLARE ARGUMENTS.
88C
89      INTEGER  NE, NEXT(2), IERR
90      REAL  X1, X2, F1, F2, D1, D2, XE(*), FE(*)
91C
92C  DECLARE LOCAL VARIABLES.
93C
94      INTEGER  I
95      REAL  C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
96      SAVE ZERO
97      DATA  ZERO /0./
98C
99C  VALIDITY-CHECK ARGUMENTS.
100C
101C***FIRST EXECUTABLE STATEMENT  CHFEV
102      IF (NE .LT. 1)  GO TO 5001
103      H = X2 - X1
104      IF (H .EQ. ZERO)  GO TO 5002
105C
106C  INITIALIZE.
107C
108      IERR = 0
109      NEXT(1) = 0
110      NEXT(2) = 0
111      XMI = MIN(ZERO, H)
112      XMA = MAX(ZERO, H)
113C
114C  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
115C
116      DELTA = (F2 - F1)/H
117      DEL1 = (D1 - DELTA)/H
118      DEL2 = (D2 - DELTA)/H
119C                                           (DELTA IS NO LONGER NEEDED.)
120      C2 = -(DEL1+DEL1 + DEL2)
121      C3 = (DEL1 + DEL2)/H
122C                               (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
123C
124C  EVALUATION LOOP.
125C
126      DO 500  I = 1, NE
127         X = XE(I) - X1
128         FE(I) = F1 + X*(D1 + X*(C2 + X*C3))
129C          COUNT EXTRAPOLATION POINTS.
130         IF ( X.LT.XMI )  NEXT(1) = NEXT(1) + 1
131         IF ( X.GT.XMA )  NEXT(2) = NEXT(2) + 1
132C        (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
133  500 CONTINUE
134C
135C  NORMAL RETURN.
136C
137      RETURN
138C
139C  ERROR RETURNS.
140C
141 5001 CONTINUE
142C     NE.LT.1 RETURN.
143      IERR = -1
144      CALL XERMSG ('SLATEC', 'CHFEV',
145     +   'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
146      RETURN
147C
148 5002 CONTINUE
149C     X1.EQ.X2 RETURN.
150      IERR = -2
151      CALL XERMSG ('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', IERR,
152     +   1)
153      RETURN
154C------------- LAST LINE OF CHFEV FOLLOWS ------------------------------
155      END
Note: See TracBrowser for help on using the repository browser.