1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
*DECK LIBA28
SUBROUTINE LIBA28(X,XP,NXP,LL,FP,L,IPROX,I0)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the interpolation weights FP for the Lagrange interpolation.
*
*Copyright:
* Copyright (C) 2002 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input
* X abscissa.
* XP tabulated abscissa.
* NXP number of tabulated points.
* LL requested interpolation order.
*
*Parameters: output
* FP weights for Lagrange interpolation.
* L interpolation limit (number of non-zero weights).
* IPROX index of closest tabulated point.
* I0 number of leading zero weights.
*
*Comments:
* Evaluation method.
* F(X) = sum for I = 1 to L of F(I+I0)*FP(I)
* for LL.le.NXP uses an LL-point "centered" Lagrange interpolation,
* otherwise it uses a linear interpolation formula.
* Attention: it is assumed that XP(I+1) > XP(I) for all I
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NXP,LL,L,IPROX,I0
REAL X
DOUBLE PRECISION XP(NXP),FP(LL)
*----
* LOCAL VARIABLES
*----
DOUBLE PRECISION PP,XI,XX
*
L=LL
IF(L.LE.NXP) GO TO 20
*----
* CONSTANT FUNCTION DEFINED BY A SINGLE POINT
*----
IF(NXP.LE.1)THEN
I0=0
L=1
IPROX=1
FP(1)=1.0D0
GO TO 100
ENDIF
*----
* NXP < L: SWITCH TO LINEAR INTERPOLATION
*----
L=2
20 L2=(L+1)/2
*----
* LOCATE FIRST POINT TO THE RIGHT OF XX (PT IXP)
*----
DO 30 IXP=1,NXP
IF(XP(IXP).GT.X)GO TO 40
30 CONTINUE
*----
* X IS TO THE RIGHT OF EVERY POINT XP
*----
IXP=1
IPROX=NXP
GO TO 60
*----
* XP(IMIN) IS THE FIRST POINT FOR THE INTERPOLATION
*----
40 IMIN=IXP-L2
IPROX=IXP
IF(IXP.GT.1)THEN
IF((X-XP(IXP-1)).LT.(XP(IXP)-X))IPROX=IXP-1
ENDIF
IF(L.EQ.1)THEN
FP(1)=1.0D0
I0=IPROX-1
GO TO 100
ENDIF
IF(IMIN.GE.1)GO TO 50
IMIN=1
IMAX=L
GO TO 70
*----
* XP(IMAX) IS THE LAST POINT FOR THE INTERPOLATION
*----
50 IMAX=IMIN+L-1
IF(IMAX.LE.NXP)GO TO 70
60 IMAX=NXP
IMIN=NXP-L+1
*----
* CENTERED POLYNOMIAL INTERPOLATION OF DEGRE L
*----
70 I0=IMIN-1
XX=X
DO 90 I=IMIN,IMAX
PP=1.0D0
XI=XP(I)
DO 80 J=IMIN,IMAX
IF(I.NE.J)PP=PP*((XX-XP(J))/(XI-XP(J)))
80 CONTINUE
FP(I+1-IMIN)=PP
90 CONTINUE
100 RETURN
END
|
