MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
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.
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
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
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
119contains
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 allocate(c1(my_lt%n_points(1)),c2(my_lt%n_points(2)))
471 c1 = lt_get_xdata(my_lt%x_min(1), my_lt%dx(1), my_lt%n_points(1))
472 c2 = lt_get_xdata(my_lt%x_min(2), my_lt%dx(2), my_lt%n_points(2))
473 allocate(tmp(my_lt%n_points(1), size(x2)))
474
475 ! Interpolate along first coordinate
476 do ix = 1, size(x2)
477 tmp(:, ix) = lt_get_spaced_data(x1, y(:, ix), c1)
478 end do
479
480 ! Interpolate along second coordinate
481 do ix = 1, my_lt%n_points(1)
482 my_lt%rows_cols(ix, :, col_ix) = &
483 lt_get_spaced_data(x2, tmp(ix, :), c2)
484 end do
485 deallocate(c1,c2)
486 end subroutine lt2_set_col
487
488 !> Get a location in the lookup table
489 elemental function lt2_get_loc(my_lt, x1, x2) result(my_loc)
490 type(lt2_t), intent(in) :: my_lt
491 real(dp), intent(in) :: x1, x2
492 type(lt2_loc_t) :: my_loc
493 real(dp) :: frac(2)
494
495 frac = ([x1, x2] - my_lt%x_min) * my_lt%inv_dx
496 my_loc%low_ix = ceiling(frac)
497 my_loc%low_frac = my_loc%low_ix - frac
498
499 ! Check bounds
500 where (my_loc%low_ix < 1)
501 my_loc%low_ix = 1
502 my_loc%low_frac = 1
503 end where
504
505 where (my_loc%low_ix >= my_lt%n_points)
506 my_loc%low_ix = my_lt%n_points - 1
507 my_loc%low_frac = 0
508 end where
509 end function lt2_get_loc
510
511 !> Get the value of a single column at x
512 elemental function lt2_get_col(my_lt, col_ix, x1, x2) result(col_value)
513 type(lt2_t), intent(in) :: my_lt
514 integer, intent(in) :: col_ix
515 real(dp), intent(in) :: x1, x2
516 real(dp) :: col_value
517 type(lt2_loc_t) :: loc
518
519 loc = lt2_get_loc(my_lt, x1, x2)
520 col_value = lt2_get_col_at_loc(my_lt, col_ix, loc)
521 end function lt2_get_col
522
523 !> Get the value of a single column at a location
524 elemental function lt2_get_col_at_loc(my_lt, col_ix, loc) result(col_value)
525 type(lt2_t), intent(in) :: my_lt
526 integer, intent(in) :: col_ix
527 type(lt2_loc_t), intent(in) :: loc
528 integer :: ix(2)
529 real(dp) :: w(2, 2)
530 real(dp) :: col_value
531
532 ! Bilinear interpolation
533 w(1, 1) = loc%low_frac(1) * loc%low_frac(2)
534 w(2, 1) = (1 - loc%low_frac(1)) * loc%low_frac(2)
535 w(1, 2) = loc%low_frac(1) * (1 - loc%low_frac(2))
536 w(2, 2) = (1 - loc%low_frac(1)) * (1 - loc%low_frac(2))
537 ix = loc%low_ix
538
539 col_value = sum(w * my_lt%rows_cols(ix(1):ix(1)+1, &
540 ix(2):ix(2)+1, col_ix))
541 end function lt2_get_col_at_loc
542
543 ! ** 3D lookup table routines **
544
545 !> This function returns a new lookup table
546 function lt3_create(x_min, x_max, n_points, n_cols) result(my_lt)
547 real(dp), intent(in) :: x_min(3) !< Minimum coordinate
548 real(dp), intent(in) :: x_max(3) !< Maximum coordinate
549 integer, intent(in) :: n_points(3) !< How many values to store
550 integer, intent(in) :: n_cols !< Number of variables that will be looked up
551 type(lt3_t) :: my_lt
552
553 if (any(x_max <= x_min)) stop "LT3_create error: x_max <= x_min"
554 if (any(n_points <= 1)) stop "LT3_create error: n_points <= 1"
555
556 my_lt%n_points = n_points
557 my_lt%x_min = x_min
558 my_lt%dx = (x_max - x_min) / (n_points - 1)
559 my_lt%inv_dx = 1 / my_lt%dx
560
561 allocate(my_lt%rows_cols(n_points(1), n_points(2), n_points(3), n_cols))
562 my_lt%rows_cols = 0
563 my_lt%n_cols = n_cols
564 end function lt3_create
565
566 !> This function returns a new lookup table from regularly spaced data
567 function lt3_create_from_data(x_min, x_max, spaced_data) result(my_lt)
568 real(dp), intent(in) :: x_min(3) !< Minimum coordinate
569 real(dp), intent(in) :: x_max(3) !< Maximum coordinate
570 real(dp), intent(in) :: spaced_data(:, :, :, :) !< Input data of shape (nx, ny, nz, n_cols)
571 integer :: n_points(3), n_cols
572 type(lt3_t) :: my_lt
573
574 n_points(1) = size(spaced_data, 1)
575 n_points(2) = size(spaced_data, 2)
576 n_points(3) = size(spaced_data, 3)
577 n_cols = size(spaced_data, 4)
578
579 if (any(x_max <= x_min)) stop "LT3_create error: x_max <= x_min"
580 if (any(n_points <= 1)) stop "LT3_create error: n_points <= 1"
581
582 my_lt%n_cols = n_cols
583 my_lt%n_points = n_points
584 my_lt%x_min = x_min
585 my_lt%dx = (x_max - x_min) / (n_points - 1)
586 my_lt%inv_dx = 1 / my_lt%dx
587 !TODO modified so that it does not rely on automatic allocation
588 allocate(my_lt%rows_cols(n_points(1), n_points(2), n_points(3), n_cols))
589 my_lt%rows_cols = spaced_data
590 end function lt3_create_from_data
591
592 !> Get a location in the lookup table
593 elemental function lt3_get_loc(my_lt, x1, x2, x3) result(my_loc)
594 type(lt3_t), intent(in) :: my_lt
595 real(dp), intent(in) :: x1, x2, x3
596 type(lt3_loc_t) :: my_loc
597 real(dp) :: frac(3)
598
599 frac = ([x1, x2, x3] - my_lt%x_min) * my_lt%inv_dx
600 my_loc%low_ix = ceiling(frac)
601 my_loc%low_frac = my_loc%low_ix - frac
602
603 ! Check bounds
604 where (my_loc%low_ix < 1)
605 my_loc%low_ix = 1
606 my_loc%low_frac = 1
607 end where
608
609 where (my_loc%low_ix >= my_lt%n_points)
610 my_loc%low_ix = my_lt%n_points - 1
611 my_loc%low_frac = 0
612 end where
613 end function lt3_get_loc
614
615 !> Get the value of a single column at x
616 elemental function lt3_get_col(my_lt, col_ix, x1, x2, x3) result(col_value)
617 type(lt3_t), intent(in) :: my_lt
618 integer, intent(in) :: col_ix
619 real(dp), intent(in) :: x1, x2, x3
620 real(dp) :: col_value
621 type(lt3_loc_t) :: loc
622
623 loc = lt3_get_loc(my_lt, x1, x2, x3)
624 col_value = lt3_get_col_at_loc(my_lt, col_ix, loc)
625 end function lt3_get_col
626
627 !> Get the value of a single column at a location
628 elemental function lt3_get_col_at_loc(my_lt, col_ix, loc) result(col_value)
629 type(lt3_t), intent(in) :: my_lt
630 integer, intent(in) :: col_ix
631 type(lt3_loc_t), intent(in) :: loc
632 integer :: ix(3)
633 real(dp) :: w(2, 2, 2)
634 real(dp) :: col_value
635
636 ! Bilinear interpolation
637 w(1, 1, 1) = loc%low_frac(1) * loc%low_frac(2)
638 w(2, 1, 1) = (1 - loc%low_frac(1)) * loc%low_frac(2)
639 w(1, 2, 1) = loc%low_frac(1) * (1 - loc%low_frac(2))
640 w(2, 2, 1) = (1 - loc%low_frac(1)) * (1 - loc%low_frac(2))
641 w(:, :, 2) = w(:, :, 1) * (1 - loc%low_frac(3))
642 w(:, :, 1) = w(:, :, 1) * loc%low_frac(3)
643 ix = loc%low_ix
644
645 col_value = sum(w * my_lt%rows_cols(ix(1):ix(1)+1, &
646 ix(2):ix(2)+1, ix(3):ix(3)+1, col_ix))
647 end function lt3_get_col_at_loc
648
649end module mod_lookup_table
A Fortran 90 module for creating 1D and 2D lookup tables. These tables can be used to efficiently int...
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.
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 lt_get_col_at_loc(my_lt, col_ix, loc)
Get the value of a single column at a location.
elemental real(dp) function, public lt_get_col(my_lt, col_ix, x)
Get the value of a single column at x.
pure subroutine, public lt_get_data(my_lt, x_data, cols_rows)
Get the x-coordinates and the columns of the lookup table.
pure subroutine, public lt2_set_col(my_lt, col_ix, x1, x2, y)
Fill the column with index col_ix using linearly interpolated data.
type(lt_t) function, public lt_create(x_min, x_max, n_points, n_cols)
This function returns a new lookup table.
elemental real(dp) function, public lt2_get_col(my_lt, col_ix, x1, x2)
Get the value of a single column at x.
type(lt3_t) function, public lt3_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.
elemental real(dp) function, public lt3_get_col(my_lt, col_ix, x1, x2, x3)
Get the value of a single column at x.
elemental type(lt3_loc_t) function, public lt3_get_loc(my_lt, x1, x2, x3)
Get a location in the lookup table.
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.
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.
elemental type(lt_loc_t) function, public lt_get_loc(my_lt, x)
Get a location in the lookup table.
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(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.
subroutine, public lt_to_file(my_lt, filename)
Write the lookup table to file (in binary, potentially unportable)
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.
subroutine, public lt_from_file(my_lt, filename)
Read the lookup table from file (in binary, potentially unportable)
pure real(dp) function, dimension(n_points), public lt_get_xdata(x_min, dx, n_points)
Returns the x-coordinates of 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.
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,...
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_linear(list, val)
Linear search of sorted list for the smallest ix such that list(ix) >= val. On failure,...
elemental type(lt2_loc_t) function, public lt2_get_loc(my_lt, x1, x2)
Get a location in the lookup table.
pure real(dp) function, dimension(my_lt%n_cols), public lt_get_mcol(my_lt, x)
Get the values of all columns at x.
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...
pure subroutine, public lt_add_col(my_lt, x, y)
Add a new column by linearly interpolating the (x, y) data.
Type to indicate a location in a 2D lookup table.
The 2D lookup table type.
Type to indicate a location in a 3D lookup table.
The 3D lookup table type.
Type to indicate a location in the lookup table, which can be used to speed up multiple lookups of di...
The lookup table type. There can be one or more columns, for which values can be looked up for a give...