How to use the lenstronomy.LensModel.lens_model.LensModel function in lenstronomy

def test_sis_ray_tracing(self):
        z_source = 1.5
        lens_model_list = ['SIS']
        redshift_list = [0.5]
        lensModelMutli = MultiPlane(z_source=z_source, lens_model_list=lens_model_list, lens_redshift_list=redshift_list)
        lensModel = LensModel(lens_model_list=lens_model_list)
        kwargs_lens = [{'theta_E': 1, 'center_x': 0, 'center_y': 0}]
        beta_x_simple, beta_y_simple = lensModel.ray_shooting(1, 0, kwargs_lens)
        beta_x_multi, beta_y_multi = lensModelMutli.ray_shooting(1, 0, kwargs_lens)
        npt.assert_almost_equal(beta_x_simple, beta_x_multi, decimal=10)
        npt.assert_almost_equal(beta_y_simple, beta_y_multi, decimal=10)
        npt.assert_almost_equal(beta_x_simple, 0, decimal=10)
        npt.assert_almost_equal(beta_y_simple, 0, decimal=10)

def test_radial_tangential_distortions(self):
        lens_model_list = ['CURVED_ARC', 'SHEAR', 'FLEXION']
        center_x, center_y = 0, 0
        curvature = 1./2
        lens = LensModel(lens_model_list=lens_model_list)
        kwargs_lens = [{'tangential_stretch': 10, 'radial_stretch': 1., 'curvature': curvature,
                        'direction': -10, 'center_x': center_x, 'center_y': center_y},
                       {'gamma1': -0., 'gamma2': -0.0},
                       {'g1': 0., 'g2': 0., 'g3': -0., 'g4': 0}]

        extensions = LensModelExtensions(lensModel=lens)

        lambda_rad, lambda_tan, orientation_angle, dlambda_tan_dtan, dlambda_tan_drad, dlambda_rad_drad, dlambda_rad_dtan, dphi_tan_dtan, dphi_tan_drad, dphi_rad_drad, dphi_rad_dtan = extensions.radial_tangential_differentials(
            x=center_x, y=center_y, kwargs_lens=kwargs_lens, smoothing_3rd=0.0001)

        l = 1. / dphi_tan_dtan
        npt.assert_almost_equal(l, 1./curvature)

{'alpha_Rs': 0.003, 'center_y': -0.1, 'center_x': -1, 'Rs': 0.13},
                          {'alpha_Rs': 0.0008, 'center_y': 0.42, 'center_x': 0.1, 'Rs': 0.13},
                          {'theta_E': 0.3, 'center_y': 0.7, 'center_x': 1.08},
                          {'alpha_Rs': 0.006, 'center_y': 0.5, 'center_x': 0.75, 'Rs': 0.16}]

        self.convention_idx = [14,18]
        #self.convention_idx = False

        lens_model_list_full = lens_model_list_simple + front_halos + main_halos + back_halos
        redshift_list_full = redshift_list_simple + front_redshifts + main_redshifts + back_redshifts
        self.kwargs_lens_full = self.kwargs_lens_simple + self.front_args + self.main_args + self.back_args

        self.lens_model_full = LensModel(lens_model_list_full, z_source=1.5, lens_redshift_list=redshift_list_full,
                                         cosmo=self.cosmo,
                                         multi_plane=True)
        self.lens_model_full_convention = LensModel(lens_model_list_full, z_source=1.5, lens_redshift_list=redshift_list_full,
                                         cosmo=self.cosmo,
                                         multi_plane=True, observed_convention_index=self.convention_idx)

        self.lens_model_front = LensModel(front_halos + main_halos, lens_redshift_list=front_redshifts + main_redshifts,
                                          z_source=1.5,
                                          cosmo=self.cosmo, multi_plane=True)
        self.kwargs_front = self.front_args + self.main_args

        self.lens_model_simple = LensModel(lens_model_list_simple, z_source=1.5, lens_redshift_list=redshift_list_simple,
                                           cosmo=self.cosmo,
                                           multi_plane=True)

        self.optimizer_simple = Optimizer(self.x_pos_simple, self.y_pos_simple,
                                          magnification_target=self.magnification_simple,
                                          redshift_list=redshift_list_simple, simplex_n_iterations=2,
                                          lens_model_list=lens_model_list_simple, kwargs_lens=self.kwargs_lens_simple,

def test_arcs_at_image_position(self):
        # lensing quantities
        kwargs_spp = {'theta_E': 1.26, 'gamma': 2., 'e1': 0.1, 'e2': -0.1, 'center_x': 0.0, 'center_y': 0.0}  # parameters of the deflector lens model

        # the lens model is a supperposition of an elliptical lens model with external shear
        lens_model_list = ['SPEP']  #, 'SHEAR']
        kwargs_lens = [kwargs_spp]  #, kwargs_shear]
        lens_model_class = LensModel(lens_model_list=lens_model_list)
        lensEquationSolver = LensEquationSolver(lens_model_class)
        source_x = 0.
        source_y = 0.05
        x_image, y_image = lensEquationSolver.findBrightImage(source_x, source_y, kwargs_lens, numImages=4,
                                                              min_distance=0.05, search_window=5)
        arc_model = LensModel(lens_model_list=['CURVED_ARC', 'SHIFT'])
        for i in range(len(x_image)):
            x0, y0 = x_image[i], y_image[i]
            print(x0, y0, i)
            ext = LensModelExtensions(lensModel=lens_model_class)
            kwargs_arc_i = ext.curved_arc_estimate(x0, y0, kwargs_lens)
            alpha_x, alpha_y = lens_model_class.alpha(x0, y0, kwargs_lens)
            kwargs_arc = [kwargs_arc_i, {'alpha_x': alpha_x, 'alpha_y': alpha_y}]
            print(kwargs_arc_i)
            direction = kwargs_arc_i['direction']
            print(np.cos(direction), np.sin(direction))

deltaPix = 0.1  # pixel size in arcsec (area per pixel = deltaPix**2)
        fwhm = 0.5  # full width half max of PSF

        # PSF specification

        kwargs_data = sim_util.data_configure_simple(numPix, deltaPix, exp_time, sigma_bkg)
        data_class = ImageData(**kwargs_data)
        kwargs_psf_gaussian = {'psf_type': 'GAUSSIAN', 'fwhm': fwhm, 'pixel_size': deltaPix}
        psf = PSF(**kwargs_psf_gaussian)
        kwargs_psf = {'psf_type': 'PIXEL', 'kernel_point_source': psf.kernel_point_source}
        psf_class = PSF(**kwargs_psf)
        kwargs_spemd = {'theta_E': 1., 'gamma': 1.8, 'center_x': 0, 'center_y': 0, 'e1': 0.1, 'e2': 0.1}

        lens_model_list = ['SPEP']
        self.kwargs_lens = [kwargs_spemd]
        lens_model_class = LensModel(lens_model_list=lens_model_list)
        kwargs_sersic = {'amp': 1., 'R_sersic': 0.1, 'n_sersic': 2, 'center_x': 0, 'center_y': 0}
        # 'SERSIC_ELLIPSE': elliptical Sersic profile
        kwargs_sersic_ellipse = {'amp': 1., 'R_sersic': .6, 'n_sersic': 3, 'center_x': 0, 'center_y': 0,
                                 'e1': 0.1, 'e2': 0.1}

        lens_light_model_list = ['SERSIC']
        self.kwargs_lens_light = [kwargs_sersic]
        lens_light_model_class = LightModel(light_model_list=lens_light_model_list)
        source_model_list = ['SERSIC_ELLIPSE']
        self.kwargs_source = [kwargs_sersic_ellipse]
        source_model_class = LightModel(light_model_list=source_model_list)

        kwargs_numerics = {'supersampling_factor': 1, 'supersampling_convolution': False, 'compute_mode': 'regular'}
        imageModel = ImageModel(data_class, psf_class, lens_model_class, source_model_class,
                                lens_light_model_class, kwargs_numerics=kwargs_numerics)
        image_sim = sim_util.simulate_simple(imageModel, self.kwargs_lens, self.kwargs_source,

psf_fwhm = 0.7  # Gaussian FWHM psf
        kwargs_cosmo = {'d_d': 1000, 'd_s': 1500, 'd_ds': 800}

        # light profile
        light_profile_list = ['HERNQUIST']
        r_eff = 1.8
        kwargs_light = [{'Rs':  r_eff, 'amp': 1.}]  # effective half light radius (2d projected) in arcsec

        # mass profile
        mass_profile_list = ['SPP']
        theta_E = 1.2
        gamma = 2.
        kwargs_profile = [{'theta_E': theta_E, 'gamma': gamma}]  # Einstein radius (arcsec) and power-law slope

        # mge of lens profile
        lensModel = LensModel(mass_profile_list)
        r_array = np.logspace(-2, 2, 100)*theta_E
        kappa_r = lensModel.kappa(r_array, 0, kwargs_profile)
        amps, sigmas, norm = mge.mge_1d(r_array, kappa_r, N=20)
        mass_profile_list_mge = ['MULTI_GAUSSIAN_KAPPA']
        kwargs_profile_mge = [{'amp': amps, 'sigma': sigmas}]
        kwargs_psf = {'psf_type': 'GAUSSIAN', 'fwhm': psf_fwhm}
        kwargs_model = {'mass_profile_list': mass_profile_list,
                        'light_profile_list': light_profile_list,
                        'anisotropy_model': anisotropy_type}
        kwargs_numerics = {'interpol_grid_num': 100, 'log_integration': True,
                           'max_integrate': 100, 'min_integrate': 0.01}

        galkin = Galkin(kwargs_model=kwargs_model, kwargs_psf=kwargs_psf, kwargs_cosmo=kwargs_cosmo,
                        kwargs_aperture=kwargs_aperture, kwargs_numerics=kwargs_numerics)
        sigma_v = galkin.dispersion(kwargs_profile, kwargs_light, kwargs_anisotropy)

Rs = 10.
        alpha_Rs = 10
        x = np.linspace(Rs, Rs, 1000)
        y = np.linspace(0.2 * Rs, 2 * Rs, 1000)
        center_x, center_y = -1.2, 0.46

        zlist_single = [0.5, 0.5]
        zlist_multi = [0.5, 0.8]
        zlist = [zlist_single, zlist_multi]
        numerical_alpha_class = TestClass()

        for i, flag in enumerate([False, True]):
            lensmodel = LensModel(lens_model_list=['NumericalAlpha', 'NFW'], z_source=1.5, z_lens=0.5, lens_redshift_list=zlist[i],
                                  multi_plane=flag, numerical_alpha_class=numerical_alpha_class)

            lensmodel_nfw = LensModel(lens_model_list=['NFW', 'NFW'], z_source=1.5, z_lens=0.5, lens_redshift_list=zlist[i],
                                  multi_plane=flag, numerical_alpha_class=numerical_alpha_class)

            keywords_num = [{'norm': alpha_Rs, 'Rs': Rs, 'center_x': center_x, 'center_y': center_y},
                            {'alpha_Rs': 0.7*alpha_Rs, 'Rs': 2*Rs, 'center_x': center_x, 'center_y': center_y}]
            keywords_nfw = [{'alpha_Rs': alpha_Rs, 'Rs': Rs, 'center_x': center_x, 'center_y': center_y},
                            {'alpha_Rs': 0.7 * alpha_Rs, 'Rs': 2 * Rs, 'center_x': center_x, 'center_y': center_y}]

            dx, dy = lensmodel.alpha(x, y, keywords_num)
            dxnfw, dynfw = lensmodel_nfw.alpha(x, y, keywords_nfw)
            npt.assert_almost_equal(dx, dxnfw)
            npt.assert_almost_equal(dy, dynfw)

:param y_values: 1d array of y coordinates to interpolate
        (e.g. np.linspace(ymin,ymax,steps))
        :param interp_lensmodel: lensmodel to interpolate
        :param interp_args: kwargs for interp_lensmodel
        :return: interpolated lensmodel
        """
        xx, yy = np.meshgrid(x_values, y_values)
        L = int(len(x_values))
        xx, yy = xx.ravel(), yy.ravel()

        f_x, f_y = interp_lensmodel.alpha(xx, yy, interp_args)

        interp_args = [{'f_x': f_x.reshape(L, L), 'f_y': f_y.reshape(L, L),
                        'grid_interp_x': x_values, 'grid_interp_y': y_values}]

        return LensModel(lens_model_list=['INTERPOL'], redshift_list=[self._z_background], cosmo=self.astropy_instance,
                         z_source=self.z_source, multi_plane=True), interp_args

else:
        lens_model_list_i = [lens_model_list[k] for k in index_lens_model_list[band_index]]
        if lens_redshift_list is not None:
            lens_redshift_list_i = [lens_redshift_list[k] for k in index_lens_model_list[band_index]]
        else:
            lens_redshift_list_i = lens_redshift_list
        if observed_convention_index is not None:
            counter = 0
            observed_convention_index_i = []
            for k in index_lens_model_list[band_index]:
                if k in observed_convention_index:
                    observed_convention_index_i.append(counter)
                counter += 1
        else:
            observed_convention_index_i = observed_convention_index
    lens_model_class = LensModel(lens_model_list=lens_model_list_i, z_lens=z_lens, z_source=z_source,
                                 lens_redshift_list=lens_redshift_list_i,
                                 multi_plane=multi_plane, cosmo=cosmo,
                                 observed_convention_index=observed_convention_index_i)

    if index_source_light_model_list is None or all_models is True:
        source_light_model_list_i = source_light_model_list
        source_deflection_scaling_list_i = source_deflection_scaling_list
        source_redshift_list_i = source_redshift_list
    else:
        source_light_model_list_i = [source_light_model_list[k] for k in index_source_light_model_list[band_index]]
        if source_deflection_scaling_list is None:
            source_deflection_scaling_list_i = source_deflection_scaling_list
        else:
            source_deflection_scaling_list_i = [source_deflection_scaling_list[k] for k in index_source_light_model_list[band_index]]
        if source_redshift_list is None:
            source_redshift_list_i = source_redshift_list

halo_args.append(lensmodel_args[i])
                if observed_convention_inds is not None:
                    if i in observed_convention_inds: convention_inds.append(count)
                count += 1

            else:

                macro_names.append(lensmodel.lens_model_list[i])
                macro_redshifts.append(z)
                macro_args.append(lensmodel_args[i])

        macromodel = LensModel(lens_model_list=macro_names, lens_redshift_list=macro_redshifts, cosmo=self._astropy_instance,
                               multi_plane=True,
                               z_source=self._z_source, numerical_alpha_class=numerical_alpha_class)

        halo_lensmodel = LensModel(lens_model_list=halo_names, lens_redshift_list=halo_redshifts,
                                   cosmo=self._astropy_instance, multi_plane=True, z_source=self._z_source,
                                   numerical_alpha_class = numerical_alpha_class, observed_convention_index=convention_inds)

        return macromodel, macro_args, halo_lensmodel, halo_args