ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
helena_spline.f90
Go to the documentation of this file.
1 module helena_spline
2 
3  use itm_types
4 
5  implicit none
6 
7  type spline_coefficients
8  sequence
9  real(r8), dimension(:), pointer :: sp1, sp2, sp3, sp4
10  end type spline_coefficients
11 
12 
13 contains
14 
15 subroutine spline(n, x, y, alfa, beta, typ, p)
16 !-----------------------------------------------------------------------------
17 ! input:
18 !
19 ! n Anzahl der Knoten
20 ! x Array der x-Werte
21 ! y Array der y-Werte
22 ! alfa Randbedingung in x(1)
23 ! beta " in x(n)
24 ! typ = 0 not-a-knot Spline
25 ! 1 alfa, beta 1. Ableitungen vorgegeben
26 ! 2 " " 2. " "
27 ! 3 " " 3. " "
28 !
29 ! Bemerkung: mit typ = 2 und alfa = beta = 0 erhaelt man
30 ! einen natuerlichen Spline
31 !
32 ! output:
33 !
34 ! p Arrays der Splinekoeffizienten
35 ! s = p%sp1(i) + p%sp2(i) * (x - x(i)) + p%sp3(i) * (x - x(i))**2 &
36 ! + p%sp4(i) * (x - x(i))**3
37 !
38 ! Bei Anwendungsfehlern wird das Programm mit entsprechender
39 ! Fehlermeldung abgebrochen
40 !
41 ! recent extension: spline now works with both monotonically increasing
42 ! and decreasing x values
43 ! simply use spline without worrying about changes in
44 ! the derivatives
45 ! sign changes to the specified derivatives are done
46 ! automatically inside spline and undone on return
47 ! warning: must be used with spwert and same interpolation basis
48 !-----------------------------------------------------------------------------
49 
50  use itm_types
51 
52  implicit none
53 
54  integer(itm_i4), intent(in) :: n, typ
55  real(r8), dimension(n), intent(in) :: x
56  real(r8), dimension(n), intent(in) :: y
57  real(r8), intent(in) :: alfa, beta
58  type(spline_coefficients), intent(out) :: p
59 
60  real(r8), dimension(n) :: xlocal
61  integer(itm_i4), parameter :: iu6 = 6
62 
63  real(r8) :: h(10001)
64  real(r8) :: alfa_loc, beta_loc
65  integer(itm_i4) :: i, ierr, nmax
66  logical :: inverted
67 
68  xlocal = x
69 !-- check whether abscissa has to be inverted
70  if (xlocal(n) < xlocal(1)) then
71  inverted = .true.
72  else
73  inverted = .false.
74  end if
75 
76  nmax = 10001
77  if (typ < 0 .or. typ > 3) then
78  write(iu6, *) 'Error in routine spline: wrong type'
79  stop
80  end if
81 
82  if (n < 3 .or. n > nmax) then
83  write(iu6, *) 'Error in routine spline: n < 3 or n > nmax'
84  write(iu6, *) 'n = ',n,' , nmax = ', nmax
85  stop
86  end if
87 
88 !-- invert abscissa if necessary and adapt derivatives
89  if (inverted) then
90  xlocal = -xlocal
91  alfa_loc = (-1)**typ * alfa
92  beta_loc = (-1)**typ * beta
93  else
94  alfa_loc = alfa
95  beta_loc = beta
96  end if
97 
98 !-- check whether abscissa is monotonic
99  do i = 1, n - 1
100  h(i) = xlocal(i + 1) - xlocal(i)
101  if (h(i) <= 0._r8) then
102  write(iu6, *) 'Error in spline: abscissa not monotonic x(i - 1) >= x(i)'
103  write(*,*) i
104  write(*,*) xlocal(1:n)
105  pause 'DPC: Error in spline: abscissa not monotonic x(i - 1) >= x(i)'
106  stop
107  end if
108  end do
109 
110 !-- build equations
111  do i = 1, n - 2
112  p%sp1(i) = 3.0_r8 * ((y(i + 2) - y(i + 1)) / h(i + 1) - (y(i + 1) &
113  - y(i)) / h(i))
114  p%sp2(i) = h(i)
115  p%sp3(i) = h(i + 1)
116  p%sp4(i) = 2.0_r8 * (h(i) + h(i + 1))
117  end do
118 
119 !-- boundary conditions
120 
121  select case (typ)
122 !-- not-a-knot spline
123  case (0)
124  p%sp1(1) = p%sp1(1) * h(2) / (h(1) + h(2))
125  p%sp1(n - 2) = p%sp1(n - 2) * h(n - 2) / (h(n - 1) + h(n - 2))
126  p%sp4(1) = p%sp4(1) - h(1)
127  p%sp4(n - 2) = p%sp4(n - 2) - h(n - 1)
128  p%sp3(1) = p%sp3(1) - h(1)
129  p%sp2(n - 2) = p%sp2(n - 2) - h(n - 1)
130 !-- first derivative given
131  case (1)
132  p%sp1(1) = p%sp1(1) - 1.5_r8 * ((y(2) - y(1)) / h(1) - alfa_loc)
133  p%sp1(n - 2) = p%sp1(n - 2) - 1.5_r8 * (beta_loc - (y(n) - y(n - 1)) &
134  / h(n - 1))
135  p%sp4(1) = p%sp4(1) - 0.5_r8 * h(1)
136  p%sp4(n - 2) = p%sp4(n - 2) - 0.5_r8 * h(n - 1)
137 !-- second derivative given
138  case (2)
139  p%sp1(1) = p%sp1(1) - 0.5_r8 * alfa_loc * h(1)
140  p%sp1(n - 2) = p%sp1(n - 2) - 0.5_r8 * beta_loc * h(n - 1)
141 !-- third derivative given
142  case (3)
143  p%sp1(1) = p%sp1(1) + 0.5_r8 * alfa_loc * h(1) * h(1)
144  p%sp1(n - 2) = p%sp1(n - 2) - 0.5_r8 * beta_loc * h(n - 1) * h(n - 1)
145  p%sp4(1) = p%sp4(1) + h(1)
146  p%sp4(n - 2) = p%sp4(n - 2) + h(n - 1)
147  end select
148 
149 !-- calculate coefficients
150  call sgtsl(n - 2, p%sp2, p%sp4, p%sp3, p%sp1, ierr)
151 
152  if (ierr /= 0) then
153  write(iu6, *) 'error in sgtsl: matrix singular'
154  stop
155  end if
156 
157 !-- write solution vector
158  call dcopy(n - 2, p%sp1, 1, p%sp3(2), 1)
159 
160 !-- correct first and last value of p%sp3 depending on boundary conditions
161  select case (typ)
162  case (0)
163  p%sp3(1) = p%sp3(2) + h(1) * (p%sp3(2) - p%sp3(3)) / h(2)
164  p%sp3(n) = p%sp3(n - 1) + h(n - 1) * (p%sp3(n - 1) - p%sp3(n - 2)) &
165  / h(n - 2)
166  case (1)
167  p%sp3(1) = 1.5_r8*((y(2) - y(1)) / h(1) - alfa_loc) / h(1) - 0.5_r8 &
168  * p%sp3(2)
169  p%sp3(n) = -1.5_r8*((y(n) - y(n - 1)) / h(n - 1) - beta_loc) / &
170  h(n - 1) - 0.5_r8 * p%sp3(n - 1)
171  case (2)
172  p%sp3(1) = 0.5_r8 * alfa_loc
173  p%sp3(n) = 0.5_r8 * beta_loc
174  case (3)
175  p%sp3(1) = p%sp3(2) - 0.5_r8 * alfa_loc * h(1)
176  p%sp3(n) = p%sp3(n - 1) + 0.5_r8 * beta_loc * h(n - 1)
177  end select
178  call dcopy(n, y, 1, p%sp1, 1)
179 
180  do i = 1, n - 1
181  p%sp2(i) = (p%sp1(i + 1) - p%sp1(i)) / h(i) - h(i) * (p%sp3(i + 1) &
182  + 2.0_r8 * p%sp3(i)) / 3.0_r8
183  p%sp4(i) = (p%sp3(i + 1) - p%sp3(i)) / (3.0_r8 * h(i))
184  end do
185 
186  p%sp2(n) = (3.0_r8 * p%sp4(n - 1) * h(n - 1) + 2.0_r8 &
187  * p%sp3(n - 1)) * h(n - 1) + p%sp2(n - 1)
188 
189  p%sp4(n) = 0.0_r8 ! for complete initialization
190 
191 !-- undo changes to abscissa and derivatives
192  if (inverted) then
193  xlocal = -xlocal
194  end if
195 
196  return
197 
198 end subroutine spline
199 
200 function spwert(n, xwert, p, x, abltg, order) result(f_spwert)
201 !-----------------------------------------------------------------------
202 ! input:
203 !
204 ! n anzahl der knotenpunkte
205 ! xwert stelle an der funktionswerte berechnet werden
206 ! p arrays der splinekoeffizienten (aus spline)
207 ! x array der knotenpunkte
208 !
209 ! output:
210 !
211 ! spwert funktionswert an der stelle xwert
212 ! abltg(i) wert der i-ten ableitung bei xwert
213 ! order specifies the maximum order of returned derivatives
214 !-----------------------------------------------------------------------
215 
216  use itm_types
217 
218  implicit none
219 
220  integer(itm_i4), intent(in) :: n
221  type(spline_coefficients), intent(in) :: p
222  real(r8), intent(in) :: xwert
223  real(r8), intent(in) :: x(n)
224  integer(itm_i4), intent(in) :: order
225  real(r8), intent(out) :: abltg(3)
226  real(r8) :: f_spwert
227  integer(itm_i4), save :: i, k
228  integer(itm_i4) :: m
229 
230  real(r8) :: x_loc(n)
231  real(r8) :: xx, xwert_loc
232  logical, save :: first = .true.
233  logical :: inverted
234 
235 !-- check whether abscissa has to be inverted
236  if (x(n) < x(1)) then
237  inverted = .true.
238  else
239  inverted = .false.
240  end if
241 
242 !-- invert abscissa if necessary
243  if (inverted) then
244  x_loc = -x
245  xwert_loc = -xwert
246  else
247  x_loc = x
248  xwert_loc = xwert
249  end if
250 
251 !-- find matching interval (binary search)
252  if (first .or. i > n .or. k > n) then
253  i = 1
254  k = n
255  first = .false.
256  else
257  if (xwert_loc < x_loc(k)) then
258  i = 1
259  else
260  k = n
261  end if
262  end if
263 
264  m = (i + k) / 2
265 
266  do
267  if (m == i) exit
268  if (xwert_loc < x_loc(m)) then
269  k = m
270  else
271  i = m
272  end if
273  m = (i + k) / 2
274  end do
275 
276  xx = xwert_loc - x_loc(i)
277  f_spwert = ((p%sp4(i) * xx + p%sp3(i)) * xx + p%sp2(i)) * xx + p%sp1(i)
278 
279  if (order <= 0) return
280 
281  abltg(1) = (3.0_r8 * p%sp4(i) * xx + 2.0_r8 * p%sp3(i)) * xx &
282  + p%sp2(i)
283 
284 !-- adapt derivative if necessary
285  if (inverted) then
286  abltg(1) = -abltg(1)
287  end if
288 
289  if (order == 1) return
290 
291  abltg(2) = 6.0_r8 * p%sp4(i) * xx + 2.0_r8 * p%sp3(i)
292 
293  if (order == 2) return
294 
295  abltg(3) = 6.0_r8 * p%sp4(i)
296 
297 !-- adapt derivative if necessary
298  if (inverted) then
299  abltg(3) = -abltg(3)
300  end if
301 
302  return
303 
304 end function spwert
305 
306 end module helena_spline
spline
subroutine spline(N, X, Y, ALFA, BETA, TYP, A, B, C, D)
Definition: solution6.f90:655
spwert
REAL *8 function spwert(N, XWERT, A, B, C, D, X, ABLTG)
Definition: solution6.f90:155
sgtsl
subroutine sgtsl(n, c, d, e, b, info)
Definition: sgtsl.f90:3
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
helena_spline::spline_coefficients
Definition: helena_spline.f90:7
helena_spline
Definition: helena_spline.f90:1