EFTCAMB  Reference documentation for version 3.0
02_equispaced_interpolation_linear_1D.f90
Go to the documentation of this file.
1 !----------------------------------------------------------------------------------------
2 !
3 ! This file is part of EFTCAMB.
4 !
5 ! Copyright (C) 2013-2016 by the EFTCAMB authors
6 !
7 ! The EFTCAMB code is free software;
8 ! You can use it, redistribute it, and/or modify it under the terms
9 ! of the GNU General Public License as published by the Free Software Foundation;
10 ! either version 3 of the License, or (at your option) any later version.
11 ! The full text of the license can be found in the file eftcamb/LICENSE at
12 ! the top level of the EFTCAMB distribution.
13 !
14 !----------------------------------------------------------------------------------------
15 
19 
20 
21 !----------------------------------------------------------------------------------------
24 
26 
28 
29  use precision
30  use amlutils
32 
33  implicit none
34 
35  private
36 
37  public equispaced_linear_interpolate_function_1d
38 
39  !----------------------------------------------------------------------------------------
41  type :: equispaced_linear_interpolate_function_1d
42 
43  ! options:
44  integer :: num_points
45 
46  ! parameters:
47  real(dl) :: x_initial
48  real(dl) :: x_final
49  real(dl) :: null_value
50  logical :: has_null_value
51  real(dl) :: grid_width
52 
53  ! arrays with the values:
54  real(dl), allocatable, dimension(:) :: x
55  real(dl), allocatable, dimension(:) :: y
56  real(dl), allocatable, dimension(:) :: yp
57  real(dl), allocatable, dimension(:) :: ypp
58  real(dl), allocatable, dimension(:) :: yppp
59  real(dl), allocatable, dimension(:) :: yint
60 
61  contains
62 
63  procedure :: initialize => equispacedlinearintepolatefunction1dinit
64  procedure :: precompute => equispacedlinearintepolatefunction1dprecompute
65  procedure :: value => equispacedlinearintepolatefunction1dvalue
66  procedure :: first_derivative => equispacedlinearintepolatefunction1dfirstderivative
67  procedure :: second_derivative => equispacedlinearintepolatefunction1dsecondderivative
68  procedure :: third_derivative => equispacedlinearintepolatefunction1dthirdderivative
69  procedure :: integral => equispacedlinearintepolatefunction1dintegral
70  procedure :: initialize_derivatives => equispacedlinearintepolatefunction1dinitderivatives
71 
72  end type equispaced_linear_interpolate_function_1d
73 
74  ! ---------------------------------------------------------------------------------------------
75 
76 contains
77 
78  ! ---------------------------------------------------------------------------------------------
80  subroutine equispacedlinearintepolatefunction1dinit( self, num_points, x_initial, x_final, null_value )
81 
82  implicit none
83 
84  class(equispaced_linear_interpolate_function_1d) :: self
85 
86  integer , intent(in) :: num_points
87  real(dl), intent(in) :: x_initial
88  real(dl), intent(in) :: x_final
89  real(dl), intent(in), optional :: null_value
90 
91  integer :: i
92 
93  ! initialize the null value:
94  if ( present(null_value) ) then
95  self%has_null_value = .true.
96  self%null_value = null_value
97  else
98  self%has_null_value = .false.
99  self%null_value = 0._dl
100  end if
101 
102  ! allocate the x vector:
103  if ( allocated(self%x) ) deallocate( self%x )
104 
105  self%num_points = num_points
106 
107  allocate( self%x( self%num_points ) )
108 
109  ! store initial and final times:
110  self%x_initial = x_initial
111  self%x_final = x_final
112 
113  ! compute the width of the interpolation:
114  self%grid_width = ( self%x_final -self%x_initial )/REAL( self%num_points -1 )
115 
116  ! fill in the equispaced x array:
117  do i=1, self%num_points
118  self%x(i) = self%x_initial + REAL(i-1)*self%grid_width
119  end do
120 
121  ! allocate the other vectors:
122  if ( allocated(self%y) ) deallocate( self%y )
123  if ( allocated(self%yp) ) deallocate( self%yp )
124  if ( allocated(self%ypp) ) deallocate( self%ypp )
125  if ( allocated(self%yppp) ) deallocate( self%yppp )
126  if ( allocated(self%yint) ) deallocate( self%yint )
127 
128  allocate( self%y( self%num_points ) )
129  allocate( self%yp( self%num_points ) )
130  allocate( self%ypp( self%num_points ) )
131  allocate( self%yppp( self%num_points ) )
132  allocate( self%yint( self%num_points ) )
133 
134  end subroutine equispacedlinearintepolatefunction1dinit
135 
136  ! ---------------------------------------------------------------------------------------------
139  subroutine equispacedlinearintepolatefunction1dprecompute( self, x, ind, mu )
140 
141  implicit none
142 
143  class(equispaced_linear_interpolate_function_1d) :: self
144  real(dl), intent(in) :: x
145  integer , intent(out) :: ind
146  real(dl), intent(out) :: mu
147 
148  real(dl) :: x1, x2
149 
150  ! compute the interpolation index:
151  ind = int( ( x-self%x_initial)/self%grid_width ) +1
152  ! store the x values:
153  x1 = self%x(ind)
154  x2 = self%x(ind+1)
155  ! compute the linear interpolation coefficient:
156  mu = (x-x1)/(x2-x1)
157 
158  end subroutine equispacedlinearintepolatefunction1dprecompute
159 
160  ! ---------------------------------------------------------------------------------------------
162  function equispacedlinearintepolatefunction1dvalue( self, x, index, coeff )
163 
164  implicit none
165 
166  class(equispaced_linear_interpolate_function_1d) :: self
167  real(dl), intent(in) :: x
168  integer , intent(in), optional :: index
169  real(dl), intent(in), optional :: coeff
170  real(dl) :: equispacedlinearintepolatefunction1dvalue
171  integer :: ind
172  real(dl) :: x1, x2, y1, y2, mu
173 
174  ! initialize to null value:
175  equispacedlinearintepolatefunction1dvalue = self%null_value
176  if ( self%has_null_value ) then
177  ! if outside the interpolation range return the null value:
178  if ( x <= self%x_initial .or. x >= self%x_final ) return
179  else
180  ! if below the interpolation range return the first value:
181  if ( x <= self%x_initial ) then
182  equispacedlinearintepolatefunction1dvalue = self%y(1)
183  return
184  end if
185  ! if above the interpolation range return the first value:
186  if ( x >= self%x_final ) then
187  equispacedlinearintepolatefunction1dvalue = self%y(self%num_points)
188  return
189  end if
190  end if
191 
192  ! return the index of the point:
193  if ( present(index) ) then
194  ind = index
195  else
196  ind = int( ( x-self%x_initial)/self%grid_width ) +1
197  end if
198 
199  ! get the interpolation coefficient:
200  if ( present(coeff) ) then
201  mu = coeff
202  else
203  ! store the x values:
204  x1 = self%x(ind)
205  x2 = self%x(ind+1)
206  ! compute the linear interpolation coefficient:
207  mu = (x-x1)/(x2-x1)
208  end if
209 
210  ! store the y values:
211  y1 = self%y(ind)
212  y2 = self%y(ind+1)
213 
214  ! compute the linear interpolation:
215  equispacedlinearintepolatefunction1dvalue = y1*( 1._dl -mu ) +y2*mu
216 
217  end function equispacedlinearintepolatefunction1dvalue
218 
219  ! ---------------------------------------------------------------------------------------------
221  function equispacedlinearintepolatefunction1dfirstderivative( self, x, index, coeff )
222 
223  implicit none
224 
225  class(equispaced_linear_interpolate_function_1d) :: self
226  real(dl), intent(in) :: x
227  integer , intent(in), optional :: index
228  real(dl), intent(in), optional :: coeff
229  real(dl) :: equispacedlinearintepolatefunction1dfirstderivative
230 
231  integer :: ind
232  real(dl) :: x1, x2, y1, y2, mu
233 
234  ! initialize to null value:
235  equispacedlinearintepolatefunction1dfirstderivative = self%null_value
236  if ( self%has_null_value ) then
237  ! if outside the interpolation range return the null value:
238  if ( x <= self%x_initial .or. x >= self%x_final ) return
239  else
240  ! if below the interpolation range return the first value:
241  if ( x <= self%x_initial ) then
242  equispacedlinearintepolatefunction1dfirstderivative = self%yp(1)
243  return
244  end if
245  ! if above the interpolation range return the first value:
246  if ( x >= self%x_final ) then
247  equispacedlinearintepolatefunction1dfirstderivative = self%yp(self%num_points)
248  return
249  end if
250  end if
251 
252  ! return the index of the point:
253  if ( present(index) ) then
254  ind = index
255  else
256  ind = int( ( x-self%x_initial)/self%grid_width ) +1
257  end if
258 
259  ! get the interpolation coefficient:
260  if ( present(coeff) ) then
261  mu = coeff
262  else
263  ! store the x values:
264  x1 = self%x(ind)
265  x2 = self%x(ind+1)
266  ! compute the linear interpolation coefficient:
267  mu = (x-x1)/(x2-x1)
268  end if
269 
270  ! store the y values:
271  y1 = self%yp(ind)
272  y2 = self%yp(ind+1)
273 
274  ! compute the linear interpolation:
275  equispacedlinearintepolatefunction1dfirstderivative = y1*( 1._dl -mu ) +y2*mu
276 
277  end function equispacedlinearintepolatefunction1dfirstderivative
278 
279  ! ---------------------------------------------------------------------------------------------
281  function equispacedlinearintepolatefunction1dsecondderivative( self, x, index, coeff )
282 
283  implicit none
284 
285  class(equispaced_linear_interpolate_function_1d) :: self
286  real(dl), intent(in) :: x
287  integer , intent(in), optional :: index
288  real(dl), intent(in), optional :: coeff
289  real(dl) :: equispacedlinearintepolatefunction1dsecondderivative
290 
291  integer :: ind
292  real(dl) :: x1, x2, y1, y2, mu
293 
294  ! initialize to null value:
295  equispacedlinearintepolatefunction1dsecondderivative = self%null_value
296  if ( self%has_null_value ) then
297  ! if outside the interpolation range return the null value:
298  if ( x <= self%x_initial .or. x >= self%x_final ) return
299  else
300  ! if below the interpolation range return the first value:
301  if ( x <= self%x_initial ) then
302  equispacedlinearintepolatefunction1dsecondderivative = self%ypp(1)
303  return
304  end if
305  ! if above the interpolation range return the first value:
306  if ( x >= self%x_final ) then
307  equispacedlinearintepolatefunction1dsecondderivative = self%ypp(self%num_points)
308  return
309  end if
310  end if
311 
312  ! return the index of the point:
313  if ( present(index) ) then
314  ind = index
315  else
316  ind = int( ( x-self%x_initial)/self%grid_width ) +1
317  end if
318 
319  ! get the interpolation coefficient:
320  if ( present(coeff) ) then
321  mu = coeff
322  else
323  ! store the x values:
324  x1 = self%x(ind)
325  x2 = self%x(ind+1)
326  ! compute the linear interpolation coefficient:
327  mu = (x-x1)/(x2-x1)
328  end if
329 
330  ! store the y values:
331  y1 = self%ypp(ind)
332  y2 = self%ypp(ind+1)
333 
334  ! compute the linear interpolation:
335  equispacedlinearintepolatefunction1dsecondderivative = y1*( 1._dl -mu ) +y2*mu
336 
337  end function equispacedlinearintepolatefunction1dsecondderivative
338 
339  ! ---------------------------------------------------------------------------------------------
341  function equispacedlinearintepolatefunction1dthirdderivative( self, x, index, coeff )
342 
343  implicit none
344 
345  class(equispaced_linear_interpolate_function_1d) :: self
346  real(dl), intent(in) :: x
347  integer , intent(in), optional :: index
348  real(dl), intent(in), optional :: coeff
349  real(dl) :: equispacedlinearintepolatefunction1dthirdderivative
350 
351  integer :: ind
352  real(dl) :: x1, x2, y1, y2, mu
353 
354  ! initialize to null value:
355  equispacedlinearintepolatefunction1dthirdderivative = self%null_value
356  if ( self%has_null_value ) then
357  ! if outside the interpolation range return the null value:
358  if ( x <= self%x_initial .or. x >= self%x_final ) return
359  else
360  ! if below the interpolation range return the first value:
361  if ( x <= self%x_initial ) then
362  equispacedlinearintepolatefunction1dthirdderivative = self%yppp(1)
363  return
364  end if
365  ! if above the interpolation range return the first value:
366  if ( x >= self%x_final ) then
367  equispacedlinearintepolatefunction1dthirdderivative = self%yppp(self%num_points)
368  return
369  end if
370  end if
371 
372  ! return the index of the point:
373  if ( present(index) ) then
374  ind = index
375  else
376  ind = int( ( x-self%x_initial)/self%grid_width ) +1
377  end if
378 
379  ! get the interpolation coefficient:
380  if ( present(coeff) ) then
381  mu = coeff
382  else
383  ! store the x values:
384  x1 = self%x(ind)
385  x2 = self%x(ind+1)
386  ! compute the linear interpolation coefficient:
387  mu = (x-x1)/(x2-x1)
388  end if
389 
390  ! store the y values:
391  y1 = self%yppp(ind)
392  y2 = self%yppp(ind+1)
393 
394  ! compute the linear interpolation:
395  equispacedlinearintepolatefunction1dthirdderivative = y1*( 1._dl -mu ) +y2*mu
396 
397  end function equispacedlinearintepolatefunction1dthirdderivative
398 
399  ! ---------------------------------------------------------------------------------------------
401  function equispacedlinearintepolatefunction1dintegral( self, x, index, coeff )
402 
403  implicit none
404 
405  class(equispaced_linear_interpolate_function_1d) :: self
406  real(dl), intent(in) :: x
407  integer , intent(in), optional :: index
408  real(dl), intent(in), optional :: coeff
409  real(dl) :: equispacedlinearintepolatefunction1dintegral
410 
411  integer :: ind
412  real(dl) :: x1, x2, y1, y2, mu
413 
414  ! initialize to null value:
415  equispacedlinearintepolatefunction1dintegral = self%null_value
416  if ( self%has_null_value ) then
417  ! if outside the interpolation range return the null value:
418  if ( x <= self%x_initial .or. x >= self%x_final ) return
419  else
420  ! if below the interpolation range return the first value:
421  if ( x <= self%x_initial ) then
422  equispacedlinearintepolatefunction1dintegral = self%yint(1)
423  return
424  end if
425  ! if above the interpolation range return the first value:
426  if ( x >= self%x_final ) then
427  equispacedlinearintepolatefunction1dintegral = self%yint(self%num_points)
428  return
429  end if
430  end if
431 
432  ! return the index of the point:
433  if ( present(index) ) then
434  ind = index
435  else
436  ind = int( ( x-self%x_initial)/self%grid_width ) +1
437  end if
438 
439  ! get the interpolation coefficient:
440  if ( present(coeff) ) then
441  mu = coeff
442  else
443  ! store the x values:
444  x1 = self%x(ind)
445  x2 = self%x(ind+1)
446  ! compute the linear interpolation coefficient:
447  mu = (x-x1)/(x2-x1)
448  end if
449 
450  ! store the y values:
451  y1 = self%yint(ind)
452  y2 = self%yint(ind+1)
453 
454  ! compute the linear interpolation:
455  equispacedlinearintepolatefunction1dintegral = y1*( 1._dl -mu ) +y2*mu
456 
457  end function equispacedlinearintepolatefunction1dintegral
458 
459  ! ---------------------------------------------------------------------------------------------
462  subroutine equispacedlinearintepolatefunction1dinitderivatives( self )
463 
464  implicit none
465 
466  class(equispaced_linear_interpolate_function_1d) :: self
467 
468  write(*,*) 'ERROR: not yet implemented (IW)'
469  stop
470 
471  end subroutine equispacedlinearintepolatefunction1dinitderivatives
472 
473  ! ---------------------------------------------------------------------------------------------
474 
476 
477 !----------------------------------------------------------------------------------------
This module contains the class that can be used for 1D linearly interpolated functions on an equispac...
This module contains various generic algorithms that are useful to EFTCAMB.