MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_lookup_table.t
Go to the documentation of this file.
1 !> A Fortran 90 module for creating 1D and 2D lookup tables. These tables can be
2 !> used to efficiently interpolate one or more values.
3 !>
4 !> Author: Jannis Teunissen
6  implicit none
7  private
8 
9  ! The precision of the real numbers used in the tables
10  integer, parameter :: dp = kind(1.0d0)
11 
12  ! ** Routines for finding indices in sorted lists **
13  public :: find_index_linear
14  public :: find_index_bsearch
15  public :: find_index_adaptive
16 
17  !> The lookup table type. There can be one or more columns, for which values
18  !> can be looked up for a given 'x-coordinate'.
19  type lt_t
20  integer :: n_points !< The number of points
21  integer :: n_cols !< The number of columns
22  real(dp) :: x_min !< The minimum lookup coordinate
23  real(dp) :: dx !< The x-spacing in the lookup coordinate
24  real(dp) :: inv_dx !< The inverse x-spacing
25 
26  ! The table is stored in two ways, to speed up different types of lookups.
27  real(dp), allocatable :: cols_rows(:, :) !< The table in column-major order
28  real(dp), allocatable :: rows_cols(:, :) !< The table in row-major order
29  end type lt_t
30 
31  !> The 2D lookup table type
32  type lt2_t
33  integer :: n_points(2) !< The size of the table
34  integer :: n_cols !< The number of columns/variables
35  real(dp) :: x_min(2) !< The minimum lookup coordinate
36  real(dp) :: dx(2) !< The x-spacing in the lookup coordinate
37  real(dp) :: inv_dx(2) !< The inverse x-spacing
38  real(dp), allocatable :: rows_cols(:, :, :)
39  end type lt2_t
40 
41  !> The 3D lookup table type
42  type lt3_t
43  integer :: n_points(3) !< The size of the table
44  integer :: n_cols !< The number of columns/variables
45  real(dp) :: x_min(3) !< The minimum lookup coordinate
46  real(dp) :: dx(3) !< The x-spacing in the lookup coordinate
47  real(dp) :: inv_dx(3) !< The inverse x-spacing
48  real(dp), allocatable :: rows_cols(:, :, :, :)
49  end type lt3_t
50 
51  !> Type to indicate a location in the lookup table, which can be used to speed
52  !> up multiple lookups of different columns.
53  type lt_loc_t
54  private
55  integer :: low_ix !< The x-value lies between low_ix and low_ix+1
56  real(dp) :: low_frac !< The distance from low_ix (up to low_ix+1), given
57  !< as a real number between 0 and 1.
58  end type lt_loc_t
59 
60  !> Type to indicate a location in a 2D lookup table
61  type lt2_loc_t
62  private
63  !> The x-value lies between low_ix and low_ix+1
64  integer :: low_ix(2)
65  !> The distance from low_ix (up to low_ix+1), given as a real number
66  !> between 0 and 1.
67  real(dp) :: low_frac(2)
68  end type lt2_loc_t
69 
70  !> Type to indicate a location in a 3D lookup table
71  type lt3_loc_t
72  private
73  !> The x-value lies between low_ix and low_ix+1
74  integer :: low_ix(3)
75  !> The distance from low_ix (up to low_ix+1), given as a real number
76  !> between 0 and 1.
77  real(dp) :: low_frac(3)
78  end type lt3_loc_t
79 
80  ! Public types
81  public :: lt_t
82  public :: lt_loc_t
83  public :: lt2_t
84  public :: lt2_loc_t
85  public :: lt3_t
86  public :: lt3_loc_t
87 
88  ! Public methods
89  public :: lt_create ! Create a new lookup table
90  public :: lt_create_from_data ! Create a new lookup table from existing data
91  public :: lt_get_xdata ! Get the x-values of a table
92  public :: lt_get_spaced_data ! Convert values to regularly spaced
93  public :: lt_set_col ! Set one table column
94  public :: lt_add_col ! Add a column
95  public :: lt_get_loc ! Get the index (row) of a value
96  public :: lt_get_col ! Interpolate one column
97  public :: lt_get_mcol ! Interpolate multiple columns
98  public :: lt_get_col_at_loc ! Get one column at location
99  public :: lt_get_mcol_at_loc ! Get multiple columns at location
100  public :: lt_get_data ! Get all the data of the table
101  public :: lt_lin_interp_list ! Linearly interpolate a list
102  public :: lt_to_file ! Store lookup table in file
103  public :: lt_from_file ! Restore lookup table from file
104 
105  ! Public methods
106  public :: lt2_create ! Create a new lookup table
107  public :: lt2_create_from_data ! Create a new lookup table from existing data
108  public :: lt2_set_col ! Set one table column
109  public :: lt2_get_loc ! Get the index (row) of a value
110  public :: lt2_get_col ! Interpolate one column
111  public :: lt2_get_col_at_loc ! Get one column at location
112 
113  public :: lt3_create ! Create a new lookup table
114  public :: lt3_create_from_data ! Create a new lookup table from existing data
115  public :: lt3_get_loc ! Get the index (row) of a value
116  public :: lt3_get_col ! Interpolate one column
117  public :: lt3_get_col_at_loc ! Get one column at location
118 
119 contains
120 
121  ! ** Routines for finding indices **
122 
123  !> Linear search of sorted list for the smallest ix such that list(ix) >= val.
124  !> On failure, returns size(list)+1
125  pure function find_index_linear(list, val) result(ix)
126  real(dp), intent(in) :: list(:) !< Sorted list
127  real(dp), intent(in) :: val !< Value to search for
128  integer :: ix
129 
130  do ix = 1, size(list)
131  if (list(ix) >= val) exit
132  end do
133  end function find_index_linear
134 
135  !> Binary search of sorted list for the smallest ix such that list(ix) >= val.
136  !> On failure, returns size(list)+1
137  pure function find_index_bsearch(list, val) result(ix)
138  real(dp), intent(in) :: list(:) !< Sorted list
139  real(dp), intent(in) :: val !< Value to search for
140  integer :: ix, i_min, i_max, i_middle
141 
142  i_min = 1
143  i_max = size(list)
144 
145  do while (i_min < i_max)
146  ! This safely performs: i_middle = (i_max + i_min) / 2
147  i_middle = i_min + ishft(i_max - i_min, -1)
148 
149  if (list(i_middle) >= val) then
150  i_max = i_middle
151  else
152  i_min = i_middle + 1
153  end if
154  end do
155 
156  ix = i_min
157  if (val > list(ix)) ix = size(list) + 1
158  end function find_index_bsearch
159 
160  !> Adaptive search (combination of linear and binary search) of sorted list
161  !> for the smallest ix such that list(ix) >= val. On failure, returns
162  !> size(list)+1
163  pure function find_index_adaptive(list, val) result(ix)
164  real(dp), intent(in) :: list(:) !< Sorted list
165  real(dp), intent(in) :: val !< Value to search for
166  integer :: ix
167  integer, parameter :: binary_search_limit = 40
168 
169  if (size(list) < binary_search_limit) then
170  ix = find_index_linear(list, val)
171  else
172  ix = find_index_bsearch(list, val)
173  end if
174  end function find_index_adaptive
175 
176  ! ** 1D lookup table routines **
177 
178  !> This function returns a new lookup table
179  function lt_create(x_min, x_max, n_points, n_cols) result(my_lt)
180  real(dp), intent(in) :: x_min !< Minimum x-coordinate
181  real(dp), intent(in) :: x_max !< Maximum x-coordinate
182  integer, intent(in) :: n_points !< How many x-values to store
183  integer, intent(in) :: n_cols !< Number of variables that will be looked up
184  type(lt_t) :: my_lt
185 
186  if (x_max <= x_min) stop "set_xdata: x_max should be > x_min"
187  if (n_points <= 1) stop "set_xdata: n_points should be bigger than 1"
188 
189  my_lt%n_points = n_points
190  my_lt%x_min = x_min
191  my_lt%dx = (x_max - x_min) / (n_points - 1)
192  my_lt%inv_dx = 1 / my_lt%dx
193 
194  allocate(my_lt%cols_rows(n_cols, n_points))
195  allocate(my_lt%rows_cols(n_points, n_cols))
196  my_lt%cols_rows = 0
197  my_lt%rows_cols = 0
198  my_lt%n_cols = n_cols
199  end function lt_create
200 
201  !> This function returns a new lookup table from regularly spaced data
202  function lt_create_from_data(x_min, x_max, spaced_data) result(my_lt)
203  real(dp), intent(in) :: x_min !< Minimum coordinate
204  real(dp), intent(in) :: x_max !< Maximum coordinate
205  real(dp), intent(in) :: spaced_data(:, :) !< Input data of shape (n_points, n_cols)
206  integer :: n_points, n_cols
207  type(lt_t) :: my_lt
208 
209  n_points = size(spaced_data, 1)
210  n_cols = size(spaced_data, 2)
211 
212  if (x_max <= x_min) stop "LT_create error: x_max <= x_min"
213  if (n_points <= 1) stop "LT_create error: n_points <= 1"
214 
215  my_lt%n_cols = n_cols
216  my_lt%n_points = n_points
217  my_lt%x_min = x_min
218  my_lt%dx = (x_max - x_min) / (n_points - 1)
219  my_lt%inv_dx = 1 / my_lt%dx
220 
221  my_lt%rows_cols = spaced_data
222  my_lt%cols_rows = transpose(spaced_data)
223  end function lt_create_from_data
224 
225  !> Returns the x-coordinates of the lookup table
226  pure function lt_get_xdata(x_min, dx, n_points) result(xdata)
227  real(dp), intent(in) :: x_min, dx
228  integer, intent(in) :: n_points
229  real(dp) :: xdata(n_points)
230  integer :: ix
231 
232  do ix = 1, n_points
233  xdata(ix) = x_min + (ix-1) * dx
234  end do
235  end function lt_get_xdata
236 
237  !> Linearly interpolate the (x, y) input data to the new_x coordinates
238  pure function lt_get_spaced_data(in_x, in_y, new_x) result(out_yy)
239  real(dp), intent(in) :: in_x(:), in_y(:), new_x(:)
240  real(dp) :: out_yy(size(new_x))
241  integer :: ix
242 
243  do ix = 1, size(new_x)
244  call lt_lin_interp_list(in_x, in_y, new_x(ix), out_yy(ix))
245  end do
246  end function lt_get_spaced_data
247 
248  !> Fill the column with index col_ix using the linearly interpolated (x, y)
249  !> data
250  pure subroutine lt_set_col(my_lt, col_ix, x, y)
251  type(lt_t), intent(inout) :: my_lt
252  integer, intent(in) :: col_ix
253  real(dp), intent(in) :: x(:), y(:)
254  my_lt%cols_rows(col_ix, :) = lt_get_spaced_data(x, y, &
255  lt_get_xdata(my_lt%x_min, my_lt%dx, my_lt%n_points))
256  my_lt%rows_cols(:, col_ix) = my_lt%cols_rows(col_ix, :)
257  end subroutine lt_set_col
258 
259  !> Add a new column by linearly interpolating the (x, y) data
260  pure subroutine lt_add_col(my_lt, x, y)
261  type(lt_t), intent(inout) :: my_lt
262  real(dp), intent(in) :: x(:), y(:)
263  type(lt_t) :: temp_lt
264 
265  temp_lt = my_lt
266  deallocate(my_lt%cols_rows)
267  deallocate(my_lt%rows_cols)
268  allocate(my_lt%cols_rows(my_lt%n_cols+1, my_lt%n_points))
269  allocate(my_lt%rows_cols(my_lt%n_points, my_lt%n_cols+1))
270 
271  my_lt%cols_rows(1:my_lt%n_cols, :) = temp_lt%cols_rows
272  my_lt%rows_cols(:, 1:my_lt%n_cols) = temp_lt%rows_cols
273  my_lt%n_cols = my_lt%n_cols + 1
274  my_lt%cols_rows(my_lt%n_cols, :) = lt_get_spaced_data(x, y, &
275  lt_get_xdata(my_lt%x_min, my_lt%dx, my_lt%n_points))
276  my_lt%rows_cols(:, my_lt%n_cols) = my_lt%cols_rows(my_lt%n_cols, :)
277  end subroutine lt_add_col
278 
279  !> Get a location in the lookup table
280  elemental function lt_get_loc(my_lt, x) result(my_loc)
281  type(lt_t), intent(in) :: my_lt
282  real(dp), intent(in) :: x
283  type(lt_loc_t) :: my_loc
284  real(dp) :: frac
285 
286  frac = (x - my_lt%x_min) * my_lt%inv_dx
287  my_loc%low_ix = ceiling(frac)
288  my_loc%low_frac = my_loc%low_ix - frac
289 
290  ! Check bounds
291  if (my_loc%low_ix < 1) then
292  my_loc%low_ix = 1
293  my_loc%low_frac = 1
294  else if (my_loc%low_ix >= my_lt%n_points) then
295  my_loc%low_ix = my_lt%n_points - 1
296  my_loc%low_frac = 0
297  end if
298  end function lt_get_loc
299 
300  !> Get the values of all columns at x
301  pure function lt_get_mcol(my_lt, x) result(col_values)
302  type(lt_t), intent(in) :: my_lt
303  real(dp), intent(in) :: x
304  real(dp) :: col_values(my_lt%n_cols)
305  type(lt_loc_t) :: loc
306 
307  loc = lt_get_loc(my_lt, x)
308  col_values = lt_get_mcol_at_loc(my_lt, loc)
309  end function lt_get_mcol
310 
311  !> Get the value of a single column at x
312  elemental function lt_get_col(my_lt, col_ix, x) result(col_value)
313  type(lt_t), intent(in) :: my_lt
314  integer, intent(in) :: col_ix
315  real(dp), intent(in) :: x
316  real(dp) :: col_value
317  type(lt_loc_t) :: loc
318 
319  loc = lt_get_loc(my_lt, x)
320  col_value = lt_get_col_at_loc(my_lt, col_ix, loc)
321  end function lt_get_col
322 
323  !> Get the values of all columns at a location
324  pure function lt_get_mcol_at_loc(my_lt, loc) result(col_values)
325  type(lt_t), intent(in) :: my_lt
326  type(lt_loc_t), intent(in) :: loc
327  real(dp) :: col_values(my_lt%n_cols)
328 
329  col_values = loc%low_frac * my_lt%cols_rows(:, loc%low_ix) + &
330  (1-loc%low_frac) * my_lt%cols_rows(:, loc%low_ix+1)
331  end function lt_get_mcol_at_loc
332 
333  !> Get the value of a single column at a location
334  elemental function lt_get_col_at_loc(my_lt, col_ix, loc) result(col_value)
335  type(lt_t), intent(in) :: my_lt
336  integer, intent(in) :: col_ix
337  type(lt_loc_t), intent(in) :: loc
338  real(dp) :: col_value
339 
340  col_value = loc%low_frac * my_lt%rows_cols(loc%low_ix, col_ix) + &
341  (1-loc%low_frac) * my_lt%rows_cols(loc%low_ix+1, col_ix)
342  end function lt_get_col_at_loc
343 
344  !> Get the x-coordinates and the columns of the lookup table
345  pure subroutine lt_get_data(my_lt, x_data, cols_rows)
346  type(lt_t), intent(in) :: my_lt
347  real(dp), intent(out) :: x_data(:), cols_rows(:, :)
348 
349  x_data = lt_get_xdata(my_lt%x_min, my_lt%dx, my_lt%n_points)
350  cols_rows = my_lt%cols_rows
351  end subroutine lt_get_data
352 
353  !> Compute by use of linear interpolation the value in the middle
354  ! of a domain D = [x_list(1) , x_list(size(x_list))].
355  ! If x_value is left of domain D,
356  ! then the value becomes the value at the left side of D,
357  ! if x_value is right of domain D,
358  ! then the value becomes the value at the rigth side of D
359  pure subroutine lt_lin_interp_list(x_list, y_list, x_value, y_value)
360  real(dp), intent(in) :: x_list(:), y_list(:)
361  real(dp), intent(in) :: x_value
362  real(dp), intent(out) :: y_value
363 
364  integer :: ix, iMin, iMax
365  real(dp) :: temp
366 
367  imin = 1
368  imax = size(x_list)
369 
370  if (x_value <= x_list(imin)) then
371  y_value = y_list(imin)
372  else if (x_value >= x_list(imax)) then
373  y_value = y_list(imax)
374  else
375  ix = find_index_adaptive(x_list, x_value)
376  temp = (x_value - x_list(ix-1)) / (x_list(ix) - x_list(ix-1))
377  y_value = (1 - temp) * y_list(ix-1) + temp * y_list(ix)
378  end if
379  end subroutine lt_lin_interp_list
380 
381  !> Write the lookup table to file (in binary, potentially unportable)
382  subroutine lt_to_file(my_lt, filename)
383  type(lt_t), intent(in) :: my_lt
384  character(len=*), intent(in) :: filename
385  integer :: my_unit
386 
387  open(newunit=my_unit, file=trim(filename), form='UNFORMATTED', &
388  access='STREAM', status='REPLACE')
389  write(my_unit) my_lt%n_points, my_lt%n_cols
390  write(my_unit) my_lt%x_min, my_lt%dx, my_lt%inv_dx
391  write(my_unit) my_lt%cols_rows
392  close(my_unit)
393  end subroutine lt_to_file
394 
395  !> Read the lookup table from file (in binary, potentially unportable)
396  subroutine lt_from_file(my_lt, filename)
397  type(lt_t), intent(inout) :: my_lt
398  character(len=*), intent(in) :: filename
399  integer :: my_unit
400 
401  open(newunit=my_unit, file=trim(filename), form='UNFORMATTED', &
402  access='STREAM', status='OLD')
403  read(my_unit) my_lt%n_points, my_lt%n_cols
404  read(my_unit) my_lt%x_min, my_lt%dx, my_lt%inv_dx
405 
406  allocate(my_lt%cols_rows(my_lt%n_cols, my_lt%n_points))
407  allocate(my_lt%rows_cols(my_lt%n_points, my_lt%n_cols))
408 
409  read(my_unit) my_lt%cols_rows
410  my_lt%rows_cols = transpose(my_lt%cols_rows)
411 
412  close(my_unit)
413  end subroutine lt_from_file
414 
415  ! ** 2D lookup table routines **
416 
417  !> This function returns a new lookup table
418  function lt2_create(x_min, x_max, n_points, n_cols) result(my_lt)
419  real(dp), intent(in) :: x_min(2) !< Minimum coordinate
420  real(dp), intent(in) :: x_max(2) !< Maximum coordinate
421  integer, intent(in) :: n_points(2) !< How many values to store
422  integer, intent(in) :: n_cols !< Number of variables that will be looked up
423  type(lt2_t) :: my_lt
424 
425  if (any(x_max <= x_min)) stop "LT2_create error: x_max <= x_min"
426  if (any(n_points <= 1)) stop "LT2_create error: n_points <= 1"
427 
428  my_lt%n_points = n_points
429  my_lt%x_min = x_min
430  my_lt%dx = (x_max - x_min) / (n_points - 1)
431  my_lt%inv_dx = 1 / my_lt%dx
432 
433  allocate(my_lt%rows_cols(n_points(1), n_points(2), n_cols))
434  my_lt%rows_cols = 0
435  my_lt%n_cols = n_cols
436  end function lt2_create
437 
438  !> This function returns a new lookup table from regularly spaced data
439  function lt2_create_from_data(x_min, x_max, spaced_data) result(my_lt)
440  real(dp), intent(in) :: x_min(2) !< Minimum coordinate
441  real(dp), intent(in) :: x_max(2) !< Maximum coordinate
442  real(dp), intent(in) :: spaced_data(:, :, :) !< Input data of shape (nx, ny, n_cols)
443  integer :: n_points(2), n_cols
444  type(lt2_t) :: my_lt
445 
446  n_points(1) = size(spaced_data, 1)
447  n_points(2) = size(spaced_data, 2)
448  n_cols = size(spaced_data, 3)
449 
450  if (any(x_max <= x_min)) stop "LT2_create error: x_max <= x_min"
451  if (any(n_points <= 1)) stop "LT2_create error: n_points <= 1"
452 
453  my_lt%n_cols = n_cols
454  my_lt%n_points = n_points
455  my_lt%x_min = x_min
456  my_lt%dx = (x_max - x_min) / (n_points - 1)
457  my_lt%inv_dx = 1 / my_lt%dx
458 
459  my_lt%rows_cols = spaced_data
460  end function lt2_create_from_data
461 
462  !> Fill the column with index col_ix using linearly interpolated data
463  pure subroutine lt2_set_col(my_lt, col_ix, x1, x2, y)
464  type(lt2_t), intent(inout) :: my_lt
465  integer, intent(in) :: col_ix
466  real(dp), intent(in) :: x1(:), x2(:), y(:, :)
467  real(dp), allocatable :: tmp(:, :), c1(:), c2(:)
468  integer :: ix
469 
470  c1 = lt_get_xdata(my_lt%x_min(1), my_lt%dx(1), my_lt%n_points(1))
471  c2 = lt_get_xdata(my_lt%x_min(2), my_lt%dx(2), my_lt%n_points(2))
472  allocate(tmp(my_lt%n_points(1), size(x2)))
473 
474  ! Interpolate along first coordinate
475  do ix = 1, size(x2)
476  tmp(:, ix) = lt_get_spaced_data(x1, y(:, ix), c1)
477  end do
478 
479  ! Interpolate along second coordinate
480  do ix = 1, my_lt%n_points(1)
481  my_lt%rows_cols(ix, :, col_ix) = &
482  lt_get_spaced_data(x2, tmp(ix, :), c2)
483  end do
484  end subroutine lt2_set_col
485 
486  !> Get a location in the lookup table
487  elemental function lt2_get_loc(my_lt, x1, x2) result(my_loc)
488  type(lt2_t), intent(in) :: my_lt
489  real(dp), intent(in) :: x1, x2
490  type(lt2_loc_t) :: my_loc
491  real(dp) :: frac(2)
492 
493  frac = ([x1, x2] - my_lt%x_min) * my_lt%inv_dx
494  my_loc%low_ix = ceiling(frac)
495  my_loc%low_frac = my_loc%low_ix - frac
496 
497  ! Check bounds
498  where (my_loc%low_ix < 1)
499  my_loc%low_ix = 1
500  my_loc%low_frac = 1
501  end where
502 
503  where (my_loc%low_ix >= my_lt%n_points)
504  my_loc%low_ix = my_lt%n_points - 1
505  my_loc%low_frac = 0
506  end where
507  end function lt2_get_loc
508 
509  !> Get the value of a single column at x
510  elemental function lt2_get_col(my_lt, col_ix, x1, x2) result(col_value)
511  type(lt2_t), intent(in) :: my_lt
512  integer, intent(in) :: col_ix
513  real(dp), intent(in) :: x1, x2
514  real(dp) :: col_value
515  type(lt2_loc_t) :: loc
516 
517  loc = lt2_get_loc(my_lt, x1, x2)
518  col_value = lt2_get_col_at_loc(my_lt, col_ix, loc)
519  end function lt2_get_col
520 
521  !> Get the value of a single column at a location
522  elemental function lt2_get_col_at_loc(my_lt, col_ix, loc) result(col_value)
523  type(lt2_t), intent(in) :: my_lt
524  integer, intent(in) :: col_ix
525  type(lt2_loc_t), intent(in) :: loc
526  integer :: ix(2)
527  real(dp) :: w(2, 2)
528  real(dp) :: col_value
529 
530  ! Bilinear interpolation
531  w(1, 1) = loc%low_frac(1) * loc%low_frac(2)
532  w(2, 1) = (1 - loc%low_frac(1)) * loc%low_frac(2)
533  w(1, 2) = loc%low_frac(1) * (1 - loc%low_frac(2))
534  w(2, 2) = (1 - loc%low_frac(1)) * (1 - loc%low_frac(2))
535  ix = loc%low_ix
536 
537  col_value = sum(w * my_lt%rows_cols(ix(1):ix(1)+1, &
538  ix(2):ix(2)+1, col_ix))
539  end function lt2_get_col_at_loc
540 
541  ! ** 3D lookup table routines **
542 
543  !> This function returns a new lookup table
544  function lt3_create(x_min, x_max, n_points, n_cols) result(my_lt)
545  real(dp), intent(in) :: x_min(3) !< Minimum coordinate
546  real(dp), intent(in) :: x_max(3) !< Maximum coordinate
547  integer, intent(in) :: n_points(3) !< How many values to store
548  integer, intent(in) :: n_cols !< Number of variables that will be looked up
549  type(lt3_t) :: my_lt
550 
551  if (any(x_max <= x_min)) stop "LT3_create error: x_max <= x_min"
552  if (any(n_points <= 1)) stop "LT3_create error: n_points <= 1"
553 
554  my_lt%n_points = n_points
555  my_lt%x_min = x_min
556  my_lt%dx = (x_max - x_min) / (n_points - 1)
557  my_lt%inv_dx = 1 / my_lt%dx
558 
559  allocate(my_lt%rows_cols(n_points(1), n_points(2), n_points(3), n_cols))
560  my_lt%rows_cols = 0
561  my_lt%n_cols = n_cols
562  end function lt3_create
563 
564  !> This function returns a new lookup table from regularly spaced data
565  function lt3_create_from_data(x_min, x_max, spaced_data) result(my_lt)
566  real(dp), intent(in) :: x_min(3) !< Minimum coordinate
567  real(dp), intent(in) :: x_max(3) !< Maximum coordinate
568  real(dp), intent(in) :: spaced_data(:, :, :, :) !< Input data of shape (nx, ny, nz, n_cols)
569  integer :: n_points(3), n_cols
570  type(lt3_t) :: my_lt
571 
572  n_points(1) = size(spaced_data, 1)
573  n_points(2) = size(spaced_data, 2)
574  n_points(3) = size(spaced_data, 3)
575  n_cols = size(spaced_data, 4)
576 
577  if (any(x_max <= x_min)) stop "LT3_create error: x_max <= x_min"
578  if (any(n_points <= 1)) stop "LT3_create error: n_points <= 1"
579 
580  my_lt%n_cols = n_cols
581  my_lt%n_points = n_points
582  my_lt%x_min = x_min
583  my_lt%dx = (x_max - x_min) / (n_points - 1)
584  my_lt%inv_dx = 1 / my_lt%dx
585 
586  my_lt%rows_cols = spaced_data
587  end function lt3_create_from_data
588 
589  !> Get a location in the lookup table
590  elemental function lt3_get_loc(my_lt, x1, x2, x3) result(my_loc)
591  type(lt3_t), intent(in) :: my_lt
592  real(dp), intent(in) :: x1, x2, x3
593  type(lt3_loc_t) :: my_loc
594  real(dp) :: frac(3)
595 
596  frac = ([x1, x2, x3] - my_lt%x_min) * my_lt%inv_dx
597  my_loc%low_ix = ceiling(frac)
598  my_loc%low_frac = my_loc%low_ix - frac
599 
600  ! Check bounds
601  where (my_loc%low_ix < 1)
602  my_loc%low_ix = 1
603  my_loc%low_frac = 1
604  end where
605 
606  where (my_loc%low_ix >= my_lt%n_points)
607  my_loc%low_ix = my_lt%n_points - 1
608  my_loc%low_frac = 0
609  end where
610  end function lt3_get_loc
611 
612  !> Get the value of a single column at x
613  elemental function lt3_get_col(my_lt, col_ix, x1, x2, x3) result(col_value)
614  type(lt3_t), intent(in) :: my_lt
615  integer, intent(in) :: col_ix
616  real(dp), intent(in) :: x1, x2, x3
617  real(dp) :: col_value
618  type(lt3_loc_t) :: loc
619 
620  loc = lt3_get_loc(my_lt, x1, x2, x3)
621  col_value = lt3_get_col_at_loc(my_lt, col_ix, loc)
622  end function lt3_get_col
623 
624  !> Get the value of a single column at a location
625  elemental function lt3_get_col_at_loc(my_lt, col_ix, loc) result(col_value)
626  type(lt3_t), intent(in) :: my_lt
627  integer, intent(in) :: col_ix
628  type(lt3_loc_t), intent(in) :: loc
629  integer :: ix(3)
630  real(dp) :: w(2, 2, 2)
631  real(dp) :: col_value
632 
633  ! Bilinear interpolation
634  w(1, 1, 1) = loc%low_frac(1) * loc%low_frac(2)
635  w(2, 1, 1) = (1 - loc%low_frac(1)) * loc%low_frac(2)
636  w(1, 2, 1) = loc%low_frac(1) * (1 - loc%low_frac(2))
637  w(2, 2, 1) = (1 - loc%low_frac(1)) * (1 - loc%low_frac(2))
638  w(:, :, 2) = w(:, :, 1) * (1 - loc%low_frac(3))
639  w(:, :, 1) = w(:, :, 1) * loc%low_frac(3)
640  ix = loc%low_ix
641 
642  col_value = sum(w * my_lt%rows_cols(ix(1):ix(1)+1, &
643  ix(2):ix(2)+1, ix(3):ix(3)+1, col_ix))
644  end function lt3_get_col_at_loc
645 
646 end module mod_lookup_table
pure real(dp) function, dimension(n_points), public lt_get_xdata(x_min, dx, n_points)
Returns the x-coordinates of the lookup table.
elemental real(dp) function, public lt_get_col(my_lt, col_ix, x)
Get the value of a single column at x.
elemental real(dp) function, public lt_get_col_at_loc(my_lt, col_ix, loc)
Get the value of a single column at a location.
pure subroutine, public lt_lin_interp_list(x_list, y_list, x_value, y_value)
Compute by use of linear interpolation the value in the middle.
type(lt3_t) function, public lt3_create(x_min, x_max, n_points, n_cols)
This function returns a new lookup table.
type(lt3_t) function, public lt3_create_from_data(x_min, x_max, spaced_data)
This function returns a new lookup table from regularly spaced data.
elemental real(dp) function, public lt3_get_col_at_loc(my_lt, col_ix, loc)
Get the value of a single column at a location.
elemental type(lt2_loc_t) function, public lt2_get_loc(my_lt, x1, x2)
Get a location in the lookup table.
pure integer function, public find_index_bsearch(list, val)
Binary search of sorted list for the smallest ix such that list(ix) >= val. On failure, returns size(list)+1.
elemental type(lt_loc_t) function, public lt_get_loc(my_lt, x)
Get a location in the lookup table.
pure real(dp) function, dimension(my_lt%n_cols), public lt_get_mcol_at_loc(my_lt, loc)
Get the values of all columns at a location.
A Fortran 90 module for creating 1D and 2D lookup tables. These tables can be used to efficiently int...
subroutine, public lt_to_file(my_lt, filename)
Write the lookup table to file (in binary, potentially unportable)
pure subroutine, public lt2_set_col(my_lt, col_ix, x1, x2, y)
Fill the column with index col_ix using linearly interpolated data.
pure subroutine, public lt_add_col(my_lt, x, y)
Add a new column by linearly interpolating the (x, y) data.
type(lt_t) function, public lt_create_from_data(x_min, x_max, spaced_data)
This function returns a new lookup table from regularly spaced data.
pure subroutine, public lt_get_data(my_lt, x_data, cols_rows)
Get the x-coordinates and the columns of the lookup table.
type(lt2_t) function, public lt2_create(x_min, x_max, n_points, n_cols)
This function returns a new lookup table.
pure integer function, public find_index_adaptive(list, val)
Adaptive search (combination of linear and binary search) of sorted list for the smallest ix such tha...
elemental real(dp) function, public lt2_get_col(my_lt, col_ix, x1, x2)
Get the value of a single column at x.
elemental real(dp) function, public lt2_get_col_at_loc(my_lt, col_ix, loc)
Get the value of a single column at a location.
elemental real(dp) function, public lt3_get_col(my_lt, col_ix, x1, x2, x3)
Get the value of a single column at x.
pure real(dp) function, dimension(my_lt%n_cols), public lt_get_mcol(my_lt, x)
Get the values of all columns at x.
elemental type(lt3_loc_t) function, public lt3_get_loc(my_lt, x1, x2, x3)
Get a location in the lookup table.
subroutine, public lt_from_file(my_lt, filename)
Read the lookup table from file (in binary, potentially unportable)
pure integer function, public find_index_linear(list, val)
Linear search of sorted list for the smallest ix such that list(ix) >= val. On failure, returns size(list)+1.
type(lt_t) function, public lt_create(x_min, x_max, n_points, n_cols)
This function returns a new lookup table.
pure subroutine, public lt_set_col(my_lt, col_ix, x, y)
Fill the column with index col_ix using the linearly interpolated (x, y) data.
pure real(dp) function, dimension(size(new_x)), public lt_get_spaced_data(in_x, in_y, new_x)
Linearly interpolate the (x, y) input data to the new_x coordinates.
type(lt2_t) function, public lt2_create_from_data(x_min, x_max, spaced_data)
This function returns a new lookup table from regularly spaced data.