module int1 { # generic 1-d interpolation integer :: nd, nf integer, dimension (nd), allocatable :: x1, x2, y1, y2 logical, dimension (nd), allocatable :: m real, dimension (nf,nd), allocatable :: w #% _init (coord, o1, d1, n1, interp, nd, nf, weight, ww) real, dimension (:), intent (in) :: coord, weight real, dimension (:,:), intent (out) :: ww real, intent (in) :: o1, d1 integer, intent (in) :: n1 optional :: weight, ww interface { integer function interp (x, w) { real, intent (in) :: x real, dimension (:) :: w } } integer :: id, ix, i1, i2, stat, j real :: rx if (present (ww)) ww = 0. do id = 1, nd { rx = (coord (id) - o1)/d1 ix = rx - 0.5*nf + 1. i1 = rx; rx -= i1 i1 = ix + 1 ; i2 = ix + nf if (i2 < 1 .or. i1 > n1) {m (id) = .true. ; w (:, id) = 0.; cycle} m (id) = .false. ; stat = interp (rx, w (:,id)) x1 (id) = max (i1,1); y1 (id) = x1 (id) - ix x2 (id) = min (i2,n1); y2 (id) = x2 (id) - ix if (y1 (id) > 1) w (1:y1(id)-1,id) = 0. if (y2 (id) < nf) w (y2(id)+1:nf,id) = 0. if (present (weight)) w (:,id) = w (:,id) * weight (id) if (present (ww)) do j = 1, nf ww (i1:i2+1-j,j) += w (1:1+nf-j,id) * w (j:nf,id) } #% _lop (mdl, ord) integer :: i do i = 1, nd { if (m (i)) cycle if (adj) mdl (x1(i):x2(i)) += w (y1(i):y2(i),i) * ord (i) else ord (i) += sum (mdl (x1(i):x2(i)) * w (y1(i):y2(i),i)) } #% _close }