code: def rfiw_phantom(gtab, snr=None, rng=None):
      """rectangle fiber immersed in water"""
      # define voxel index
      slice_ind = np.zeros((10, 10, 8))
      slice_ind[4:7, 4:7, :] = 1
      slice_ind[4:7, 7, :] = 2
      slice_ind[7, 7, :] = 3
      slice_ind[7, 4:7, :] = 4
      slice_ind[7, 3, :] = 5
      slice_ind[4:7, 3, :] = 6
      slice_ind[3, 3, :] = 7
      slice_ind[3, 4:7, :] = 8
      slice_ind[3, 7, :] = 9
  
      # Define tissue diffusion parameters
      # Restricted diffusion
      ADr = 0.99e-3
      RDr = 0.0
      # Hindered diffusion
      ADh = 2.26e-3
      RDh = 0.87e-3
      # S0 value for tissue
      S1 = 50
      # Fraction between Restricted and Hindered diffusion
      fia = 0.51
  
      # Define water diffusion
      Dwater = 3e-3
      S2 = 100  # S0 value for water
  
      # Define tissue volume fraction for each voxel type (in index order)
      f = np.array(
          [0.0, 1.0, 0.6, 0.18, 0.30, 0.15, 0.50, 0.35, 0.70, 0.42]
      )
  
      # Define S0 for each voxel (in index order)
      S0 = S1 * f + S2 * (1 - f)
  
      # multi tensor simulations assume that each water pull as constant S0
      # since I am assuming that tissue and water voxels have different S0,
      # tissue volume fractions have to be adjusted to the measured f values when
      # constant S0 are assumed constant. Doing this correction, simulations will
      # be analogous to simulates that S0 are different for each media. (For more
      # details on this contact the phantom designer)
      f1 = f * S1 / S0
  
      mevals = np.array(
          [[ADr, RDr, RDr], [ADh, RDh, RDh], [Dwater, Dwater, Dwater]]
      )
      angles = [(0, 0, 1), (0, 0, 1), (0, 0, 1)]
      DWI = np.zeros(slice_ind.shape + (gtab.bvals.size,))
      for i in range(10):
          fractions = [
              f1[i] * fia * 100,
              f1[i] * (1 - fia) * 100,
              (1 - f1[i]) * 100,
          ]
          sig, direction = multi_tensor(
              gtab,
              mevals,
              S0=S0[i],
              angles=angles,
              fractions=fractions,
              snr=None,
          )
          DWI[slice_ind == i, :] = sig
  
      if snr is None:
          return DWI
      else:
          if rng is None:
              rng = np.random.default_rng()
          sigma = S2 * 1.0 / snr
          n1 = rng.normal(0, sigma, size=DWI.shape)
          n2 = rng.normal(0, sigma, size=DWI.shape)
          return [
              np.sqrt(
                  (DWI / np.sqrt(2) + n1) ** 2
                  + (DWI / np.sqrt(2) + n2) ** 2
              ),
              sigma,
          ]
  
imports: import warnings
  import numpy
  import scipy.special
  from numpy.testing import assert_, assert_equal, assert_raises, assert_array_almost_equal, assert_warns
  from dipy.denoise.localpca import dimensionality_problem_message, create_patch_radius_arr, compute_patch_size, compute_num_samples, compute_suggested_patch_radius, localpca, mppca, genpca, _pca_classifier
  from dipy.sims.voxel import multi_tensor
  from dipy.core.gradients import gradient_table, generate_bvecs
  from dipy.testing.decorators import set_random_number_generator
label_desc: Use image node_id 2 for calling a global function with example usage: rfiw_phantom(gtab, snr, rng) and returns: DWI, list