10 integer,
parameter :: dp = kind(1.0d0)
27 real(dp),
allocatable :: cols_rows(:, :)
28 real(dp),
allocatable :: rows_cols(:, :)
33 integer :: n_points(2)
38 real(dp),
allocatable :: rows_cols(:, :, :)
43 integer :: n_points(3)
48 real(dp),
allocatable :: rows_cols(:, :, :, :)
67 real(dp) :: low_frac(2)
77 real(dp) :: low_frac(3)
126 real(dp),
intent(in) :: list(:)
127 real(dp),
intent(in) :: val
130 do ix = 1,
size(list)
131 if (list(ix) >= val)
exit
138 real(dp),
intent(in) :: list(:)
139 real(dp),
intent(in) :: val
140 integer :: ix, i_min, i_max, i_middle
145 do while (i_min < i_max)
147 i_middle = i_min + ishft(i_max - i_min, -1)
149 if (list(i_middle) >= val)
then
157 if (val > list(ix)) ix =
size(list) + 1
164 real(dp),
intent(in) :: list(:)
165 real(dp),
intent(in) :: val
167 integer,
parameter :: binary_search_limit = 40
169 if (
size(list) < binary_search_limit)
then
179 function lt_create(x_min, x_max, n_points, n_cols)
result(my_lt)
180 real(dp),
intent(in) :: x_min
181 real(dp),
intent(in) :: x_max
182 integer,
intent(in) :: n_points
183 integer,
intent(in) :: n_cols
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"
189 my_lt%n_points = n_points
191 my_lt%dx = (x_max - x_min) / (n_points - 1)
192 my_lt%inv_dx = 1 / my_lt%dx
194 allocate(my_lt%cols_rows(n_cols, n_points))
195 allocate(my_lt%rows_cols(n_points, n_cols))
198 my_lt%n_cols = n_cols
203 real(dp),
intent(in) :: x_min
204 real(dp),
intent(in) :: x_max
205 real(dp),
intent(in) :: spaced_data(:, :)
206 integer :: n_points, n_cols
209 n_points =
size(spaced_data, 1)
210 n_cols =
size(spaced_data, 2)
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"
215 my_lt%n_cols = n_cols
216 my_lt%n_points = n_points
218 my_lt%dx = (x_max - x_min) / (n_points - 1)
219 my_lt%inv_dx = 1 / my_lt%dx
221 my_lt%rows_cols = spaced_data
222 my_lt%cols_rows = transpose(spaced_data)
227 real(dp),
intent(in) :: x_min, dx
228 integer,
intent(in) :: n_points
229 real(dp) :: xdata(n_points)
233 xdata(ix) = x_min + (ix-1) * dx
239 real(dp),
intent(in) :: in_x(:), in_y(:), new_x(:)
240 real(dp) :: out_yy(size(new_x))
243 do ix = 1,
size(new_x)
251 type(
lt_t),
intent(inout) :: my_lt
252 integer,
intent(in) :: col_ix
253 real(dp),
intent(in) :: x(:), y(:)
256 my_lt%rows_cols(:, col_ix) = my_lt%cols_rows(col_ix, :)
261 type(
lt_t),
intent(inout) :: my_lt
262 real(dp),
intent(in) :: x(:), y(:)
263 type(
lt_t) :: temp_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))
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
276 my_lt%rows_cols(:, my_lt%n_cols) = my_lt%cols_rows(my_lt%n_cols, :)
281 type(
lt_t),
intent(in) :: my_lt
282 real(dp),
intent(in) :: x
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
291 if (my_loc%low_ix < 1)
then
294 else if (my_loc%low_ix >= my_lt%n_points)
then
295 my_loc%low_ix = my_lt%n_points - 1
302 type(
lt_t),
intent(in) :: my_lt
303 real(dp),
intent(in) :: x
304 real(dp) :: col_values(my_lt%n_cols)
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
325 type(
lt_t),
intent(in) :: my_lt
327 real(dp) :: col_values(my_lt%n_cols)
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)
335 type(
lt_t),
intent(in) :: my_lt
336 integer,
intent(in) :: col_ix
338 real(dp) :: col_value
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)
346 type(
lt_t),
intent(in) :: my_lt
347 real(dp),
intent(out) :: x_data(:), cols_rows(:, :)
349 x_data =
lt_get_xdata(my_lt%x_min, my_lt%dx, my_lt%n_points)
350 cols_rows = my_lt%cols_rows
360 real(dp),
intent(in) :: x_list(:), y_list(:)
361 real(dp),
intent(in) :: x_value
362 real(dp),
intent(out) :: y_value
364 integer :: ix, imin, imax
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)
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)
383 type(
lt_t),
intent(in) :: my_lt
384 character(len=*),
intent(in) :: filename
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
397 type(
lt_t),
intent(inout) :: my_lt
398 character(len=*),
intent(in) :: filename
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
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))
409 read(my_unit) my_lt%cols_rows
410 my_lt%rows_cols = transpose(my_lt%cols_rows)
418 function lt2_create(x_min, x_max, n_points, n_cols)
result(my_lt)
419 real(dp),
intent(in) :: x_min(2)
420 real(dp),
intent(in) :: x_max(2)
421 integer,
intent(in) :: n_points(2)
422 integer,
intent(in) :: n_cols
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"
428 my_lt%n_points = n_points
430 my_lt%dx = (x_max - x_min) / (n_points - 1)
431 my_lt%inv_dx = 1 / my_lt%dx
433 allocate(my_lt%rows_cols(n_points(1), n_points(2), n_cols))
435 my_lt%n_cols = n_cols
440 real(dp),
intent(in) :: x_min(2)
441 real(dp),
intent(in) :: x_max(2)
442 real(dp),
intent(in) :: spaced_data(:, :, :)
443 integer :: n_points(2), n_cols
446 n_points(1) =
size(spaced_data, 1)
447 n_points(2) =
size(spaced_data, 2)
448 n_cols =
size(spaced_data, 3)
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"
453 my_lt%n_cols = n_cols
454 my_lt%n_points = n_points
456 my_lt%dx = (x_max - x_min) / (n_points - 1)
457 my_lt%inv_dx = 1 / my_lt%dx
459 my_lt%rows_cols = spaced_data
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(:)
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)))
481 do ix = 1, my_lt%n_points(1)
482 my_lt%rows_cols(ix, :, col_ix) = &
490 type(
lt2_t),
intent(in) :: my_lt
491 real(dp),
intent(in) :: x1, x2
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
500 where (my_loc%low_ix < 1)
505 where (my_loc%low_ix >= my_lt%n_points)
506 my_loc%low_ix = my_lt%n_points - 1
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
525 type(
lt2_t),
intent(in) :: my_lt
526 integer,
intent(in) :: col_ix
530 real(dp) :: col_value
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))
539 col_value = sum(w * my_lt%rows_cols(ix(1):ix(1)+1, &
540 ix(2):ix(2)+1, col_ix))
546 function lt3_create(x_min, x_max, n_points, n_cols)
result(my_lt)
547 real(dp),
intent(in) :: x_min(3)
548 real(dp),
intent(in) :: x_max(3)
549 integer,
intent(in) :: n_points(3)
550 integer,
intent(in) :: n_cols
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"
556 my_lt%n_points = n_points
558 my_lt%dx = (x_max - x_min) / (n_points - 1)
559 my_lt%inv_dx = 1 / my_lt%dx
561 allocate(my_lt%rows_cols(n_points(1), n_points(2), n_points(3), n_cols))
563 my_lt%n_cols = n_cols
568 real(dp),
intent(in) :: x_min(3)
569 real(dp),
intent(in) :: x_max(3)
570 real(dp),
intent(in) :: spaced_data(:, :, :, :)
571 integer :: n_points(3), n_cols
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)
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"
582 my_lt%n_cols = n_cols
583 my_lt%n_points = n_points
585 my_lt%dx = (x_max - x_min) / (n_points - 1)
586 my_lt%inv_dx = 1 / my_lt%dx
588 allocate(my_lt%rows_cols(n_points(1), n_points(2), n_points(3), n_cols))
589 my_lt%rows_cols = spaced_data
594 type(
lt3_t),
intent(in) :: my_lt
595 real(dp),
intent(in) :: x1, x2, x3
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
604 where (my_loc%low_ix < 1)
609 where (my_loc%low_ix >= my_lt%n_points)
610 my_loc%low_ix = my_lt%n_points - 1
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
629 type(
lt3_t),
intent(in) :: my_lt
630 integer,
intent(in) :: col_ix
633 real(dp) :: w(2, 2, 2)
634 real(dp) :: col_value
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)
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))
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...