forked from MDSplus/mdsplus
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ppval.f
197 lines (196 loc) · 7.59 KB
/
ppval.f
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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
! Ken Klare, LANL P-4 (c)1990,1991
!***********************************************************************
!********************* 3.5 Piecewise Polynomial ************************
!-----------------------------------------------------------------------
Real Function PPVAL(x,korder,nintv,break,c)
Real x !input, the point to evaluate
Integer korder !input, the order of the polynomial (degree+1)
Integer nintv !input, the number of polynomial pieces
Real break(*)!input, [nintv+1] breakpoints
Real c(*) !input, [korder,nintv] polynomial coefficients
!+Evaluate a piecewise polynomial
Real PPDER
PPVAL = PPDER(0,x,korder,nintv,break,c)
End
!-----------------------------------------------------------------------
Real Function PPDER(ideriv,x,korder,nintv,break,c)
Integer ideriv !input, the order of deriavtive to be evaluated
Real x !input, the point to evaluate
Integer korder !input, the order of the polynomial (degree+1)
Integer nintv !input, the number of polynomial pieces
Real break(*)!input, [nintv+1] breakpoints
Real c(korder,*)!input, [korder,nintv] polynomial coefficients
!+Evaluate derivative of a piecewise polynomial
Integer j,m,keep
Real h
Data keep/1000000/
! can only take so many derivatives
If (ideriv.ge.korder) Then
PPDER = 0
Else
! look up the place
Call SRCH(nintv,x,break,1,keep)
If (keep.gt.0) Then
PPDER = c(ideriv+1,keep)
Else
j = -1-keep
if (j.le.0) j = 1
h = x - break(j)
PPDER = c(korder,j)
Do m = korder-ideriv-1, 1, -1
PPDER = PPDER*h/m + c(m+ideriv,j)
Enddo
Endif
Endif
End
!-----------------------------------------------------------------------
Real Function PPITG(a,b,korder,nintv,break,c)
Real a,b !input, the range to evaluate
Integer korder !input, the order of the polynomial (degree+1)
Integer nintv !input, the number of polynomial pieces
Real break(*)!input, [nintv] breakpoints
Real c(korder,*)!input, [korder,nintv] polynomial coefficients
!+Evaluate the integral of a piecewise polynomial
Integer ia, ib, j
Real P9ITG
data ia /999999/
data ib /999999/
Call SRCH(nintv,b,break,1,ib)
If (ib.gt.0) Then
PPITG = 0
Else
ib = -1-ib
if (ib.le.0) ib = 1
PPITG = P9ITG(b-break(ib),korder,c(1,ib))
Endif
Call SRCH(nintv,a,break,1,ia)
If (ia.le.0) Then
ia = -1-ia
if (ia.le.0) ia = 1
PPITG = PPITG - P9ITG(a-break(ia),korder,c(1,ia))
Endif
Do j = ia,ib-1
PPITG = PPITG + P9ITG(break(j+1)-break(j),korder,c(1,j))
Enddo
Do j = ib,ia-1
PPITG = PPITG - P9ITG(break(j+1)-break(j),korder,c(1,j))
Enddo
End
!-----------------------------------------------------------------------
Real Function P9ITG(h,korder,c)
Real h !input, offset point
Integer korder !input, number of terms, degree+1
Real c(*) !input, [korder] value and derivatives
Integer j
!+Evaluate integral of single polynomial
P9ITG = 0
Do j = korder,1,-1
P9ITG = (P9ITG*h + c(j))/j
Enddo
P9ITG = P9ITG*h
End
!-----------------------------------------------------------------------
Subroutine PCOEF(ndata,xdata,fdata,cof)
Integer ndata !input, polynomial order (degree+1)
Real xdata(*)!input, [ndata] the ordinates, must be distinct
Real fdata(*)!input, [ndata] the abscissae
Real cof(*) !input, [ndata] the polynomial coeff
!+Create single piecewise polynomial, value and derivatives at first
! Not an IMSL routine. Method of G.Rybicki from Numerical Recipes POLCOE
Integer nmax
Parameter (nmax=15)
Integer i,j,k,n
Real xx,s(nmax),b,ff,phi,fact(3:nmax)
Data fact/2,6,24,120,720,5040,40320,362880,3628800,39916800,
1 479001600,6227020800E0,87178291200E0/
n=ndata
Do i=1,n
cof(i)=0
s(i)=0
Enddo
! s(n)=-(xdata(1)-xdata(1))
! recurrence
Do i=2,n
xx=xdata(i)-xdata(1)
Do j=n+1-i,n-1
s(j)=s(j)-xx*s(j+1)
Enddo
s(n)=s(n)-xx
Enddo
Do j=1,n
xx=xdata(j)-xdata(1)
phi=n
! phi=product(xj-xk) for j.ne.k
Do k=n-1,1,-1
phi=k*s(k+1)+xx*phi
Enddo
ff=fdata(j)/phi
b=1
Do k=n,1,-1
cof(k)=cof(k)+b*ff
b=s(k)+xx*b
Enddo
Enddo
Do k=3,n
cof(k)=cof(k)*fact(k)
Enddo
End
!***********************************************************************
!********************* Extension for vectors ***************************
!-----------------------------------------------------------------------
Subroutine PPVALV(nx,x,korder,nintv,break,c,y)
Integer nx !input, the number of points
Real x(*) !input, the [nx] points to evaluate
Integer korder !input, the order of the polynomial (degree+1)
Integer nintv !input, the number of polynomial pieces
Real break(*)!input, [nintv+1] breakpoints
Real c(*) !input, [korder,nintv] polynomial coefficients
Real y(*) !output, the evaluated points
!+Evaluate a piecewise polynomial
Call PPDERV(0,nx,x,korder,nintv,break,c,y)
End
!-----------------------------------------------------------------------
Subroutine PPDERV(ideriv,nx,x,korder,nintv,break,c,y)
Integer ideriv !input, the order of deriavtive to be evaluated
Integer nx !input, the number of points
Real x(*) !input, the [nx] points to evaluate
Integer korder !input, the order of the polynomial (degree+1)
Integer nintv !input, the number of polynomial pieces
Real break(*)!input, [nintv+1] breakpoints
Real c(korder,*)!input, [korder,nintv] polynomial coefficients
Real y(*) !output, the evaluated points
!+Evaluate derivative of a piecewise polynomial
Integer j
Real PPDER
Do j=1,nx
y(j) = PPDER(ideriv,x(j),korder,nintv,break,c)
Enddo
End
!-----------------------------------------------------------------------
Subroutine PPITGV(na,a,nb,b,korder,nintv,break,c,y)
Integer na !input, the number of points
Real a(*) !input, the [na] points to evaluate from first
Integer nb !input, the number of points
Real b(*) !input, the [nb] points to evaluate from first
Integer korder !input, the order of the polynomial (degree+1)
Integer nintv !input, the number of polynomial pieces
Real break(*)!input, [nintv] breakpoints
Real c(korder,*)!input, [korder,nintv] polynomial coefficients
Real y(*) !output, the evaluated points
!+Evaluate the integral of a piecewise polynomial
Real PPITG
Integer j,ja,jb,n,inca,incb
n = min(na,nb)
If (na.eq.1) n = nb
inca = 0
incb = 0
If (na.gt.1) inca = 1
If (nb.gt.1) incb = 1
ja = 1
jb = 1
Do j=1,n
y(j) = PPITG(a(ja),b(jb),korder,nintv,break,c)
ja = ja + inca
jb = jb + incb
Enddo
End