code: def interpolate_1d(x, xp, *args, **kwargs):
      r
      fill_value = kwargs.pop('fill_value', np.nan)
      axis = kwargs.pop('axis', 0)
      x = np.asanyarray(x).reshape(-1)
      ndim = xp.ndim
      sort_args = np.argsort(xp, axis=axis)
      sort_x = np.argsort(x)
      sorter = broadcast_indices(xp, sort_args, ndim, axis)
      xp = xp[sorter]
      variables = [arr[sorter] for arr in args]
      x_array = x[sort_x]
      expand = [np.newaxis] * ndim
      expand[axis] = slice(None)
      x_array = x_array[tuple(expand)]
      minv = np.apply_along_axis(np.searchsorted, axis, xp, x[sort_x])
      minv2 = np.copy(minv)
      if ((np.max(minv) == xp.shape[axis]) or (np.min(minv) == 0)) and fill_value is None:
          raise ValueError('Interpolation point out of data bounds encountered')
      if np.max(minv) == xp.shape[axis]:
          warnings.warn('Interpolation point out of data bounds encountered')
          minv2[minv == xp.shape[axis]] = xp.shape[axis] - 1
      if np.min(minv) == 0:
          minv2[minv == 0] = 1
      above = broadcast_indices(xp, minv2, ndim, axis)
      below = broadcast_indices(xp, minv2 - 1, ndim, axis)
      if np.any(x_array < xp[below]):
          warnings.warn('Interpolation point out of data bounds encountered')
      ret = []
      for var in variables:
          var_interp = var[below] + (var[above] - var[below]) * ((x_array - xp[below])
                                                                 / (xp[above] - xp[below]))
          var_interp[minv == xp.shape[axis]] = fill_value
          var_interp[x_array < xp[below]] = fill_value
          if x[0] > x[-1]:
              var_interp = np.swapaxes(np.swapaxes(var_interp, 0, axis)[::-1], 0, axis)
          ret.append(var_interp)
      if len(ret) == 1:
          return ret[0]
      else:
          return ret