docs: def velocity_dispersion_numerical(self, kwargs_lens, kwargs_lens_light, kwargs_anisotropy, kwargs_aperture, psf_fwhm,
                                        aperture_type, anisotropy_model, r_eff=None, kwargs_numerics={}, MGE_light=False,
                                        MGE_mass=False, lens_model_kinematics_bool=None, light_model_kinematics_bool=None,
                                        Hernquist_approx=False):
          """
          Computes the LOS velocity dispersion of the deflector galaxy with arbitrary combinations of light and mass models.
          For a detailed description, visit the description of the Galkin() class.
          Additionaly to executing the Galkin routine, it has an optional Multi-Gaussian-Expansion decomposition of lens
          and light models that do not have a three-dimensional distribution built in, such as Sersic profiles etc.
  
          The center of all the lens and lens light models that are part of the kinematic estimate must be centered on the
          same point.
  
          :param kwargs_lens: lens model parameters
          :param kwargs_lens_light: lens light parameters
          :param kwargs_anisotropy: anisotropy parameters (see Galkin module)
          :param kwargs_aperture: aperture parameters (see Galkin module)
          :param psf_fwhm: full width at half maximum of the seeing (Gaussian form)
          :param aperture_type: type of aperture (see Galkin module
          :param anisotropy_model: stellar anisotropy model (see Galkin module)
          :param r_eff: a rough estimate of the half light radius of the lens light in case of computing the MGE of the
           light profile
          :param kwargs_numerics: keyword arguments that contain numerical options (see Galkin module)
          :param MGE_light: bool, if true performs the MGE for the light distribution
          :param MGE_mass: bool, if true performs the MGE for the mass distribution
          :param lens_model_kinematics_bool: bool list of length of the lens model. Only takes a subset of all the models
              as part of the kinematics computation (can be used to ignore substructure, shear etc that do not describe the
              main deflector potential
          :param light_model_kinematics_bool: bool list of length of the light model. Only takes a subset of all the models
              as part of the kinematics computation (can be used to ignore light components that do not describe the main
              deflector
          :return: LOS velocity dispersion [km/s]
          """
  
          kwargs_cosmo = {'D_d': self.lensCosmo.D_d, 'D_s': self.lensCosmo.D_s, 'D_ds': self.lensCosmo.D_ds}
          mass_profile_list = []
          kwargs_profile = []
          if lens_model_kinematics_bool is None:
              lens_model_kinematics_bool = [True] * len(kwargs_lens)
          for i, lens_model in enumerate(self.kwargs_options['lens_model_list']):
              if lens_model_kinematics_bool[i] is True:
                  mass_profile_list.append(lens_model)
                  if lens_model in ['INTERPOL', 'INTERPOL_SCLAED']:
                      center_x, center_y = self._lensModelExt.lens_center(kwargs_lens, k=i)
                      kwargs_lens_i = copy.deepcopy(kwargs_lens[i])
                      kwargs_lens_i['grid_interp_x'] -= center_x
                      kwargs_lens_i['grid_interp_y'] -= center_y
                  else:
                      kwargs_lens_i = {k: v for k, v in kwargs_lens[i].items() if not k in ['center_x', 'center_y']}
                  kwargs_profile.append(kwargs_lens_i)
  
          if MGE_mass is True:
              lensModel = LensModel(lens_model_list=mass_profile_list)
              massModel = LensModelExtensions(lensModel)
              theta_E = massModel.effective_einstein_radius(kwargs_profile)
              r_array = np.logspace(-4, 2, 200) * theta_E
              mass_r = lensModel.kappa(r_array, np.zeros_like(r_array), kwargs_profile)
              amps, sigmas, norm = mge.mge_1d(r_array, mass_r, N=20)
              mass_profile_list = ['MULTI_GAUSSIAN_KAPPA']
              kwargs_profile = [{'amp': amps, 'sigma': sigmas}]
  
          light_profile_list = []
          kwargs_light = []
          if light_model_kinematics_bool is None:
              light_model_kinematics_bool = [True] * len(kwargs_lens_light)
          for i, light_model in enumerate(self.kwargs_options['lens_light_model_list']):
              if light_model_kinematics_bool[i]:
                  light_profile_list.append(light_model)
                  kwargs_lens_light_i = {k: v for k, v in kwargs_lens_light[i].items() if not k in ['center_x', 'center_y']}
                  if 'q' in kwargs_lens_light_i:
                      kwargs_lens_light_i['q'] = 1
                  kwargs_light.append(kwargs_lens_light_i)
          if r_eff is None:
              lensAnalysis = LensAnalysis({'lens_light_model_list': light_profile_list})
              r_eff = lensAnalysis.half_light_radius_lens(kwargs_light, model_bool_list=light_model_kinematics_bool)
          if Hernquist_approx is True:
              light_profile_list = ['HERNQUIST']
              kwargs_light = [{'Rs':  r_eff, 'amp': 1.}]
          else:
              if MGE_light is True:
                  lightModel = LightModel(light_profile_list)
                  r_array = np.logspace(-3, 2, 200) * r_eff * 2
                  flux_r = lightModel.surface_brightness(r_array, 0, kwargs_light)
                  amps, sigmas, norm = mge.mge_1d(r_array, flux_r, N=20)
                  light_profile_list = ['MULTI_GAUSSIAN']
                  kwargs_light = [{'amp': amps, 'sigma': sigmas}]
  
          galkin = Galkin(mass_profile_list, light_profile_list, aperture_type=aperture_type,
                          anisotropy_model=anisotropy_model, fwhm=psf_fwhm, kwargs_cosmo=kwargs_cosmo, **kwargs_numerics)
          sigma2 = galkin.vel_disp(kwargs_profile, kwargs_light, kwargs_anisotropy, kwargs_aperture)
          return sigma2