<filename>Dans_Diffraction/classes_scattering.py
# -*- coding: utf-8 -*-
"""
Scattering Class "classes_scattering.py"
    Scattering functions for Crystal class.

By <NAME>, PhD
Diamond
2017

Version 2.0.1
Last updated: 28/09/21

Version History:
10/09/17 0.1    Program created
30/10/17 1.0    Main functions finshed, some testing done
06/01/18 1.1    Renamed classes_scattering.py
31/10/18 1.2    Added print_symmetric_contributions
21/01/19 1.3    Added non-resonant diffraction, corrected resonant diffraction
16/12/19 1.4    Added multiple_scattering code, print_all_reflections updated with units
18/02/20 1.5    Added tensor_scattering code
20/03/20 1.6    Increased powder gauss width from 2fwhm to 6fwhm, added powder averaging
14/04/20 1.6    Added powder_correction
26/05/20 1.7    Removed tensor_scattering
16/06/20 1.7.1  Added output option of setup_scatter
04/01/21 1.7.2  Added structure_factor function
21/01/21 1.8    Added xray_dispersion scattering function
10/06/21 1.9    Added x_ray calculation using Waasmaier and Kirfel scattering factors.
09/07/21 1.9    Added new scattering factors as option on normal scattering functions
20/08/21 2.0    Switched over to new scattering module, added self.powder()
28/09/21 2.0.1  Added __repr__

@author: DGPorter
"""

import sys,os
import numpy as np
import datetime

from . import functions_general as fg
from . import functions_crystallography as fc
from . import functions_scattering as fs
from . import multiple_scattering as ms
# from . import tensor_scattering as ts  # Removed V1.7

__version__ = '2.0.1'
__scattering_types__ = {'xray': ['xray','x','x-ray','thomson','charge'],
                        'neutron': ['neutron','n','nuclear'],
                        'xray magnetic': ['xray magnetic','magnetic xray','spin xray','xray spin'],
                        'neutron magnetic': ['neutron magnetic','magnetic neutron','magnetic'],
                        'xray resonant': ['xray resonant','resonant','resonant xray','rxs'],
                        'xray dispersion': ['dispersion', 'xray dispersion']}


class Scattering:
    """
    Simulate diffraction from Crystal
    Useage:
        xtl = Crystal()
        xtl.Scatter.setup_scatter(type='x-ray',energy_keV=8.0)
        xtl.Scatter.intensity([h,k,l]) # Returns intensity
        print(xtl.Scatter.print_all_refelctions()) # Returns formated string of all allowed reflections
        
        Allowed radiation types:
            'xray','neutron','xray magnetic','neutron magnetic','xray resonant'
    """
    
    #------Options-------
    # Standard Options
    _hkl = None  # added so recalculation not required
    _scattering_type = 'xray'  # 'xray','neutron','xray magnetic','neutron magnetic','xray resonant'
    _scattering_specular_direction = [0,0,1]  # reflection
    _scattering_parallel_direction = [0,0,1]  # transmission
    _scattering_theta_offset = 0.0
    _scattering_min_theta = -180.0
    _scattering_max_theta = 180.0
    _scattering_min_twotheta = -180.0
    _scattering_max_twotheta = 180.0
    _time_report = True
    _debug_mode = False

    # powder options
    _powder_units = 'tth' # tth (two theta), Q, d
    _powder_background = 0.0
    _powder_peak_width = 0.01  # in Deg
    _powder_average = True
    _powder_pixels = 2000
    
    # Complex Structure factor
    _return_structure_factor = False

    # Uses the coefficients for analytical approximation to the scattering factors from:
    #        "<NAME>, Acta Cryst. (1995) A51, 416-431"
    _use_waaskirf_scattering_factor = False
    
    # Thermal Factors
    _use_isotropic_thermal_factor = True
    _use_anisotropic_thermal_factor = False
    
    # Magnetic Options
    _calclate_magnetic_component = True
    _use_magnetic_form_factor = True
    
    # Polarisation Options
    _polarised = False
    _polarisation = 'sp'
    _polarisation_vector_incident = [0, 1, 0]
    
    # Radiation energy
    _energy_kev = fg.Cu
    
    # Resonant X-ray Options
    _azimuthal_angle = 0
    _azimuthal_reference = [1,0,0]
    _resonant_flm = (0, 1, 0)
    _resonant_approximation_e1e1 = True
    _resonant_approximation_e2e2 = False
    _resonant_approximation_e1e2 = False
    _resonant_approximation_m1m1 = False
    
    def __init__(self, xtl):
        self.xtl = xtl

        # Initialise the scattering type container
        self.Type = ScatteringTypes(self, __scattering_types__)

    def __repr__(self):
        return 'Scattering(%s, %s)' % (self.xtl.name, self._scattering_type)

    def __str__(self):
        out = 'Scatter(%r)\n' % self.xtl
        out += '       Type: %s\n' % self._scattering_type
        out += '     Energy: %s keV\n' % self._energy_kev
        out += ' Wavelength: %s A\n' % fc.energy2wave(self._energy_kev)
        out += ' ---Settings---\n'
        out += '      Powder units: %s\n' % self._powder_units
        out += '    Isotropic ADPs: %s\n' % self._use_isotropic_thermal_factor
        #out += '  Anisotropic ADPs: %s\n' % self._use_anisotropic_thermal_factor
        out += '  Specular Direction (reflection): (%2.0f,%2.0f,%2.0f)\n' % (
            self._scattering_specular_direction[0], self._scattering_specular_direction[1],
            self._scattering_specular_direction[2])
        out += 'Parallel Direction (transmission): (%2.0f,%2.0f,%2.0f)\n' % (
            self._scattering_parallel_direction[0], self._scattering_parallel_direction[1],
            self._scattering_parallel_direction[2])
        out += '      theta offset: %s\n' % self._scattering_theta_offset
        out += '         min theta: %s\n' % self._scattering_min_theta
        out += '         max theta: %s\n' % self._scattering_max_theta
        out += '      min twotheta: %s\n' % self._scattering_min_twotheta
        out += '      max twotheta: %s\n' % self._scattering_max_twotheta
        out += ' ---X-Ray Settings---\n'
        out += ' Waasmaier scattering factor: %s\n' % self._use_waaskirf_scattering_factor
        out += ' ---Magnetic Settings---\n'
        out += '   Mag. scattering: %s\n' % self._calclate_magnetic_component
        out += '  Mag. form factor: %s\n' % self._use_magnetic_form_factor
        out += '         Polarised: %s\n' % self._polarised
        out += '      Polarisation: %s\n' % self._polarisation
        out += '       Pol. vector: (%s,%s,%s)\n' % (self._polarisation_vector_incident[0],
                                                     self._polarisation_vector_incident[1],
                                                     self._polarisation_vector_incident[2])
        out += ' ---Resonant Settings---\n'
        out += '   Azimuthal angle: %s\n' % self._azimuthal_angle
        out += '    Azimuthal ref.: (%s,%s,%s)\n' % (self._azimuthal_reference[0], self._azimuthal_reference[1],
                                                     self._azimuthal_reference[2])
        out += '               flm: (%s,%s,%s)\n' % (self._resonant_flm[0], self._resonant_flm[1],
                                                     self._resonant_flm[2])
        out += '  use e1e1 approx.: %s\n' % self._resonant_approximation_e1e1
        #out += '  use e2e2 approx.: %s\n' % self._resonant_approximation_e2e2
        #out += '  use e1e2 approx.: %s\n' % self._resonant_approximation_e1e2
        #out += '  use m1m1 approx.: %s\n' % self._resonant_approximation_m1m1
        return out

    def setup_scatter(self, scattering_type=None, energy_kev=None, wavelength_a=None, powder_units=None,
                      specular=None, parallel=None, theta_offset=None,
                      min_theta=None, max_theta=None, min_twotheta=None, max_twotheta=None,
                      output=True, scattering_factors=None, magnetic_formfactor=None,
                      polarisation=None, polarisation_vector=None, azimuthal_reference=None, azimuth=None, flm=None):
        """
        Simple way to set scattering parameters, each parameter is internal to xtl (self)

        scattering_type: self._scattering type            :  'xray','neutron','xray magnetic','neutron magnetic','xray resonant', 'xray dispersion'
        energy_kev  : self._energy_kev                    :  radiation energy in keV
        wavelength_a: self._wavelength_a                  :  radiation wavelength in Angstrom
        powder_units: self._powder_units                  :  units to use when displaying/ plotting ['twotheta', 'd',' 'q']
        min_twotheta: self._scattering_min_two_theta      :  minimum detector (two-theta) angle
        max_twotheta: self._scattering_max_two_theta      :  maximum detector (two-theta) angle
        min_theta   : self._scattering_min_theta          :  minimum sample angle = -opening angle
        max_theta   : self._scattering_max_theta          :  maximum sample angle = opening angle
        theta_offset: self._scattering_theta_offset       :  sample offset angle
        specular    : self._scattering_specular_direction : [h,k,l] : reflections normal to sample surface
        parallel    : self._scattering_parallel_direction : [h,k,l] : reflections normal to sample surface
        scattering_factors: self._use_waaskirf_scattering_factor : xray scattering factor ['waaskirf', 'itc']
        magnetic_formfactor: self._use_magnetic_form_factor: True/False magnetic form factor for magnetic SF
        polarisation: self._polarisation                  : beam polarisation setting ['ss', 'sp'*, 'sp', 'pp']
        polarisation_vector: _polarisation_vector_incident: [x,y,z] incident polarisation vector
        azimuthal_reference: self._azimuthal_reference    : [h,k,l] direction of azimuthal zero angle
        azimuth    : self._azimuthal_angle                : azimuthal angle in deg
        flm        : self._resonant_flm                   : Resonant settings (flm1, flm2, flm3)
        """

        if scattering_type is not None:
            self._scattering_type = scattering_type

        if energy_kev is not None:
            self._energy_kev = energy_kev

        if wavelength_a is not None:
            self._energy_kev = fc.wave2energy(wavelength_a)

        if powder_units is not None:
            self._powder_units = powder_units

        if specular is not None:
            self._scattering_specular_direction = specular

        if parallel is not None:
            self._scattering_parallel_direction = parallel

        if theta_offset is not None:
            self._scattering_theta_offset = theta_offset

        if min_theta is not None:
            self._scattering_min_theta = min_theta

        if max_theta is not None:
            self._scattering_max_theta = max_theta

        if min_twotheta is not None:
            self._scattering_min_twotheta = min_twotheta

        if max_twotheta is not None:
            self._scattering_max_twotheta = max_twotheta

        if scattering_factors is not None:
            if scattering_factors.lower() in ['ws', 'waaskirf', 'alternate', 'alt']:
                print('Using scattering factors from: "Waasmaier and Kirfel, Acta Cryst. (1995) A51, 416-431"')
                self._use_waaskirf_scattering_factor = True
            else:
                print('Using scattering factors from: International Tables of Crystallography Vol. C, Table 6.1.1.4')
                self._use_waaskirf_scattering_factor = False

        if magnetic_formfactor is not None:
            self._use_magnetic_form_factor = magnetic_formfactor

        if polarisation_vector is not None:
            self._polarisation_vector_incident = np.array(polarisation_vector, dtype=np.float).reshape(3)

        if polarisation is not None:
            self._polarisation = polarisation

        if azimuthal_reference is not None:
            self._azimuthal_reference = np.array(azimuthal_reference, dtype=np.float).reshape(3)

        if azimuth is not None:
            self._azimuthal_angle = azimuth

        if flm is not None:
            self._resonant_flm = np.array(flm).reshape(3)

        if output:
            print(self)
            # print('Scattering Options:')
            # print('                            Type : %s' % self._scattering_type)
            # print('                  Default Energy : %6.3f keV' % self._energy_kev)
            # print('                    Powder Units : %s' % self._powder_units)
            # print('  Specular Direction (reflection): (%2.0f,%2.0f,%2.0f)' % (
            #     self._scattering_specular_direction[0], self._scattering_specular_direction[1],
            #     self._scattering_specular_direction[2]))
            # print('Parallel Direction (transmission): (%2.0f,%2.0f,%2.0f)' % (
            #     self._scattering_parallel_direction[0], self._scattering_parallel_direction[1],
            #     self._scattering_parallel_direction[2]))
            # print('                   Sample Offset : %5.2f' % self._scattering_theta_offset)
            # print('             Minimum Theta angle : %5.2f' % self._scattering_min_theta)
            # print('             Maximum Theta angle : %5.2f' % self._scattering_max_theta)
            # print('         Minimum Two-Theta angle : %5.2f' % self._scattering_min_twotheta)
            # print('         Maximum Two-Theta angle : %5.2f' % self._scattering_max_twotheta)

    def get_hkl(self, regenerate=True, remove_symmetric=False, reflection=False, transmission=False,  **kwargs):
        """
        Return stored hkl or generate
        :param regenerate: if True, hkl list will be regenerated, if False - previous list will be returned
        :param remove_symmetric: generate only non-symmetric hkl values
        :param reflection: generate only reflections possible in reflection geometry
        :param transmission: generate only reflections possible in transmission geometry
        :param kwargs: additional options to pass to setup_scatter()
        :return: array
        """
        if not regenerate and self._hkl is not None:
            return self._hkl
        self.setup_scatter(output=False, **kwargs)

        en = self._energy_kev
        max_tth = self._scattering_max_twotheta

        hkl = self.xtl.Cell.all_hkl(en, max_tth)
        if remove_symmetric:
            hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)
        hkl = self.xtl.Cell.sort_hkl(hkl)[1:]  # remove (0, 0, 0)

        if reflection:
            tth = self.xtl.Cell.tth(hkl, en)
            hkl = hkl[tth > self._scattering_min_twotheta, :]
            tth = tth[tth > self._scattering_min_twotheta]
            theta = self.xtl.Cell.theta_reflection(hkl, en, self._scattering_specular_direction,
                                                   self._scattering_theta_offset)
            p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)
            p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))
            hkl = hkl[p1 * p2]
        elif transmission:
            tth = self.xtl.Cell.tth(hkl, en)
            hkl = hkl[tth > self._scattering_min_twotheta, :]
            tth = tth[tth > self._scattering_min_twotheta]
            theta = self.xtl.Cell.theta_transmission(hkl, en, self._scattering_parallel_direction)

            p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)
            p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))
            hkl = hkl[p1 * p2]
        self._hkl = hkl
        return self._hkl

    def _debug(self):
        """Toggle on debug mode"""
        if self._debug_mode:
            print('Debug mode: off')
            self._debug_mode = False
            fs.DEBUG_MODE = False
        else:
            print('Debug mode: on')
            self._debug_mode = True
            fs.DEBUG_MODE = True

    def structure_factor(self, hkl=None, scattering_type=None, **kwargs):
        """
        Calculate the structure factor at reflection indexes (h,k,l)

        Notes:
        - Uses x-ray atomic form factors, calculated from approximated tables in the ITC
        - This may be a little slow for large numbers of reflections, as it is not currently
         possible to use accelerated calculation methods in Jython.
        - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs
        - Crystal.scale is used to scale the complex structure factor, so the intensity is
         reduced by (Crystal.scale)^2
        - Testing against structure factors calculated by Vesta.exe is exactly the same when using Waasmaier structure factors.

        :param hkl: array[n,3] : reflection indexes (h, k, l)
        :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']
        :param kwargs: additional options to pass to scattering function
        :return: complex array[n] : structure factors
        """

        if hkl is None:
            hkl = self.get_hkl()
        if scattering_type is None:
            scattering_type = self._scattering_type

        hkl = np.asarray(np.rint(hkl), dtype=np.float).reshape([-1, 3])
        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()

        q = self.xtl.Cell.calculateQ(hkl)
        r = self.xtl.Cell.calculateR(uvw)

        moment = fc.euler_moment(mxmymz, self.xtl.Cell.UV())
        azi_ref_q = self.xtl.Cell.calculateQ(self._azimuthal_reference)

        if 'energy_kev' in kwargs:
            energy_kev = kwargs.pop('energy_kev')
        elif 'wavelength_a' in kwargs:
            energy_kev = fc.wave2energy(kwargs.pop('wavelength_a'))
        else:
            energy_kev = self._energy_kev
        energy_kev = np.asarray(energy_kev, dtype=np.float).reshape(-1)
        nenergy = len(energy_kev)

        if 'psi' in kwargs:
            psi = kwargs.pop('psi')
        else:
            psi = self._azimuthal_angle
        psi = np.asarray(psi, dtype=np.float).reshape(-1)
        npsi = len(psi)

        options = fs.options(
            occ=occ,
            moment=moment,
            incident_polarisation_vector=self._polarisation_vector_incident,
            polarisation=self._polarisation,
            azi_ref_q=azi_ref_q,
            f0=self._resonant_flm[0],
            f1=self._resonant_flm[1],
            f2=self._resonant_flm[2],
        )
        options.update(kwargs)
        scattering_fun = fs.get_scattering_function(scattering_type)
        # print('Scattering function: %s' % scattering_fun.__name__)

        # Break up long lists of HKLs
        nref, natom = len(q), len(r)
        n_arrays = np.ceil(nref * natom / fs.MAX_QR_ARRAY)
        if n_arrays > 1:
            print('Splitting %d reflections (%d atoms) into %1.0f parts' % (nref, natom, n_arrays))
        q_array = np.array_split(q, n_arrays)
        sf = np.zeros([nref, nenergy, npsi], dtype=np.complex)
        start_time = datetime.datetime.now()
        for e, enval in enumerate(energy_kev):
            for p, psival in enumerate(psi):
                ls = 0
                for n, _q in enumerate(q_array):
                    if n_arrays > 1:
                        print(' Starting %2.0f/%2.0f: %d:%d' % (n + 1, n_arrays, ls, ls + len(_q)))
                    qmag = fg.mag(_q)
                    # Scattering factors
                    if scattering_type in fs.SCATTERING_TYPES['neutron']:
                        ff = fc.atom_properties(atom_type, 'Coh_b')
                    elif scattering_type in fs.SCATTERING_TYPES['xray fast']:
                        ff = fc.atom_properties(atom_type, 'Z')
                    elif scattering_type in fs.SCATTERING_TYPES['xray dispersion']:
                        ff = fc.xray_scattering_factor_resonant(atom_type, qmag, enval)
                    elif self._use_waaskirf_scattering_factor:
                        ff = fc.xray_scattering_factor_WaasKirf(atom_type, qmag)
                    else:
                        ff = fc.xray_scattering_factor(atom_type, qmag)

                    # Get Debye-Waller factor
                    if self._use_isotropic_thermal_factor:
                        dw = fc.debyewaller(uiso, qmag)
                    elif self._use_anisotropic_thermal_factor:
                        raise Exception('anisotropic thermal factor calcualtion not implemented yet')
                    else:
                        dw = None

                    # Get magnetic form factors
                    if self._use_magnetic_form_factor:
                        mf = fc.magnetic_form_factor(atom_type, qmag)
                    else:
                        mf = None

                    options['scattering_factor'] = ff
                    options['debyewaller'] = dw
                    options['magnetic_formfactor'] = mf
                    options['energy_kev'] = enval
                    options['psi'] = psival

                    sf[ls: ls + len(_q), e, p] = scattering_fun(_q, r, **options)
                    ls = ls + len(_q)

        end_time = datetime.datetime.now()
        time_difference = end_time - start_time
        if self._time_report and time_difference.total_seconds() > 10:
            print('Calculated %d structure factors in %s' % (nref, time_difference))
        sf = sf / self.xtl.scale
        if nenergy == 1 and npsi == 1:
            return sf[:, 0, 0]  # shape(nref)
        if nenergy == 1:
            return sf[:, 0, :]  # shape(nref, nenergy)
        if npsi == 1:
            return sf[:, :, 0]  # shape(nref, nspi)
        return sf
    new_structure_factor = structure_factor

    def intensity(self, hkl=None, scattering_type=None, **options):
        """

        :param hkl:
        :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']
        :return:
        """
        return fs.intensity(self.new_structure_factor(hkl, scattering_type, **options))
    new_intensity = intensity

    def powder(self, scattering_type=None, units=None, peak_width=None, background=None,
               pixels=None, powder_average=None, **options):
        """

        Generates array of intensities along a spaced grid, equivalent to a powder pattern.
          tth, iten = generate_powder('xray', units='tth', energy_kev=8)

        :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']
        :param units: str : one of ['tth', 'dspace', 'q']
        :param peak_width: float : Peak with in units of inverse wavevector (Q)
        :param background: float : if >0, a normal background around this value will be added
        :param pixels: int : number of pixels to add to the resulting mesh
        :param powder_average: Bool : if True, intensities will be reduced for the powder average
        :param options: additional arguments to pass to intensity calculation
        :return xval: arrray : x-axis of powder scan (units)
        :return yval: array :  intensity values at each point in x-axis
        """
        if scattering_type is None:
            scattering_type = self._scattering_type
        if units is None:
            units = self._powder_units
        if peak_width is None:
            peak_width = self._powder_peak_width
        if background is None:
            background = self._powder_background
        if pixels is None:
            pixels = self._powder_pixels
        if powder_average is None:
            powder_average = self._powder_average
        if 'energy_kev' in options:
            energy_kev = options['energy_kev']
        else:
            energy_kev = self._energy_kev

        # Units
        min_twotheta = self._scattering_min_twotheta
        if min_twotheta <= 0: min_twotheta = 1.0
        max_twotheta = self._scattering_max_twotheta
        q_min = fc.calqmag(min_twotheta, energy_kev)
        q_max = fc.calqmag(max_twotheta, energy_kev)

        # Get reflections
        hmax, kmax, lmax = fc.maxHKL(q_max, self.xtl.Cell.UVstar())
        HKL = fc.genHKL([hmax, -hmax], [kmax, -kmax], [lmax, -lmax])
        HKL = self.xtl.Cell.sort_hkl(HKL)  # required for labels
        Qmag = self.xtl.Cell.Qmag(HKL)
        HKL = HKL[Qmag < q_max, :]
        Qmag = self.xtl.Cell.Qmag(HKL)

        # Calculate intensities
        I = self.intensity(HKL, scattering_type, **options)

        if powder_average:
            # Apply powder averging correction, I0/|Q|**2
            I = I/(Qmag+0.001)**2

        # create plotting mesh
        pixels = int(pixels * q_max)  # reduce this to make convolution faster
        pixel_size = q_max / (1.0 * pixels)
        peak_width_pixels = peak_width / (1.0 * pixel_size)
        mesh = np.zeros([pixels])
        mesh_q = np.linspace(q_min, q_max, pixels)

        # add reflections to background
        pixel_coord = Qmag / (1.0 * q_max)
        pixel_coord = (pixel_coord * (pixels - 1)).astype(int)

        for n in range(1, len(I)):
            mesh[pixel_coord[n]] = mesh[pixel_coord[n]] + I[n]

        # Convolve with a gaussian (if >0 or not None)
        if peak_width:
            gauss_x = np.arange(-3*peak_width_pixels, 3*peak_width_pixels + 1)  # gaussian width = 2*FWHM
            G = fg.gauss(gauss_x, None, height=1, cen=0, fwhm=peak_width_pixels, bkg=0)
            mesh = np.convolve(mesh, G, mode='same')

        # Add background (if >0 or not None)
        if background:
            bkg = np.random.normal(background, np.sqrt(background), [pixels])
            mesh = mesh + bkg

        # Change output units
        if units.lower() in ['tth', 'angle', 'two-theta', 'twotheta', 'theta']:
            xval = fc.cal2theta(mesh_q, energy_kev)
        elif units.lower() in ['d', 'dspace', 'd-spacing', 'dspacing']:
            xval = fc.q2dspace(mesh_q)
        else:
            xval = mesh_q
        return xval, mesh

    def x_ray(self, HKL):
        """
        Calculate the squared structure factor for the given HKL, using x-ray scattering factors
          Scattering.x_ray([1,0,0])
          Scattering.x_ray([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1, 3])
        Nref = len(HKL)
        
        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        Nat = len(uvw)
        
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get atomic form factors
        if self._use_waaskirf_scattering_factor:
            ff = fc.xray_scattering_factor_WaasKirf(atom_type, Qmag)
        else:
            ff = fc.xray_scattering_factor(atom_type, Qmag)
        
        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso, Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calcualtion not implemented yet')
        else:
            dw = 1
        
        # Calculate dot product
        dot_KR = np.dot(HKL, uvw.T)
        
        # Calculate structure factor
        # Broadcasting used on 2D ff
        SF =  np.sum(ff*dw*occ*np.exp(1j*2*np.pi*dot_KR),axis=1)
        #SF = np.zeros(Nref,dtype=np.complex)
        #for ref in range(Nref):
        #    for at in range(Nat): 
        #        SF[ref] += ff[ref,at]*dw[ref,at]*occ[at]*np.exp(1j*2*np.pi*dot_KR[ref,at])
        
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def x_ray_fast(self, HKL):
        """
        Calculate the squared structure factor for the given HKL, using atomic number as scattering length
          Scattering.x_ray_fast([1,0,0])
          Scattering.x_ray_fast([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        
        uvw,type,label,occ,uiso,mxmymz = self.xtl.Structure.get()
        
        # Get atomic form factors
        ff = fc.atom_properties(type, 'Z')
        
        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            Qmag = self.xtl.Cell.Qmag(HKL)
            dw = fc.debyewaller(uiso,Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        # Calculate structure factor
        SF = np.sum(ff * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)
        
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)

    def xray_dispersion(self, HKL, energy_kev):
        """
        Calculate the squared structure factor for the given HKL,
        using x-ray scattering factors with dispersion corrections
          Scattering.xray_dispersion([1,0,0], 9.71)
          Scattering.xray_dispersion([[1,0,0],[2,0,0],[3,0,0]], 2.838)
          Scattering.xray_dispersion([[1,0,0],[2,0,0],[3,0,0]], np.arange(2.83, 2.86, 0.001))
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        if energy_kev is an array, the returned array will have shape [len(HKL), len(energy_kev)]
        """

        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])

        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()

        Qmag = self.xtl.Cell.Qmag(HKL)

        # Get atomic form factors
        ff = fc.xray_scattering_factor_resonant(atom_type, Qmag, energy_kev)  # shape (len(HKL), len(type), len(en))

        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso, Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calcualtion not implemented yet')
        else:
            dw = 1

        # Calculate dot product
        dot_KR = np.dot(HKL, uvw.T)

        energy_kev = np.asarray(energy_kev).reshape(-1)
        sf = np.zeros([len(HKL), len(energy_kev)], dtype=np.complex)
        for n, en in enumerate(energy_kev):
            # Broadcasting used on 2D ff
            sf[:, n] = np.sum(ff[:, :, n] * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)
        if len(energy_kev) == 1:
            sf = sf[:, 0]
        elif len(HKL) == 1:
            sf = sf[0, :]

        sf = sf / self.xtl.scale

        if self._return_structure_factor: return sf

        # Calculate intensity
        I = sf * np.conj(sf)
        return np.real(I)
    
    def neutron(self, HKL):
        """
        Calculate the squared structure factor for the given HKL, using neutron scattering length
          Scattering.neutron([1,0,0])
          Scattering.neutron([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        """
        
        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])

        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        
        # Get atomic form factors
        ff = fc.atom_properties(atom_type, 'Coh_b')
        
        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            Qmag = self.xtl.Cell.Qmag(HKL)
            dw = fc.debyewaller(uiso,Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        # Calculate structure factor
        SF = np.sum(ff * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)
        
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def magnetic_neutron(self, HKL):
        """
        Calculate the magnetic component of the structure factor for the given HKL, using neutron rules and form factor
          Scattering.magnetic_neutron([1,0,0])
          Scattering.magnetic_neutron([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)

        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        Nat = len(uvw)
        
        Q = self.xtl.Cell.calculateQ(HKL)
        Qhat = fg.norm(Q).reshape([-1,3])
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get magnetic form factors
        if self._use_magnetic_form_factor:
            ff = fc.magnetic_form_factor(atom_type, Qmag)
        else:
            ff = np.ones([len(HKL), Nat])
        
        # Calculate moment
        momentmag = fg.mag(mxmymz).reshape([-1,1])
        momentxyz = self.xtl.Cell.calculateR(mxmymz)
        moment = momentmag*fg.norm(momentxyz) # broadcast n*1 x n*3 = n*3
        moment[np.isnan(moment)] = 0.
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)

        # Calculate structure factor
        SF = np.zeros(Nref,dtype=np.complex)
        for n,Qh in enumerate(Qhat):
            SFm = [0.,0.,0.]
            for m,mom in enumerate(moment):
                # Calculate Magnetic part
                QM = mom - np.dot(Qh,mom)*Qh

                # Calculate structure factor
                SFm = SFm + ff[n, m] * np.exp(1j * 2 * np.pi * dot_KR[n, m]) * QM
            
            # Calculate polarisation with incident neutron
            if self._polarised:
                SF[n] = np.dot(SFm, self._polarisation_vector_incident)
            else:
                #SF[n] = np.dot(SFm,SFm) # maximum possible
                SF[n] = (np.dot(SFm,[1,0,0]) + np.dot(SFm,[0,1,0]) + np.dot(SFm,[0,0,1]))/3 # average polarisation

        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def xray_magnetic(self, HKL):
        """
        Calculate the non-resonant magnetic component of the structure factor 
        for the given HKL, using x-ray rules and form factor
          Scattering.xray_magnetic([1,0,0])
          Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        
        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)
        Book: "X-ray Scattering and Absorption by Magnetic Materials" by <NAME> Collins. Ch 2. Eqn.2.21+1
        No orbital component assumed
        magnetic moments assumed to be in the same reference frame as the polarisation
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)
        
        uvw,type,label,occ,uiso,mxmymz = self.xtl.Structure.get()
        Nat = len(uvw)
        
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get magnetic form factors
        if self._use_magnetic_form_factor:
            ff = fc.magnetic_form_factor(type,Qmag)
        else:
            ff = np.ones([len(HKL),Nat])
        
        # Calculate moment
        momentmag = fg.mag(mxmymz).reshape([-1,1])
        momentxyz = self.xtl.Cell.calculateR(mxmymz) # moment direction in cartesian reference frame
        moment = momentmag*fg.norm(momentxyz) # broadcast n*1 x n*3 = n*3
        moment[np.isnan(moment)] = 0.
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        # Calculate structure factor
        SF = np.zeros(Nref,dtype=np.complex)
        for n in range(Nref):
            # Calculate vector structure factor
            SFm = [0.,0.,0.]
            for m,mom in enumerate(moment):
                SFm = SFm + ff[n,m]*np.exp(1j*2*np.pi*dot_KR[n,m])*mom
            
            # Calculate polarisation with incident x-ray
            # The reference frame of the x-ray and the crystal are assumed to be the same
            # i.e. pol=[1,0,0] || mom=[1,0,0] || (1,0,0)
            if self._polarised:
                SF[n] = np.dot(SFm,self._polarisation_vector_incident)
            else:
                #SF[n] = np.dot(SFm,SFm) # maximum possible
                SF[n] = (np.dot(SFm,[1,0,0]) + np.dot(SFm,[0,1,0]) + np.dot(SFm,[0,0,1]))/3 # average polarisation
        
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def xray_resonant(self, HKL, energy_kev=None, polarisation='sp', F0=1, F1=1, F2=1, azim_zero=[1,0,0], PSI=[0], disp=False):
        """
        Calculate structure factors using resonant scattering factors in the dipolar approximation
          I = Scattering.xray_resonant(HKL,energy_kev,polarisation,F0,F1,F2)
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
            energy_kev = x-ray energy in keV
            polarisation = x-ray polarisation: 'ss',{'sp'},'ps','pp'
            F0/1/2 = Resonance factor Flm
            azim_zero = [h,k,l] vector parallel to the origin of the azimuth
            psi = azimuthal angle defining the scattering plane
        
        Uses the E1E1 resonant x-ray scattering amplitude:
            fxr_n = (ef.ei)*F0 -i(ef X ei).z_n*F1 + (ef.z_n)(ei.z_n)F2
        
        Where ei and ef are the initial and final polarisation states, respectively,
        and z_n is a unit vector in the direction of the magnetic moment of the nth ion.
        The polarisation states are determined to be one of the natural synchrotron 
        states, where sigma (s) is perpendicular to the scattering plane and pi (p) is
        parallel to it.
                ( s-s  s-p )
                ( p-s  p-p )
        
        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (15)
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)
        
        PSI = np.asarray(PSI,dtype=np.float).reshape([-1])
        Npsi = len(PSI)

        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso,Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        SF = np.zeros([Nref,Npsi],dtype=np.complex)
        for psival in range(Npsi):
            # Get resonant form factor
            fxres = self.xray_resonant_scattering_factor(HKL,energy_kev,polarisation,F0,F1,F2,azim_zero,PSI[psival],disp=disp)
            
            # Calculate structure factor
            # Broadcasting used on 2D fxres
            SF[:,psival] = np.sum(fxres*dw*occ*np.exp(1j*2*np.pi*dot_KR),axis=1)
            
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def xray_resonant_scattering_factor(self, HKL, energy_kev=None, polarisation='sp', F0=1, F1=1, F2=1, azim_zero=[1,0,0], psi=0, disp=False):
        """
        Calcualte fxres, the resonant x-ray scattering factor
          fxres = Scattering.xray_resonant_scattering_factor(HKL,energy_kev,polarisation,F0,F1,F2,azim_zero,psi)
        energy_kev = x-ray energy in keV
            polarisation = x-ray polarisation: 'ss',{'sp'},'ps','pp'
            F0/1/2 = Resonance factor Flm
            azim_zero = [h,k,l] vector parallel to the origin of the azimuth {[1,0,0]}
            psi = azimuthal angle defining the scattering plane {0}
        
        Uses the E1E1 resonant x-ray scattering amplitude:
            fxr_n = (ef.ei)*F0 -i(ef X ei).z_n*F1 + (ef.z_n)(ei.z_n)F2
        
        Where ei and ef are the initial and final polarisation states, respectively,
        and z_n is a unit vector in the direction of the magnetic moment of the nth ion.
        The polarisation states are determined to be one of the natural synchrotron 
        states, where sigma (s) is perpendicular to the scattering plane and pi (p) is
        parallel to it.
                ( s-s  s-p )
                ( p-s  p-p )
        
        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (15)
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)
        Qmag = self.xtl.Cell.Qmag(HKL)
        tth = fc.cal2theta(Qmag, energy_kev)
        
        mxmymz = self.xtl.Structure.mxmymz()
        Nat = len(mxmymz)
        
        fxres = np.zeros([Nref,Nat],dtype=np.complex)
        for ref in range(Nref):
            # Resonant scattering factor
            # Electric Dipole transition at 3d L edge
            #F0,F1,F2 = 1,1,1 # Flm (form factor?)
            z1,z2,z3 = self.scatteringcomponents(mxmymz, HKL[ref], azim_zero, psi).T
            tthr = np.deg2rad(tth[ref])/2.0

            polarisation = polarisation.replace('-', '').replace(' ', '')
            if polarisation in ['sigmasigma', 'sigsig', 'ss']:    # Sigma-Sigma
                f0 = 1*np.ones(Nat)
                f1 = 0*np.ones(Nat)
                f2 = z2**2
            elif polarisation in ['sigmapi', 'sigpi', 'sp']:  # Sigma-Pi
                f0 = 0 * np.ones(Nat)
                f1 = z3 * np.sin(tthr) + z1 * np.cos(tthr)
                f2 = z2 * (z1 * np.sin(tthr) + z3 * np.cos(tthr))
            elif polarisation in ['pisigma', 'pisig', 'ps']:  # Pi-Sigma
                f0 = 0*np.ones(Nat)
                f1 = z1*np.cos(tthr) - z3*np.sin(tthr)
                f2 = -z2*(z1*np.sin(tthr)-z3*np.cos(tthr))
            elif polarisation in ['pipi', 'pp']:  # Pi-Pi
                f0 = np.cos(2*tthr)*np.ones(Nat)
                f1 = -z2*np.sin(2*tthr)
                f2 = -(np.cos(tthr)**2)*(z1**2*np.tan(tthr)**2 + z3**2)
            else:
                raise ValueError('Incorrect polarisation. pol should be e.g. ''ss'' or ''sp''')
            fxres[ref,:] = F0*f0 -1j*F1*f1 + F2*f2
            if disp:
                print('( h, k, l)   TTH  (    mx,    my,    mz)  (    z1,    z2,    z3)')
                fmt = '(%2.0f,%2.0f,%2.0f) %6.2f  (%6.3f,%6.3f,%6.3f)  (%6.3f,%6.3f,%6.3f)  f0=%8.4f  f1=%8.4f  f2=%8.4f fxres= (%8.4f + %8.4fi)'
                for at in range(Nat):
                    vals = (HKL[ref,0],HKL[ref,1],HKL[ref,2],tth[ref],
                            mxmymz[at,0],mxmymz[at,1],mxmymz[at,2],
                            z1[at],z2[at],z3[at],
                            f0[at],f1[at],f2[at],
                            fxres[ref,at].real,fxres[ref,at].imag)
                    print(fmt%vals)
        return fxres

    def xray_nonresonant_magnetic(self, HKL, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p', disp=False):
        """
        Calculate the non-resonant magnetic component of the structure factor
        for the given HKL, using x-ray rules and form factor
          Scattering.xray_magnetic([1,0,0])
          Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.

        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)
        Book: "X-ray Scattering and Absorption by Magnetic Materials" by Loevesy and Collins. Ch 2. Eqn.2.21+1
        No orbital component assumed
        magnetic moments assumed to be in the same reference frame as the polarisation
        """

        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])

        uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()

        kin, kout, ein, eout = self.scatteringvectors(HKL, energy_kev, azim_zero, psi, polarisation)

        Qmag = self.xtl.Cell.Qmag(HKL)

        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso, Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1

        # Get magnetic form factors
        if self._use_magnetic_form_factor:
            ff = fc.magnetic_form_factor(type, Qmag)
        else:
            ff = np.ones([len(HKL), len(uvw)])

        # Calculate moment
        momentmag = fg.mag(mxmymz).reshape([-1, 1])
        momentxyz = self.xtl.Cell.calculateR(mxmymz)  # moment direction in cartesian reference frame
        moment = momentmag * fg.norm(momentxyz)  # broadcast n*1 x n*3 = n*3
        moment[np.isnan(moment)] = 0.

        # Magnetic form factor
        # f_non-res_mag = i.r0.(hw/mc^2).fD.[.5.L.A + S.B] #equ 2 Hill+McMorrow 1996
        # ignore orbital moment L
        B = np.zeros([len(HKL), 3])
        for n in range(len(HKL)):
            #print(n,HKL[n],kin[n],kout[n],ein[n],eout[n])
            B[n, :] = np.cross(eout[n], ein[n]) + \
                np.cross(kout[n], eout[n]) * np.dot(kout[n], ein[n]) - \
                np.cross(kin[n], ein[n]) * np.dot(kin[n], eout[n]) - \
                np.cross(np.cross(kout[n], eout[n]), np.cross(kin[n], ein[n]))
        fspin = 1j * ff * np.dot(moment, B.T).T

        # Calculate dot product
        dot_KR = np.dot(HKL, uvw.T)

        # Calculate structure factor
        SF = np.sum(fspin * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)

        if disp:
            for n in range(len(HKL)):
                pr_h = '(%2.0f,%2.0f,%2.0f)'%tuple(HKL[n,:])
                pr_b = '(%6.2g,%6.2g,%6.2g)'%tuple(B[n,:])
                print('psi=%3.0f  %3d hkl=%10s B=%22s'%(psi, n, pr_h, pr_b))
                ctot = 0j
                for m in range(len(moment)):
                    if np.sum(moment[m,:]**2) < 0.01: continue
                    pr_m = '(%6.2g,%6.2g,%6.2g)'%tuple(moment[m,:])
                    dot = np.dot(moment[m],B[n,:])
                    pdot = '%6.3f'%dot
                    phase = np.exp(1j * 2 * np.pi * dot_KR[n, m])
                    prph = '%5.2f+i%5.2f'%(np.real(phase),np.imag(phase))
                    prod = dot*phase
                    pprd = '%5.2f+i%5.2f'%(np.real(prod),np.imag(prod))
                    ctot += prod
                    ptot = '%5.2f+i%5.2f'%(np.real(ctot),np.imag(ctot))
                    print('\t%3d mom=%22s dot(mom,B)=%6s   exp(ik.r)=%12s   sum=%12s   tot=%12s'%(m,pr_m,pdot,prph,pprd,ptot))

        SF = SF / self.xtl.scale

        if self._return_structure_factor: return SF

        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)

    def xray_resonant_magnetic(self, HKL, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p', F0=0, F1=1, F2=0, disp=True):
        """
        Calculate the non-resonant magnetic component of the structure factor
        for the given HKL, using x-ray rules and form factor
          Scattering.xray_magnetic([1,0,0])
          Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.

        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)
        Book: "X-ray Scattering and Absorption by Magnetic Materials" by <NAME> Collins. Ch 2. Eqn.2.21+1
        No orbital component assumed
        magnetic moments assumed to be in the same reference frame as the polarisation
        """

        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])

        uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()

        kin, kout, ein, eout = self.scatteringvectors(HKL, energy_kev, azim_zero, psi, polarisation)

        Qmag = self.xtl.Cell.Qmag(HKL)

        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso, Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1

        # Calculate moment
        momentmag = fg.mag(mxmymz).reshape([-1, 1])
        momentxyz = self.xtl.Cell.calculateR(mxmymz)  # moment direction in cartesian reference frame
        moment = momentmag * fg.norm(momentxyz)  # broadcast n*1 x n*3 = n*3
        moment[np.isnan(moment)] = 0.
        fe1e1 = np.zeros([len(HKL), len(uvw)], dtype=np.complex)
        for ref in range(len(HKL)):
            # Magnetic form factor
            # f_res_mag = [(e'.e)F0 - i(e'xe).Z*F1 + (e'.Z)*(e.Z)*F2] #equ 7 Hill+McMorrow 1996
            f0 = np.dot(eout[ref], ein[ref])
            f1 = np.dot(np.cross(eout[ref], ein[ref]), moment.T)
            f2 = np.dot(eout[ref], moment.T) * np.dot(ein[ref], moment.T)
            fe1e1[ref, :] = f0*F0 - 1j*f1*F1 + f2*F2

        # Calculate dot product
        dot_KR = np.dot(HKL, uvw.T)

        if disp:
            for ref in range(len(HKL)):
                pr_h = '(%2.0f,%2.0f,%2.0f)'%tuple(HKL[ref,:])
                print('psi=%3.0f  %3d hkl=%10s'%(psi, ref, pr_h))
                ctot = 0j
                for m in range(len(moment)):
                    if np.sum(moment[m,:]**2) < 0.01: continue
                    pr_m = '(%6.2g,%6.2g,%6.2g)'%tuple(moment[m,:])
                    f0 = np.dot(eout[ref], ein[ref])
                    f1 = np.dot(np.cross(eout[ref], ein[ref]), moment[m,:])
                    f2 = np.dot(eout[ref], moment[m,:]) * np.dot(ein[ref], moment[m,:])
                    dot = f0*F0 - 1j*f1*F1 + f2*F2
                    phase = np.exp(1j * 2 * np.pi * dot_KR[ref, m])
                    prph = '%5.2f+i%5.2f'%(np.real(phase),np.imag(phase))
                    prod = dot*phase
                    pprd = '%5.2f+i%5.2f'%(np.real(prod),np.imag(prod))
                    ctot += prod
                    ptot = '%5.2f+i%5.2f'%(np.real(ctot),np.imag(ctot))
                    print('\t%3d mom=%22s f0=%5.2f f1=%5.2f f2=%5.2f   exp(ik.r)=%12s   sum=%12s   tot=%12s'%(m,pr_m,f0,f1,f2,prph,pprd,ptot))

        # Calculate structure factor
        SF = np.sum(fe1e1 * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)

        SF = SF / self.xtl.scale

        if self._return_structure_factor: return SF

        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)

    def scatteringvectors(self, hkl, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p'):
        """
        Determine the scattering and polarisation vectors of a reflection based on energy, azimuth and polarisation.
        :param xtl: Crystal Class
        :param hkl: [n,3] array of reflections
        :param energy_kev: x-ray scattering energy in keV
        :param azim_zero: [1,3] direction along which the azimuthal zero angle is determind
        :param psi: float angle in degrees about the azimuth
        :param polarisation: polarisation with respect to the scattering plane, options:
                    'ss' : sigma-sigma polarisation
                    'sp' : sigma-pi polarisation
                    'ps' : pi-sigma polarisation
                    'pp' : pi-pi polarisation
                or: polarisation: float polarisation angle of scattered vector in degrees
        :return: kin, kout, ein, eout
        Returned values are [n,3] arrays
            kin : [n,3] array of incident wavevectors
            kout: [n,3] array of scattered wavevectors
            ein : [n,3] array of incident polarisation
            eout: [n,3] array of scattered polarisation

        The basis is chosen such that Q defines the scattering plane, sigma and pi directions are normal to this plane.
        Q is defined as Q = kout - kin, with kout +ve along the projection of azim_zero
        """

        hkl = np.asarray(hkl).reshape([-1, 3])

        if energy_kev is None:
            energy_kev = self._energy_kev

        out_kin = np.zeros([len(hkl), 3])
        out_kout = np.zeros([len(hkl), 3])
        out_ein = np.zeros([len(hkl), 3])
        out_eout = np.zeros([len(hkl), 3])
        for n in range(len(hkl)):
            # Define coordinate system I,J,Q (U1,U2,U3)
            Ihat, Jhat, Qhat = self.scatteringbasis(hkl[n, :], azim_zero, psi)
            Ihat = Ihat.reshape([-1, 3])
            Jhat = Jhat.reshape([-1, 3])
            Qhat = Qhat.reshape([-1, 3])

            # Determine wavevectors
            bragg = self.xtl.Cell.tth(hkl[n, :], energy_kev) / 2.
            rb = np.deg2rad(bragg).reshape([-1, 1])
            kin = np.cos(rb) * Ihat - np.sin(rb) * Qhat
            kout = np.cos(rb) * Ihat + np.sin(rb) * Qhat
            esig = Jhat  # sigma polarisation (in or out)
            piin = np.cross(kin, esig)  # pi polarisation in
            piout = np.cross(kout, esig)  # pi polarisation out

            # Polarisation
            try:
                # polarisation = 'ss' or 's-s'
                polarisation = polarisation.replace('-', '').replace(' ', '')
            except AttributeError:
                # polarisation = angle in deg from sigma' to pi'
                ein = 1.0 * esig
                pol = np.deg2rad(polarisation)
                eout = np.cos(pol)*esig + np.sin(pol)*piout
            if polarisation in ['sigmasigma', 'sigsig', 'ss']:
                ein = 1.0 * esig
                eout = 1.0 * esig
            elif polarisation in ['sigmapi', 'sigpi', 'sp']:
                ein = 1.0 * esig
                eout = 1.0 * piout
            elif polarisation in ['pisigma', 'pisig', 'ps']:
                ein = 1.0 * piin
                eout = 1.0 * esig
            elif polarisation in ['pipi', 'pp']:
                ein = 1.0 * piin
                eout = 1.0 * piout
            out_kin[n, :] = kin
            out_kout[n, :] = kout
            out_ein[n, :] = ein
            out_eout[n, :] = eout

        return out_kin, out_kout, out_ein, out_eout

    def scatteringcomponents(self, mxmymz, hkl, azim_zero=[1,0,0], psi=0):
        """
        Transform magnetic vector into components within the scattering plane
            ***warning - may not be correct for non-cubic systems***
        """

        # Define coordinate system I,J,Q (U1,U2,U3)
        U = self.scatteringbasis(hkl, azim_zero, psi)
        
        # Determine components of the magnetic vector
        z1z2z3 = np.dot(mxmymz, U.T) # [mxmymz.I, mxmymz.J, mxmymz.Q]
        return fg.norm(z1z2z3)

    def scatteringbasis(self, hkl, azim_zero=[1, 0, 0], psi=0):
        """
        Determine the scattering and polarisation vectors of a reflection based on energy, azimuth and polarisation.
        :param hkl: [n,3] array of reflections
        :param azim_zero: [1,3] direction along which the azimuthal zero angle is determind
        :param psi: float azimuthal angle about U3 in degrees
        :return: U1, U2, U3
        The basis is chosen such that Q defines the scattering plane, the sigma direction is normal to this plane,
        the pi direction is always within this plane.
        The azimuthal angle defines a rotation about the Q axis in a clockwise mannor, matching I16.
        At an azimuth of 0degrees, U1 is perpendicular to Q, along the direction of azim_zero.
        """

        # Define coordinate system I,J,Q (U1,U2,U3)
        # See FDMNES User's Guide p20 'II-11) Anomalous or resonant diffraction'
        # U1 || projection of azim_zero
        # U2 _|_ U1,U3
        # U3 || Q = kf-ki
        azim_zero = fg.norm(self.xtl.Cell.calculateQ(azim_zero)) # put in orthogonal basis
        Qhat = fg.norm(self.xtl.Cell.calculateQ(hkl)).reshape([-1,3])  # || Q
        AxQ = fg.norm(np.cross(azim_zero, Qhat))
        Ihat = fg.norm(np.cross(Qhat, AxQ)).reshape([-1,3])  # || to projection of azim_zero
        Jhat = fg.norm(np.cross(Qhat, Ihat)).reshape([-1,3])  # _|_ to I and Q

        # Rotate psi about Qhat
        rpsi = np.deg2rad(psi)
        # -ve sin makes clockwise rotation
        # This was checked on 21/1/19 vs CRO paper + sergio's calculations and seems to agree with experiment,
        # however we never did an azimuthal scan of the (103) which would have distinguished this completely.
        Ihat_psi = fg.norm(np.cos(rpsi) * Ihat - np.sin(rpsi) * Jhat)
        Jhat_psi = fg.norm(np.cross(Qhat, Ihat_psi))
        return np.vstack([Ihat_psi, Jhat_psi, Qhat])

    def print_scattering_coordinates(self,hkl,azim_zero=[1,0,0],psi=0):
        """
        Transform magnetic vector into components within the scattering plane
            ***warning - may not be correct for non-cubic systems***
        """
        
        # Define coordinate system I,J,Q (U1,U2,U3)
        Qhat = fg.norm(self.xtl.Cell.calculateQ(hkl)) # || Q
        AxQ = fg.norm(np.cross(azim_zero,Qhat))
        Ihat = fg.norm(np.cross(Qhat,AxQ)) # || to azim_zero
        Jhat = fg.norm(np.cross(Qhat,Ihat)) # -| to I and Q
        
        # Rotate coordinate system by azimuth
        Ihat_psi = fg.norm(np.cos(np.deg2rad(psi))*Ihat + np.sin(np.deg2rad(psi))*Jhat)
        Jhat_psi = fg.norm(np.cross(Qhat,Ihat_psi))
        
        # Determine components of the magnetic vector
        U=np.vstack([Ihat_psi,Jhat_psi,Qhat])
        print('U1 = (%5.2f,%5.2f,%5.2f)'%(U[0,0],U[0,1],U[0,2]))
        print('U2 = (%5.2f,%5.2f,%5.2f)'%(U[1,0],U[1,1],U[1,2]))
        print('U3 = (%5.2f,%5.2f,%5.2f)'%(U[2,0],U[2,1],U[2,2]))
    
    def print_intensity(self, HKL):
        """
        Print intensities calcualted in different ways
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Qmag =  self.xtl.Cell.Qmag(HKL)
        srt = np.argsort(Qmag)
        HKL = HKL[srt,:]
        
        IN=self.neutron(HKL)
        IX=self.x_ray(HKL)
        INM=self.magnetic_neutron(HKL)*1e4
        IXM=self.xray_magnetic(HKL)*1e4
        IXRss=self.xray_resonant(HKL, None, 'ss')
        IXRsp=self.xray_resonant(HKL, None, 'sp')
        IXRps=self.xray_resonant(HKL, None, 'ps')
        IXRpp=self.xray_resonant(HKL, None, 'pp')
        
        fmt = '(%2.0f,%2.0f,%2.0f)  %8.1f  %8.1f  %8.2f  %8.2f  ss=%8.2f  sp=%8.2f  ps=%8.2f  pp=%8.2f'
        print('( h, k, l)   Neutron      xray   Magn. N  Magn. XR   sig-sig    sig-pi    pi-sig     pi-pi')
        for n in range(len(HKL)):
            vals=(HKL[n][0],HKL[n][1],HKL[n][2],IN[n],IX[n],INM[n],IXM[n],IXRss[n],IXRsp[n],IXRps[n],IXRpp[n])
            print(fmt%vals)
    
    def old_intensity(self, HKL, scattering_type=None):
        """
        Calculate the squared structure factor for the given HKL
          Crystal.intensity([1,0,0])
          Crystal.intensity([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        
        Notes:
        - Uses x-ray atomic form factors, calculated from approximated tables in the ITC
        - This may be a little slow for large numbers of reflections, as it is not currently
         possible to use accelerated calculation methods in Jython.
        - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs
        - Crystal.scale is used to scale the complex structure factor, so the intensity is  
         reduced by (Crystal.scale)^2
        - Testing against structure factors calculated by Vesta.exe is very close, though there
          are some discrepancies, probably due to the method of calculation of the form factor.
        """

        if scattering_type is None:
            scattering_type = self._scattering_type
        scattering_type = scattering_type.lower()

        # Break up long lists of HKLs
        n_arrays = np.ceil(len(HKL)*len(self.xtl.Structure.u)/10000.)
        hkl_array = np.array_split(HKL, n_arrays)

        intensity = []
        for _hkl in hkl_array:
            if scattering_type in ['xray','x','x-ray','thomson','charge']:
                intensity += self.x_ray(_hkl).tolist()
            elif scattering_type in ['neutron','n','nuclear']:
                intensity += self.neutron(_hkl).tolist()
            elif scattering_type in ['xray magnetic','magnetic xray','spin xray','xray spin']:
                intensity += list(self.xray_magnetic(_hkl)*1e4)
            elif scattering_type in ['neutron magnetic','magnetic neutron','magnetic']:
                intensity += list(self.magnetic_neutron(_hkl)*1e4)
            elif scattering_type in ['xray dispersion']:
                intensity += self.xray_dispersion(_hkl, self._energy_kev).tolist()
            elif scattering_type in ['xray resonant','resonant','resonant xray','rxs']:
                intensity += self.xray_resonant(_hkl).tolist()
            elif scattering_type in ['xray resonant magnetic', 'xray magnetic resonant',
                                     'resonant magnetic', 'magnetic resonant']:
                intensity += self.xray_resonant_magnetic(
                    _hkl,
                    self._energy_kev,
                    self._azimuthal_reference,
                    self._azimuthal_angle,
                    self._polarisation,
                    F0=0, F1=1, F2=0).tolist()
            elif scattering_type in ['xray nonresonant magnetic', 'xray magnetic nonresonant',
                                     'nonresonant magnetic', 'magnetic nonresonant',
                                     'xray non-resonant magnetic', 'xray magnetic non-resonant',
                                     'non-resonant magnetic', 'magnetic non-resonant']:
                intensity += self.xray_resonant_magnetic(
                    _hkl,
                    self._energy_kev,
                    self._azimuthal_reference,
                    self._azimuthal_angle,
                    self._polarisation).tolist()
            else:
                print('Scattering type not defined')
        return np.array(intensity)

    def old_structure_factor(self, HKL, scattering_type=None):
        """
        Calculate the complex structure factor for the given HKL
          Crystal.structure_factor([1,0,0])
          Crystal.structure_factor([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the complex structure factor at each reflection.

        Notes:
        - Uses x-ray atomic form factors, calculated from approximated tables in the ITC
        - This may be a little slow for large numbers of reflections, as it is not currently
         possible to use accelerated calculation methods in Jython.
        - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs
        - Crystal.scale is used to scale the complex structure factor, so the intensity is
         reduced by (Crystal.scale)^2
        - Testing against structure factors calculated by Vesta.exe is very close, though there
          are some discrepancies, probably due to the method of calculation of the form factor.
        """
        prev_sf_setting = self._return_structure_factor
        self._return_structure_factor = True
        sf = self.intensity(HKL, scattering_type)
        self._return_structure_factor = prev_sf_setting
        return sf
    
    def hkl(self, HKL, energy_kev=None):
        """ Calculate the two-theta and intensity of the given HKL, display the result"""
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        tth = self.xtl.Cell.tth(HKL,energy_kev)
        inten = self.intensity(HKL)
        
        print('Energy = %6.3f keV' % energy_kev)
        print('( h, k, l) TwoTheta  Intensity')
        for n in range(len(tth)):
            print('(%2.0f,%2.0f,%2.0f) %8.2f  %9.2f' % (HKL[n,0],HKL[n,1],HKL[n,2],tth[n],inten[n]))
    
    def hkl_reflection(self, HKL, energy_kev=None):
        """
        Calculate the theta, two-theta and intensity of the given HKL in reflection geometry, display the result
        Uses sample orientation set up in setup_scatter
        :param HKL: [h,k,l] or list of hkl
        :param energy_kev: None or float
        :return: str
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        tth = self.xtl.Cell.tth(HKL,energy_kev)
        theta = self.xtl.Cell.theta_reflection(HKL, energy_kev, self._scattering_specular_direction, self._scattering_theta_offset)
        inten = self.intensity(HKL)
        
        print('Energy = %6.3f keV' % energy_kev)
        print('Specular Direction = (%1.0g,%1.0g,%1.0g)' %
              (self._scattering_specular_direction[0],
               self._scattering_specular_direction[1],
               self._scattering_specular_direction[2]))
        print('( h, k, l)    Theta TwoTheta  Intensity')
        for n in range(len(tth)):
            print('(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f' %
                  (HKL[n, 0], HKL[n, 1], HKL[n, 2], theta[n], tth[n], inten[n]))
    
    def hkl_transmission(self,HKL,energy_kev=None):
        " Calculate the theta, two-theta and intensity of the given HKL in transmission geometry, display the result"
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        tth = self.xtl.Cell.tth(HKL,energy_kev)
        theta = self.xtl.Cell.theta_transmission(HKL, energy_kev, self._scattering_specular_direction,self._scattering_theta_offset)
        inten = self.intensity(HKL)
        
        print('Energy = %6.3f keV' % energy_kev)
        print('Direction parallel to beam  = (%1.0g,%1.0g,%1.0g)' %(self._scattering_parallel_direction[0],self._scattering_parallel_direction[1],self._scattering_parallel_direction[2]))
        print('( h, k, l)    Theta TwoTheta  Intensity')
        for n in range(len(tth)):
            print('(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f' % (HKL[n,0],HKL[n,1],HKL[n,2],theta[n],tth[n],inten[n]))

    def generate_powder(self, q_max=8, peak_width=0.01, background=0, powder_average=True):
        """
        Generates array of intensities along a spaced grid, equivalent to a powder pattern.
          Q,I = generate_powder(energy_kev=8.0,peak_width=0.05,background=0)
            q_max = maximum Q, in A-1
            peak_width = width of convolution, in A-1
            background = average of normal background
            powder_average = True*/False, apply the powder averaging correction
          Returns:
            Q = [1000x1] array of wave-vector values
            I = [1000x1] array of intensity values

        Note: To get two-theta values use:
            tth = fc.cal2theta(Q, energy_kev)
        Note: To get d-spacing values use:
            dspace = fc.caldspace(Q)
        """

        # Get reflections
        hmax, kmax, lmax = fc.maxHKL(q_max, self.xtl.Cell.UVstar())
        HKL = fc.genHKL([hmax, -hmax], [kmax, -kmax], [lmax, -lmax])
        HKL = self.xtl.Cell.sort_hkl(HKL)  # required for labels
        Qmag = self.xtl.Cell.Qmag(HKL)
        HKL = HKL[Qmag < q_max, :]
        Qmag = self.xtl.Cell.Qmag(HKL)
        # Qmag = Qmag[Qmag<q_max]

        # Calculate intensities
        I = self.intensity(HKL)

        if powder_average:
            # Apply powder averging correction, I0/|Q|**2
            I = I/(Qmag+0.001)**2

        # create plotting mesh
        pixels = int(2000 * q_max)  # reduce this to make convolution faster
        pixel_size = q_max / (1.0 * pixels)
        peak_width_pixels = peak_width / (1.0 * pixel_size)
        mesh = np.zeros([pixels])
        mesh_q = np.linspace(0, q_max, pixels)

        # add reflections to background
        pixel_coord = Qmag / (1.0 * q_max)
        pixel_coord = (pixel_coord * (pixels - 1)).astype(int)

        for n in range(1, len(I)):
            mesh[pixel_coord[n]] = mesh[pixel_coord[n]] + I[n]

        # Convolve with a gaussian (if >0 or not None)
        if peak_width:
            gauss_x = np.arange(-3*peak_width_pixels, 3*peak_width_pixels + 1)  # gaussian width = 2*FWHM
            G = fg.gauss(gauss_x, None, height=1, cen=0, fwhm=peak_width_pixels, bkg=0)
            mesh = np.convolve(mesh, G, mode='same')

            # Add background (if >0 or not None)
        if background:
            bkg = np.random.normal(background, np.sqrt(background), [pixels])
            mesh = mesh + bkg
        return mesh_q, mesh

    def powder_correction(self, HKL, intensities, symmetric_multiplyer=True, powder_average=True):
        """
        Averages symmetric reflections and applies symmetry multipliers and 1/q^2 correction
        Ic = I0*C
        :param HKL: [nx3] array of [h,k,l] reflections
        :param intensities: [nx1] array of reflection intensities
        :return: [mx3], [mx1] arrays of averaged, corrected reflections + intensity
        """
        # Average symmetric reflections
        rhkl, rinten = self.xtl.Symmetry.average_symmetric_intensity(HKL, intensities)

        if symmetric_multiplyer:
            multiplyer = self.xtl.Symmetry.reflection_multiplyer(rhkl)
            rinten = rinten*multiplyer

        if powder_average:
            q = self.xtl.Cell.Qmag(rhkl)
            rinten = rinten/(q+0.001)**2

        return rhkl, rinten

    def print_all_reflections(self, energy_kev=None, print_symmetric=False,
                              min_intensity=0.01, max_intensity=None, units=None):
        """
        Prints a list of all allowed reflections at this energy
            energy_kev = energy in keV
            print_symmetric = False*/True : omits reflections with the same intensity at the same angle
            min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)
            max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)
            units = None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        if min_intensity is None: min_intensity = -1
        if max_intensity is None: max_intensity = np.inf
        
        hkl = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        if not print_symmetric:
            hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)
        hkl = self.xtl.Cell.sort_hkl(hkl)

        tth = self.xtl.Cell.tth(hkl, energy_kev)
        inrange = np.all([tth < self._scattering_max_twotheta, tth > self._scattering_min_twotheta], axis=0)
        hkl = hkl[inrange, :]
        tth = tth[inrange]
        #inten = np.sqrt(self.intensity(hkl)) # structure factor
        inten = self.intensity(hkl)

        if units is None:
            units = self._powder_units
        units = units.lower()
        if units in ['d', 'dspc', 'dspace', 'd space', 'd-space', 'dspacing', 'd spacing', 'd-spacing']:
            unit_str = 'd-spacing'
            unit = fc.caldspace(tth, energy_kev)
        elif units in ['q', 'wavevector']:
            unit_str = 'Q'
            unit = fc.calqmag(tth, energy_kev)
        else:
            unit_str = 'TwoTheta'
            unit = tth
        
        fmt = '(%3.0f,%3.0f,%3.0f) %10.2f  %9.2f\n'
        outstr = ''
        
        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Radiation: %s\n' % self._scattering_type
        outstr += '( h, k, l)    %10s  Intensity\n' % unit_str
        #outstr+= fmt % (hkl[0,0],hkl[0,1],hkl[0,2],tth[0],inten[0]) # hkl(0,0,0)
        count = 0
        for n in range(1, len(tth)):
            if inten[n] < min_intensity: continue
            if inten[n] > max_intensity: continue
            count += 1
            outstr += fmt % (hkl[n,0], hkl[n,1], hkl[n,2],unit[n],inten[n])
        outstr += 'Reflections: %1.0f\n' % count
        return outstr
    
    def print_ref_reflections(self, energy_kev=None, min_intensity=0.01, max_intensity=None):
        """
        Prints a list of all allowed reflections at this energy in the reflection geometry
            energy = energy in keV
            min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)
            max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)
                       |
                    // |\
                    \\ |
                     \\|___/____
                      \\   \
                       \\__
        
        Note, to change min/max theta values or specular direciton, change the following attributres of 
        the crystal object:
            self._scattering_max_two_theta   :  maximum detector (two-theta) angle
            self._scattering_min_theta       :  minimum sample angle = -opening angle
            self._scattering_max_theta       :  maximum sample angle = opening angle
            self._scattering_theta_offset    :  sample offset angle
            self._scattering_specular_direction : [h,k,l] : reflections normal to sample surface
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev

        if min_intensity is None: min_intensity = -1
        if max_intensity is None: max_intensity = np.inf

        HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        HKL = self.xtl.Cell.sort_hkl(HKL)
        tth = self.xtl.Cell.tth(HKL, energy_kev)
        HKL = HKL[tth > self._scattering_min_twotheta, :]
        tth = tth[tth > self._scattering_min_twotheta]
        theta = self.xtl.Cell.theta_reflection(HKL, energy_kev, self._scattering_specular_direction,
                                               self._scattering_theta_offset)
        # inten = np.sqrt(self.intensity(HKL)) # structure factor
        inten = self.intensity(HKL)

        p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)
        p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))
        pos_theta = p1 * p2
        
        fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'
        outstr = ''

        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Radiation: %s\n' % self._scattering_type
        outstr += 'Specular Direction = (%1.0g,%1.0g,%1.0g)\n' % (
            self._scattering_specular_direction[0],
            self._scattering_specular_direction[1],
            self._scattering_specular_direction[2]
        )
        outstr += '( h, k, l) TwoTheta    Theta  Intensity\n'
        # outstr+= fmt % (HKL[0,0],HKL[0,1],HKL[0,2],tth[0],theta[0],inten[0]) # hkl(0,0,0)
        count = 0
        for n in range(1, len(tth)):
            if inten[n] < min_intensity: continue
            if inten[n] > max_intensity: continue
            if not pos_theta[n]: continue
            # if not print_symmetric and np.abs(tth[n]-tth[n-1]) < 0.01: continue # only works if sorted
            count += 1
            outstr += fmt % (HKL[n, 0], HKL[n, 1], HKL[n, 2], tth[n], theta[n], inten[n])
        outstr += 'Reflections: %1.0f\n' % count
        return outstr
    
    def print_tran_reflections(self, energy_kev=None, min_intensity=0.01, max_intensity=None):
        """
        Prints a list of all allowed reflections at this energy in the transmission geometry
            energy = energy in keV
            min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)
            max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)
                   \ /      
             --<-- || --<-- 
                  / \       
        
        Note, to change min/max theta values or specular direciton, change the following attributres of 
        the crystal object:
            self._scattering_max_two_theta   :  maximum detector (two-theta) angle
            self._scattering_min_theta       :  minimum sample angle = -opening angle
            self._scattering_max_theta       :  maximum sample angle = opening angle
            self._scattering_theta_offset    :  sample offset angle
            self._scattering_parallel_direction : [h,k,l] : reflections normal to sample surface, parallel to beam at theta = 0
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev

        if min_intensity is None: min_intensity = -1
        if max_intensity is None: max_intensity = np.inf

        HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        HKL = self.xtl.Cell.sort_hkl(HKL)
        tth = self.xtl.Cell.tth(HKL, energy_kev)
        HKL = HKL[tth > self._scattering_min_twotheta, :]
        tth = tth[tth > self._scattering_min_twotheta]
        theta = self.xtl.Cell.theta_transmission(HKL, energy_kev, self._scattering_parallel_direction)
        # inten = np.sqrt(self.intensity(HKL)) # structure factor
        inten = self.intensity(HKL)

        p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)
        p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))
        pos_theta = p1 * p2

        fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'
        outstr = ''

        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Radiation: %s\n' % self._scattering_type
        outstr += 'Direction parallel to beam  = (%1.0g,%1.0g,%1.0g)\n' % (
            self._scattering_parallel_direction[0],
            self._scattering_parallel_direction[1],
            self._scattering_parallel_direction[2]
        )
        outstr += '( h, k, l) TwoTheta    Theta  Intensity\n'
        # outstr+= fmt % (HKL[0,0],HKL[0,1],HKL[0,2],tth[0],theta[0],inten[0]) # hkl(0,0,0)
        count = 0
        for n in range(1, len(tth)):
            if inten[n] < min_intensity: continue
            if inten[n] > max_intensity: continue
            if not pos_theta[n]: continue
            # if not print_symmetric and np.abs(tth[n]-tth[n-1]) < 0.01: continue # only works if sorted
            count += 1
            outstr += fmt % (HKL[n, 0], HKL[n, 1], HKL[n, 2], tth[n], theta[n], inten[n])
        outstr += ('Reflections: %1.0f\n' % count)
        return outstr
    
    def print_symmetric_reflections(self,HKL):
        "Prints equivalent reflections"
        
        symHKL = self.xtl.Symmetry.symmetric_reflections(HKL)
        Nsyms = len(symHKL)
        outstr = ''
        outstr+= 'Equivalent reflections: %d\n' % Nsyms
        for n in range(Nsyms):
            outstr+= '(%5.3g,%5.3g,%5.3g)\n' % (symHKL[n,0],symHKL[n,1],symHKL[n,2])
        return outstr
    
    def print_atomic_contributions(self,HKL):
        """
        Prints the atomic contributions to the structure factor
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)
        
        # Calculate the full intensity
        I = self.intensity(HKL)
        
        # Calculate the structure factors of the symmetric atomic sites
        base_label = self.xtl.Atoms.label
        uvw,type,label,occ,uiso,mxmymz = self.xtl.Structure.get()
        
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get atomic form factors
        ff = fc.xray_scattering_factor(type,Qmag)
        
        # Get Debye-Waller factor
        dw = fc.debyewaller(uiso,Qmag)
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        # Calculate structure factor
        SF =  ff*dw*occ*np.exp(1j*2*np.pi*dot_KR)
        
        # Sum structure factors of each base label in atoms
        SFbase = np.zeros([len(HKL),len(base_label)],dtype=np.complex128)
        for n in range(len(base_label)):
            label_idx = label == base_label[n]
            SFbase[:,n] = np.sum(SF[:,label_idx],axis=1)
        
        # Get the real part of the structure factor
        #SFtot = np.sqrt(np.real(SF * np.conj(SF)))
        SFrel = np.real(SFbase)
        SFimg = np.imag(SFbase)
        
        # Generate the results
        outstr = ''
        outstr+= '( h, k, l) Intensity' + ' '.join(['%12s    '%x for x in base_label])+'\n'
        for n in range(Nref):
            ss = ' '.join(['%6.1f + i%-6.1f' % (x,y) for x,y in zip(SFrel[n],SFimg[n])])
            outstr+= '(%2.0f,%2.0f,%2.0f) %9.2f    %s\n' % (HKL[n,0],HKL[n,1],HKL[n,2],I[n],ss)
        return outstr

    def print_symmetry_contributions(self,HKL):
        """
        Prints the symmetry contributions to the structure factor for each atomic site
        """

        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])
        Nref = len(HKL)

        # Calculate the full intensity
        I = self.intensity(HKL)

        # Calculate the structure factors of the symmetric atomic sites
        buvw, btype, base_label, bocc, buiso, bmxmymz = self.xtl.Atoms.get()
        operations = np.hstack([self.xtl.Symmetry.symmetric_coordinate_operations(buvw[n]) for n in range(len(buvw))])
        rotations = np.hstack(
            [self.xtl.Symmetry.symmetric_coordinate_operations(buvw[n], True)[1] for n in range(len(buvw))])

        # Calculate the structure factors
        uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        Qmag = self.xtl.Cell.Qmag(HKL)
        ff = fc.xray_scattering_factor(type, Qmag)
        dw = fc.debyewaller(uiso, Qmag)
        dot_KR = np.dot(HKL, uvw.T)
        phase = np.exp(1j * 2 * np.pi * dot_KR)
        sf = ff * dw * occ * phase

        # Generate the results
        outstr = ''
        for n in range(Nref):
            ss = '\n'
            all_phase = 0j
            all_sf = 0j
            for lab in base_label:
                label_idx = np.argwhere(label == lab)
                ss += '  %s\n' % lab
                tot_phase = 0j
                tot_sf = 0j
                for a in label_idx:
                    uvwstr = '(%-7.3g,%-7.3g,%-7.3g)' % (uvw[a, 0], uvw[a, 1], uvw[a, 2])
                    phstr = fg.complex2str(phase[n, a])
                    sfstr = fg.complex2str(sf[n, a])
                    val = (a, uvwstr, operations[a[0]], rotations[a[0]], phstr, sfstr)
                    ss += '    %3d %s %25s %20s  %s  %s\n' % val
                    tot_phase += phase[n, a]
                    tot_sf += sf[n, a]
                ss += '%74sTotal:  %s  %s\n' % (' ', fg.complex2str(tot_phase), fg.complex2str(tot_sf))
                all_phase += tot_phase
                all_sf += tot_sf
            ss += '%62s Reflection Total:  %s  %s\n' % (' ', fg.complex2str(all_phase), fg.complex2str(all_sf))
            outstr += '(%2.0f,%2.0f,%2.0f) I = %9.2f    %s\n' % (HKL[n, 0], HKL[n, 1], HKL[n, 2], I[n], ss)
        return outstr
    
    def find_close_reflections(self,HKL,energy_kev=None,max_twotheta=2,max_angle=10):
        """
        Find and print list of reflections close to the given one
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        all_HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        all_HKL = self.xtl.Cell.sort_hkl(all_HKL)
        
        all_tth = self.xtl.Cell.tth(all_HKL,energy_kev)
        tth1 = self.xtl.Cell.tth(HKL, energy_kev)
        tth_dif = np.abs(all_tth-tth1)
        
        all_Q = self.xtl.Cell.calculateQ(all_HKL)
        Q1 = self.xtl.Cell.calculateQ(HKL)
        all_angles = np.abs([fg.ang(Q1,Q2,'deg') for Q2 in all_Q])
        
        selected = (tth_dif < max_twotheta)*(all_angles < max_angle)
        sel_HKL = all_HKL[selected,:]
        sel_tth = all_tth[selected]
        sel_angles = all_angles[selected]
        sel_intensity = self.intensity(sel_HKL)
        
        # Generate Results
        fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'
        outstr = '' 
        
        outstr+= 'Energy = %6.3f keV\n' % energy_kev
        outstr+= 'Close to Reflection (%2.0f,%2.0f,%2.0f)\n' %(HKL[0],HKL[1],HKL[2])
        outstr+= '( h, k, l) TwoTheta    Angle  Intensity\n'
        count = 0
        for n in range(0,len(sel_HKL)):
            count += 1
            outstr+= fmt % (sel_HKL[n,0],sel_HKL[n,1],sel_HKL[n,2],sel_tth[n],sel_angles[n],sel_intensity[n])
        outstr+= 'Reflections: %1.0f\n'%count
        return outstr

    def multiple_scattering(self, hkl, azir=[0, 0, 1], pv=[1, 0], energy_range=[7.8, 8.2], numsteps=60,
                            full=False, pv1=False, pv2=False, sfonly=True, pv1xsf1=False):
        """
        Run multiple scattering code, plot result.

        mslist = xtl.Scatter.multiple_scattering([h,k,l], energy_range=[7.8, 8.2])

        See multiple_scattering.py for more details. Code created by Dr <NAME>, DLS
        :param hkl: [h,k,l] principle reflection
        :param azir: [h,k,l] reference of azimuthal 0 angle
        :param pv: [s,p] polarisation vector
        :param energy_range: [min, max] energy range in keV
        :param numsteps: int: number of calculation steps from energy min to max
        :param full: True/False: calculation type: full
        :param pv1: True/False: calculation type: pv1
        :param pv2: True/False: calculation type: pv2
        :param sfonly: True/False: calculation type: sfonly *default
        :param pv1xsf1: True/False: calculation type: pv1xsf1?
        :return: array
        """
        return ms.run_calcms(self.xtl, hkl, azir, pv, energy_range, numsteps,
                             full=full, pv1=pv1, pv2=pv2, sfonly=sfonly, pv1xsf1=pv1xsf1)

    def ms_azimuth(self, hkl, energy_kev, azir=[0, 0, 1], pv=[1, 0], numsteps=3, peak_width=0.1,
                   full=False, pv1=False, pv2=False, sfonly=True, pv1xsf1=False):
        """
        Returns an azimuthal dependence at a particular energy

        :param xtl: Crystal structure from Dans_Diffraction
        :param hkl: [h,k,l] principle reflection
        :param energy_kev: calculation energy
        :param azir: [h,k,l] reference of azimuthal 0 angle
        :param pv: [s,p] polarisation vector
        :param numsteps: int: number of calculation steps from energy min to max
        :param peak_width: convolution width
        :param full: True/False: calculation type: full
        :param pv1: True/False: calculation type: pv1
        :param pv2: True/False: calculation type: pv2
        :param sfonly: True/False: calculation type: sfonly *default
        :param pv1xsf1: True/False: calculation type: pv1xsf1?
        :return: None
        """

        energy_range = [energy_kev-0.001, energy_kev+0.001]
        mslist = self.xtl.Scatter.multiple_scattering(hkl, azir, pv, energy_range, numsteps,
                                                      full=full, pv1=pv1, pv2=pv2, sfonly=sfonly, pv1xsf1=pv1xsf1)

        if pv1 + pv2 + sfonly + full + pv1xsf1 != 0:
            azimuth = np.concatenate(mslist[:, [3, 4]])
            intensity = np.concatenate(mslist[:, [-1, -1]])
        else:
            azimuth = np.concatenate(mslist[:, [3, 4]])
            intensity = np.ones(azimuth)

        # create plotting mesh
        peak_width_pixels = 10
        pixels = int(360*(peak_width_pixels/ peak_width))
        mesh_azi = np.linspace(-180, 180, pixels)
        mesh = np.zeros(mesh_azi.shape)

        # add reflections to background
        pixel_coord = (azimuth-180) / 360.
        pixel_coord = (pixel_coord * pixels).astype(int)

        for n in range(len(intensity)):
            mesh[pixel_coord[n]] = mesh[pixel_coord[n]] + intensity[n]

        # Convolve with a gaussian (if >0 or not None)
        if peak_width:
            gauss_x = np.arange(-3 * peak_width_pixels, 3 * peak_width_pixels + 1)  # gaussian width = 2*FWHM
            G = fg.gauss(gauss_x, None, height=1, cen=0, fwhm=peak_width_pixels, bkg=0)
            mesh = np.convolve(mesh, G, mode='same')
        return mesh_azi, mesh

    '''  Removed tensor scattering 26/05/20 V1.7
    def tensor_scattering(self, atom_label, hkl, energy_kev=None, azir=[0, 0, 1], psideg=0, process='E1E1',
                          rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None):
        """
        Return tensor scattering intensities
          ss, sp, ps, pp = tensor_scattering('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param hkl: list/array, [h,k,l] reflection to calculate
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param psideg: float/array, azimuthal angle/ range
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :return: ss, sp, ps, pp intensity values
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        # Calculate tensor scattering
        ss, sp, ps, pp = ts.CalculateIntensityInPolarizationChannels(
            process, B, sitevec, sglist, lam, hkl, azir, psideg,
            K=rank, Time=time, Parity=parity, mk=mk, lk=lk, sk=sk
        )
        return ss, sp, ps, pp

    def tensor_scattering_stokes(self, atom_label, hkl, energy_kev=None, azir=[0, 0, 1], psideg=0, stokes=0,
                                 pol_theta=45, process='E1E1', rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None):
        """
        Return tensor scattering intensities for non-standard polarisation
          pol = tensor_scattering_stokes('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90, stokes=45)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param hkl: list/array, [h,k,l] reflection to calculate
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param psideg: float, azimuthal angle
        :param stokes: float/array, rotation of polarisation analyser (0=sigma), degrees
        :param pol_theta: float, scattering angle of polarisation analyser, degrees
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :return: array of intensity values
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        # Calculate tensor scattering
        stokesvec_swl = [0, 0, 1]
        pol = ts.CalculateIntensityFromPolarizationAnalyser(
            process, B, sitevec, sglist, lam, hkl, azir, psideg, stokes, pol_theta,
            stokesvec_swl, rank, time, parity, mk, lk, sk
        )
        return pol

    def print_tensor_scattering(self, atom_label, hkl, energy_kev=None, azir=[0, 0, 1], psideg=0, process='E1E1',
                                rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None):
        """
        Return tensor scattering intensities
          ss, sp, ps, pp = tensor_scattering('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param hkl: list/array, [h,k,l] reflection to calculate
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param psideg: float, azimuthal angle
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :return: str
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        # Calculate tensor scattering
        ss, sp, ps, pp = ts.CalculateIntensityInPolarizationChannels(
            process, B, sitevec, sglist, lam, hkl, azir, psideg,
            K=rank, Time=time, Parity=parity, mk=mk, lk=lk, sk=sk
        )
        outstr1 = ts.tensorproperties(sitevec, sglist, hkl, Parity=parity, Time=time)
        outstr2 = ts.print_tensors(B, sitevec, sglist, hkl, K=rank, Parity=parity, Time=time)
        outstr3 = "\nScattering Tensor:\n\n    [ss, sp] = [%5.2f, %5.2f]\n    [ps, pp]   [%5.2f, %5.2f]"
        outstr3 = outstr3 % (ss, sp, ps, pp)
        return outstr1 + outstr2 + outstr3

    def print_tensor_scattering_refs(self, atom_label, energy_kev=None, azir=[0, 0, 1], psideg=0, process='E1E1',
                                     rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None,
                                     print_symmetric=False, units=None):
        """
        Return tensor scattering intensities for all reflections at given azimuth and energy
          ss, sp, ps, pp = tensor_scattering('Ru1', 2.838, [0,1,0], psideg=90)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param psideg: float, azimuthal angle
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :param print_symmetric: False*/True : omits reflections with the same intensity at the same angle
        :param units: None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display
        :return: str
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        hkl = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        if not print_symmetric:
            hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)
        hkl = self.xtl.Cell.sort_hkl(hkl)

        tth = self.xtl.Cell.tth(hkl, energy_kev)
        inrange = np.all([tth < self._scattering_max_twotheta, tth > self._scattering_min_twotheta], axis=0)
        hkl = hkl[inrange, :]
        tth = tth[inrange]

        if units is None:
            units = self._powder_units
        units = units.lower()
        if units in ['d', 'dspc', 'dspace', 'd space', 'd-space', 'dspacing', 'd spacing', 'd-spacing']:
            unit_str = 'd-spacing'
            unit = fc.caldspace(tth, energy_kev)
        elif units in ['q', 'wavevector']:
            unit_str = 'Q'
            unit = fc.calqmag(tth, energy_kev)
        else:
            unit_str = 'TwoTheta'
            unit = tth

        fmt = '(%3.0f,%3.0f,%3.0f) %10.2f  %11.2f %11.2f %11.2f %11.2f\n'
        outstr = 'Tensor Scattering %s\n' % self.xtl.name
        outstr += 'Process: %s, site: %s\n' % (process, atom_label)
        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Psi_0 = (%.2g,%.2g,%.2g)  Psi = %3.3g\n' % (azir[0], azir[1], azir[2], psideg)
        outstr += '( h, k, l)    %10s  Sigma-Sigma    Sigma-Pi    Pi-Sigma       Pi-Pi\n' % unit_str

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        for n in range(1, len(tth)):
            ss, sp, ps, pp = ts.CalculateIntensityInPolarizationChannels(
                process, B, sitevec, sglist, lam, hkl[n, :], azir, psideg,
                K=rank, Time=time, Parity=parity, mk=mk, lk=lk, sk=sk
            )
            outstr += fmt % (hkl[n, 0], hkl[n, 1], hkl[n, 2], unit[n], ss, sp, ps, pp)
        outstr += 'Reflections: %1.0f\n' % len(tth)
        return outstr

    def print_tensor_scattering_refs_max(self, atom_label, energy_kev=None, azir=[0, 0, 1], process='E1E1',
                                     rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None,
                                     print_symmetric=False, units=None):
        """
        Return tensor scattering intensities for all reflections at given energy at maximum intensity psi
          ss, sp, ps, pp = tensor_scattering('Ru1', 2.838, [0,1,0], psideg=90)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :param print_symmetric: False*/True : omits reflections with the same intensity at the same angle
        :param units: None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display
        :return: str
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        hkl = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        if not print_symmetric:
            hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)
        hkl = self.xtl.Cell.sort_hkl(hkl)

        tth = self.xtl.Cell.tth(hkl, energy_kev)
        inrange = np.all([tth < self._scattering_max_twotheta, tth > self._scattering_min_twotheta], axis=0)
        hkl = hkl[inrange, :]
        tth = tth[inrange]
        # Caluclate structure factor **2
        inten = self.intensity(hkl)

        if units is None:
            units = self._powder_units
        units = units.lower()
        if units in ['d', 'dspc', 'dspace', 'd space', 'd-space', 'dspacing', 'd spacing', 'd-spacing']:
            unit_str = 'd-spacing'
            unit = fc.caldspace(tth, energy_kev)
        elif units in ['q', 'wavevector']:
            unit_str = 'Q'
            unit = fc.calqmag(tth, energy_kev)
        else:
            unit_str = 'TwoTheta'
            unit = tth

        fmt = '(%3.0f,%3.0f,%3.0f) %10.2f %10.2f  %5.2f (%3.0f) %5.2f (%3.0f) %5.2f (%3.0f) %5.2f (%3.0f)\n'
        outstr = 'Tensor Scattering %s\n' % self.xtl.name
        outstr += 'Process: %s, site: %s\n' % (process, atom_label)
        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Psi_0 = (%.2g,%.2g,%.2g)\n' % (azir[0], azir[1], azir[2])
        outstr += '( h, k, l)    %10s         I0  Sigma-Sigma    Sigma-Pi    Pi-Sigma       Pi-Pi\n' % unit_str

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        psideg = np.arange(0, 361, 5)
        for n in range(1, len(tth)):
            ss, sp, ps, pp = ts.CalculateIntensityInPolarizationChannels(
                process, B, sitevec, sglist, lam, hkl[n, :], azir, psideg,
                K=rank, Time=time, Parity=parity, mk=mk, lk=lk, sk=sk
            )
            iss = np.argmax(ss)
            isp = np.argmax(sp)
            ips = np.argmax(ps)
            ipp = np.argmax(pp)
            ss = ss[iss]
            sp = sp[isp]
            ps = ps[ips]
            pp = pp[ipp]
            psi_ss = psideg[iss]
            psi_sp = psideg[isp]
            psi_ps = psideg[ips]
            psi_pp = psideg[ipp]
            outstr += fmt % (hkl[n, 0], hkl[n, 1], hkl[n, 2], unit[n], inten[n],
                             ss, psi_ss, sp, psi_sp, ps, psi_ps, pp, psi_pp)
        outstr += 'Reflections: %1.0f\n' % len(tth)
        return outstr
    '''


class ScatteringTypes:
    """
    Container for available scattering types
    """
    def __init__(self, parent, typedict):
        self.parent = parent
        self.typedict = typedict

        typenames = typedict.keys()

        for name in typenames:
            attrname = name.replace(' ','_')
            setattr(self, attrname, self.ScatteringType(parent, name))

    def __call__(self, type):
        """
        Calling this container with a scattering type will set the scattering type in the parent class
        :param type: str scattering type
        :return: None
        """
        for typename in self.typedict:
            if type.lower() in self.typedict[typename]:
                self.parent._scattering_type = typename
                print('Changed scattering type to: %s'%typename)

    def i16(self):
        """
        Set max/min angles consistent with reflection geometry on beamline I16
        :return: None
        """
        self.parent.setup_scatter(type='xray', energy_kev=8.0, min_theta=-20, max_theta=140, min_twotheta=0, max_twotheta=140)

    def wish(self):
        """
        Set max/min angles consistent with reflection geometry on beamline I16
        :return: None
        """
        self.parent.setup_scatter(type='neutron', wavelength_a=1.5, min_theta=-180, max_theta=180, min_twotheta=-180, max_twotheta=180)

    class ScatteringType:
        """
        Container for scattering type switcher
        """
        def __init__(self, parent, typename):
            self.parent = parent
            self.typename = typename

        def __call__(self):
            self.parent._scattering_type = self.typename

<filename>Dans_Diffraction/classes_scattering.py
# -*- coding: utf-8 -*-
"""
Scattering Class "classes_scattering.py"
    Scattering functions for Crystal class.

By <NAME>, PhD
Diamond
2017

Version 2.0.1
Last updated: 28/09/21

Version History:
10/09/17 0.1    Program created
30/10/17 1.0    Main functions finshed, some testing done
06/01/18 1.1    Renamed classes_scattering.py
31/10/18 1.2    Added print_symmetric_contributions
21/01/19 1.3    Added non-resonant diffraction, corrected resonant diffraction
16/12/19 1.4    Added multiple_scattering code, print_all_reflections updated with units
18/02/20 1.5    Added tensor_scattering code
20/03/20 1.6    Increased powder gauss width from 2fwhm to 6fwhm, added powder averaging
14/04/20 1.6    Added powder_correction
26/05/20 1.7    Removed tensor_scattering
16/06/20 1.7.1  Added output option of setup_scatter
04/01/21 1.7.2  Added structure_factor function
21/01/21 1.8    Added xray_dispersion scattering function
10/06/21 1.9    Added x_ray calculation using Waasmaier and Kirfel scattering factors.
09/07/21 1.9    Added new scattering factors as option on normal scattering functions
20/08/21 2.0    Switched over to new scattering module, added self.powder()
28/09/21 2.0.1  Added __repr__

@author: DGPorter
"""

import sys,os
import numpy as np
import datetime

from . import functions_general as fg
from . import functions_crystallography as fc
from . import functions_scattering as fs
from . import multiple_scattering as ms
# from . import tensor_scattering as ts  # Removed V1.7

__version__ = '2.0.1'
__scattering_types__ = {'xray': ['xray','x','x-ray','thomson','charge'],
                        'neutron': ['neutron','n','nuclear'],
                        'xray magnetic': ['xray magnetic','magnetic xray','spin xray','xray spin'],
                        'neutron magnetic': ['neutron magnetic','magnetic neutron','magnetic'],
                        'xray resonant': ['xray resonant','resonant','resonant xray','rxs'],
                        'xray dispersion': ['dispersion', 'xray dispersion']}


class Scattering:
    """
    Simulate diffraction from Crystal
    Useage:
        xtl = Crystal()
        xtl.Scatter.setup_scatter(type='x-ray',energy_keV=8.0)
        xtl.Scatter.intensity([h,k,l]) # Returns intensity
        print(xtl.Scatter.print_all_refelctions()) # Returns formated string of all allowed reflections
        
        Allowed radiation types:
            'xray','neutron','xray magnetic','neutron magnetic','xray resonant'
    """
    
    #------Options-------
    # Standard Options
    _hkl = None  # added so recalculation not required
    _scattering_type = 'xray'  # 'xray','neutron','xray magnetic','neutron magnetic','xray resonant'
    _scattering_specular_direction = [0,0,1]  # reflection
    _scattering_parallel_direction = [0,0,1]  # transmission
    _scattering_theta_offset = 0.0
    _scattering_min_theta = -180.0
    _scattering_max_theta = 180.0
    _scattering_min_twotheta = -180.0
    _scattering_max_twotheta = 180.0
    _time_report = True
    _debug_mode = False

    # powder options
    _powder_units = 'tth' # tth (two theta), Q, d
    _powder_background = 0.0
    _powder_peak_width = 0.01  # in Deg
    _powder_average = True
    _powder_pixels = 2000
    
    # Complex Structure factor
    _return_structure_factor = False

    # Uses the coefficients for analytical approximation to the scattering factors from:
    #        "<NAME>, Acta Cryst. (1995) A51, 416-431"
    _use_waaskirf_scattering_factor = False
    
    # Thermal Factors
    _use_isotropic_thermal_factor = True
    _use_anisotropic_thermal_factor = False
    
    # Magnetic Options
    _calclate_magnetic_component = True
    _use_magnetic_form_factor = True
    
    # Polarisation Options
    _polarised = False
    _polarisation = 'sp'
    _polarisation_vector_incident = [0, 1, 0]
    
    # Radiation energy
    _energy_kev = fg.Cu
    
    # Resonant X-ray Options
    _azimuthal_angle = 0
    _azimuthal_reference = [1,0,0]
    _resonant_flm = (0, 1, 0)
    _resonant_approximation_e1e1 = True
    _resonant_approximation_e2e2 = False
    _resonant_approximation_e1e2 = False
    _resonant_approximation_m1m1 = False
    
    def __init__(self, xtl):
        self.xtl = xtl

        # Initialise the scattering type container
        self.Type = ScatteringTypes(self, __scattering_types__)

    def __repr__(self):
        return 'Scattering(%s, %s)' % (self.xtl.name, self._scattering_type)

    def __str__(self):
        out = 'Scatter(%r)\n' % self.xtl
        out += '       Type: %s\n' % self._scattering_type
        out += '     Energy: %s keV\n' % self._energy_kev
        out += ' Wavelength: %s A\n' % fc.energy2wave(self._energy_kev)
        out += ' ---Settings---\n'
        out += '      Powder units: %s\n' % self._powder_units
        out += '    Isotropic ADPs: %s\n' % self._use_isotropic_thermal_factor
        #out += '  Anisotropic ADPs: %s\n' % self._use_anisotropic_thermal_factor
        out += '  Specular Direction (reflection): (%2.0f,%2.0f,%2.0f)\n' % (
            self._scattering_specular_direction[0], self._scattering_specular_direction[1],
            self._scattering_specular_direction[2])
        out += 'Parallel Direction (transmission): (%2.0f,%2.0f,%2.0f)\n' % (
            self._scattering_parallel_direction[0], self._scattering_parallel_direction[1],
            self._scattering_parallel_direction[2])
        out += '      theta offset: %s\n' % self._scattering_theta_offset
        out += '         min theta: %s\n' % self._scattering_min_theta
        out += '         max theta: %s\n' % self._scattering_max_theta
        out += '      min twotheta: %s\n' % self._scattering_min_twotheta
        out += '      max twotheta: %s\n' % self._scattering_max_twotheta
        out += ' ---X-Ray Settings---\n'
        out += ' Waasmaier scattering factor: %s\n' % self._use_waaskirf_scattering_factor
        out += ' ---Magnetic Settings---\n'
        out += '   Mag. scattering: %s\n' % self._calclate_magnetic_component
        out += '  Mag. form factor: %s\n' % self._use_magnetic_form_factor
        out += '         Polarised: %s\n' % self._polarised
        out += '      Polarisation: %s\n' % self._polarisation
        out += '       Pol. vector: (%s,%s,%s)\n' % (self._polarisation_vector_incident[0],
                                                     self._polarisation_vector_incident[1],
                                                     self._polarisation_vector_incident[2])
        out += ' ---Resonant Settings---\n'
        out += '   Azimuthal angle: %s\n' % self._azimuthal_angle
        out += '    Azimuthal ref.: (%s,%s,%s)\n' % (self._azimuthal_reference[0], self._azimuthal_reference[1],
                                                     self._azimuthal_reference[2])
        out += '               flm: (%s,%s,%s)\n' % (self._resonant_flm[0], self._resonant_flm[1],
                                                     self._resonant_flm[2])
        out += '  use e1e1 approx.: %s\n' % self._resonant_approximation_e1e1
        #out += '  use e2e2 approx.: %s\n' % self._resonant_approximation_e2e2
        #out += '  use e1e2 approx.: %s\n' % self._resonant_approximation_e1e2
        #out += '  use m1m1 approx.: %s\n' % self._resonant_approximation_m1m1
        return out

    def setup_scatter(self, scattering_type=None, energy_kev=None, wavelength_a=None, powder_units=None,
                      specular=None, parallel=None, theta_offset=None,
                      min_theta=None, max_theta=None, min_twotheta=None, max_twotheta=None,
                      output=True, scattering_factors=None, magnetic_formfactor=None,
                      polarisation=None, polarisation_vector=None, azimuthal_reference=None, azimuth=None, flm=None):
        """
        Simple way to set scattering parameters, each parameter is internal to xtl (self)

        scattering_type: self._scattering type            :  'xray','neutron','xray magnetic','neutron magnetic','xray resonant', 'xray dispersion'
        energy_kev  : self._energy_kev                    :  radiation energy in keV
        wavelength_a: self._wavelength_a                  :  radiation wavelength in Angstrom
        powder_units: self._powder_units                  :  units to use when displaying/ plotting ['twotheta', 'd',' 'q']
        min_twotheta: self._scattering_min_two_theta      :  minimum detector (two-theta) angle
        max_twotheta: self._scattering_max_two_theta      :  maximum detector (two-theta) angle
        min_theta   : self._scattering_min_theta          :  minimum sample angle = -opening angle
        max_theta   : self._scattering_max_theta          :  maximum sample angle = opening angle
        theta_offset: self._scattering_theta_offset       :  sample offset angle
        specular    : self._scattering_specular_direction : [h,k,l] : reflections normal to sample surface
        parallel    : self._scattering_parallel_direction : [h,k,l] : reflections normal to sample surface
        scattering_factors: self._use_waaskirf_scattering_factor : xray scattering factor ['waaskirf', 'itc']
        magnetic_formfactor: self._use_magnetic_form_factor: True/False magnetic form factor for magnetic SF
        polarisation: self._polarisation                  : beam polarisation setting ['ss', 'sp'*, 'sp', 'pp']
        polarisation_vector: _polarisation_vector_incident: [x,y,z] incident polarisation vector
        azimuthal_reference: self._azimuthal_reference    : [h,k,l] direction of azimuthal zero angle
        azimuth    : self._azimuthal_angle                : azimuthal angle in deg
        flm        : self._resonant_flm                   : Resonant settings (flm1, flm2, flm3)
        """

        if scattering_type is not None:
            self._scattering_type = scattering_type

        if energy_kev is not None:
            self._energy_kev = energy_kev

        if wavelength_a is not None:
            self._energy_kev = fc.wave2energy(wavelength_a)

        if powder_units is not None:
            self._powder_units = powder_units

        if specular is not None:
            self._scattering_specular_direction = specular

        if parallel is not None:
            self._scattering_parallel_direction = parallel

        if theta_offset is not None:
            self._scattering_theta_offset = theta_offset

        if min_theta is not None:
            self._scattering_min_theta = min_theta

        if max_theta is not None:
            self._scattering_max_theta = max_theta

        if min_twotheta is not None:
            self._scattering_min_twotheta = min_twotheta

        if max_twotheta is not None:
            self._scattering_max_twotheta = max_twotheta

        if scattering_factors is not None:
            if scattering_factors.lower() in ['ws', 'waaskirf', 'alternate', 'alt']:
                print('Using scattering factors from: "Waasmaier and Kirfel, Acta Cryst. (1995) A51, 416-431"')
                self._use_waaskirf_scattering_factor = True
            else:
                print('Using scattering factors from: International Tables of Crystallography Vol. C, Table 6.1.1.4')
                self._use_waaskirf_scattering_factor = False

        if magnetic_formfactor is not None:
            self._use_magnetic_form_factor = magnetic_formfactor

        if polarisation_vector is not None:
            self._polarisation_vector_incident = np.array(polarisation_vector, dtype=np.float).reshape(3)

        if polarisation is not None:
            self._polarisation = polarisation

        if azimuthal_reference is not None:
            self._azimuthal_reference = np.array(azimuthal_reference, dtype=np.float).reshape(3)

        if azimuth is not None:
            self._azimuthal_angle = azimuth

        if flm is not None:
            self._resonant_flm = np.array(flm).reshape(3)

        if output:
            print(self)
            # print('Scattering Options:')
            # print('                            Type : %s' % self._scattering_type)
            # print('                  Default Energy : %6.3f keV' % self._energy_kev)
            # print('                    Powder Units : %s' % self._powder_units)
            # print('  Specular Direction (reflection): (%2.0f,%2.0f,%2.0f)' % (
            #     self._scattering_specular_direction[0], self._scattering_specular_direction[1],
            #     self._scattering_specular_direction[2]))
            # print('Parallel Direction (transmission): (%2.0f,%2.0f,%2.0f)' % (
            #     self._scattering_parallel_direction[0], self._scattering_parallel_direction[1],
            #     self._scattering_parallel_direction[2]))
            # print('                   Sample Offset : %5.2f' % self._scattering_theta_offset)
            # print('             Minimum Theta angle : %5.2f' % self._scattering_min_theta)
            # print('             Maximum Theta angle : %5.2f' % self._scattering_max_theta)
            # print('         Minimum Two-Theta angle : %5.2f' % self._scattering_min_twotheta)
            # print('         Maximum Two-Theta angle : %5.2f' % self._scattering_max_twotheta)

    def get_hkl(self, regenerate=True, remove_symmetric=False, reflection=False, transmission=False,  **kwargs):
        """
        Return stored hkl or generate
        :param regenerate: if True, hkl list will be regenerated, if False - previous list will be returned
        :param remove_symmetric: generate only non-symmetric hkl values
        :param reflection: generate only reflections possible in reflection geometry
        :param transmission: generate only reflections possible in transmission geometry
        :param kwargs: additional options to pass to setup_scatter()
        :return: array
        """
        if not regenerate and self._hkl is not None:
            return self._hkl
        self.setup_scatter(output=False, **kwargs)

        en = self._energy_kev
        max_tth = self._scattering_max_twotheta

        hkl = self.xtl.Cell.all_hkl(en, max_tth)
        if remove_symmetric:
            hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)
        hkl = self.xtl.Cell.sort_hkl(hkl)[1:]  # remove (0, 0, 0)

        if reflection:
            tth = self.xtl.Cell.tth(hkl, en)
            hkl = hkl[tth > self._scattering_min_twotheta, :]
            tth = tth[tth > self._scattering_min_twotheta]
            theta = self.xtl.Cell.theta_reflection(hkl, en, self._scattering_specular_direction,
                                                   self._scattering_theta_offset)
            p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)
            p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))
            hkl = hkl[p1 * p2]
        elif transmission:
            tth = self.xtl.Cell.tth(hkl, en)
            hkl = hkl[tth > self._scattering_min_twotheta, :]
            tth = tth[tth > self._scattering_min_twotheta]
            theta = self.xtl.Cell.theta_transmission(hkl, en, self._scattering_parallel_direction)

            p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)
            p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))
            hkl = hkl[p1 * p2]
        self._hkl = hkl
        return self._hkl

    def _debug(self):
        """Toggle on debug mode"""
        if self._debug_mode:
            print('Debug mode: off')
            self._debug_mode = False
            fs.DEBUG_MODE = False
        else:
            print('Debug mode: on')
            self._debug_mode = True
            fs.DEBUG_MODE = True

    def structure_factor(self, hkl=None, scattering_type=None, **kwargs):
        """
        Calculate the structure factor at reflection indexes (h,k,l)

        Notes:
        - Uses x-ray atomic form factors, calculated from approximated tables in the ITC
        - This may be a little slow for large numbers of reflections, as it is not currently
         possible to use accelerated calculation methods in Jython.
        - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs
        - Crystal.scale is used to scale the complex structure factor, so the intensity is
         reduced by (Crystal.scale)^2
        - Testing against structure factors calculated by Vesta.exe is exactly the same when using Waasmaier structure factors.

        :param hkl: array[n,3] : reflection indexes (h, k, l)
        :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']
        :param kwargs: additional options to pass to scattering function
        :return: complex array[n] : structure factors
        """

        if hkl is None:
            hkl = self.get_hkl()
        if scattering_type is None:
            scattering_type = self._scattering_type

        hkl = np.asarray(np.rint(hkl), dtype=np.float).reshape([-1, 3])
        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()

        q = self.xtl.Cell.calculateQ(hkl)
        r = self.xtl.Cell.calculateR(uvw)

        moment = fc.euler_moment(mxmymz, self.xtl.Cell.UV())
        azi_ref_q = self.xtl.Cell.calculateQ(self._azimuthal_reference)

        if 'energy_kev' in kwargs:
            energy_kev = kwargs.pop('energy_kev')
        elif 'wavelength_a' in kwargs:
            energy_kev = fc.wave2energy(kwargs.pop('wavelength_a'))
        else:
            energy_kev = self._energy_kev
        energy_kev = np.asarray(energy_kev, dtype=np.float).reshape(-1)
        nenergy = len(energy_kev)

        if 'psi' in kwargs:
            psi = kwargs.pop('psi')
        else:
            psi = self._azimuthal_angle
        psi = np.asarray(psi, dtype=np.float).reshape(-1)
        npsi = len(psi)

        options = fs.options(
            occ=occ,
            moment=moment,
            incident_polarisation_vector=self._polarisation_vector_incident,
            polarisation=self._polarisation,
            azi_ref_q=azi_ref_q,
            f0=self._resonant_flm[0],
            f1=self._resonant_flm[1],
            f2=self._resonant_flm[2],
        )
        options.update(kwargs)
        scattering_fun = fs.get_scattering_function(scattering_type)
        # print('Scattering function: %s' % scattering_fun.__name__)

        # Break up long lists of HKLs
        nref, natom = len(q), len(r)
        n_arrays = np.ceil(nref * natom / fs.MAX_QR_ARRAY)
        if n_arrays > 1:
            print('Splitting %d reflections (%d atoms) into %1.0f parts' % (nref, natom, n_arrays))
        q_array = np.array_split(q, n_arrays)
        sf = np.zeros([nref, nenergy, npsi], dtype=np.complex)
        start_time = datetime.datetime.now()
        for e, enval in enumerate(energy_kev):
            for p, psival in enumerate(psi):
                ls = 0
                for n, _q in enumerate(q_array):
                    if n_arrays > 1:
                        print(' Starting %2.0f/%2.0f: %d:%d' % (n + 1, n_arrays, ls, ls + len(_q)))
                    qmag = fg.mag(_q)
                    # Scattering factors
                    if scattering_type in fs.SCATTERING_TYPES['neutron']:
                        ff = fc.atom_properties(atom_type, 'Coh_b')
                    elif scattering_type in fs.SCATTERING_TYPES['xray fast']:
                        ff = fc.atom_properties(atom_type, 'Z')
                    elif scattering_type in fs.SCATTERING_TYPES['xray dispersion']:
                        ff = fc.xray_scattering_factor_resonant(atom_type, qmag, enval)
                    elif self._use_waaskirf_scattering_factor:
                        ff = fc.xray_scattering_factor_WaasKirf(atom_type, qmag)
                    else:
                        ff = fc.xray_scattering_factor(atom_type, qmag)

                    # Get Debye-Waller factor
                    if self._use_isotropic_thermal_factor:
                        dw = fc.debyewaller(uiso, qmag)
                    elif self._use_anisotropic_thermal_factor:
                        raise Exception('anisotropic thermal factor calcualtion not implemented yet')
                    else:
                        dw = None

                    # Get magnetic form factors
                    if self._use_magnetic_form_factor:
                        mf = fc.magnetic_form_factor(atom_type, qmag)
                    else:
                        mf = None

                    options['scattering_factor'] = ff
                    options['debyewaller'] = dw
                    options['magnetic_formfactor'] = mf
                    options['energy_kev'] = enval
                    options['psi'] = psival

                    sf[ls: ls + len(_q), e, p] = scattering_fun(_q, r, **options)
                    ls = ls + len(_q)

        end_time = datetime.datetime.now()
        time_difference = end_time - start_time
        if self._time_report and time_difference.total_seconds() > 10:
            print('Calculated %d structure factors in %s' % (nref, time_difference))
        sf = sf / self.xtl.scale
        if nenergy == 1 and npsi == 1:
            return sf[:, 0, 0]  # shape(nref)
        if nenergy == 1:
            return sf[:, 0, :]  # shape(nref, nenergy)
        if npsi == 1:
            return sf[:, :, 0]  # shape(nref, nspi)
        return sf
    new_structure_factor = structure_factor

    def intensity(self, hkl=None, scattering_type=None, **options):
        """

        :param hkl:
        :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']
        :return:
        """
        return fs.intensity(self.new_structure_factor(hkl, scattering_type, **options))
    new_intensity = intensity

    def powder(self, scattering_type=None, units=None, peak_width=None, background=None,
               pixels=None, powder_average=None, **options):
        """

        Generates array of intensities along a spaced grid, equivalent to a powder pattern.
          tth, iten = generate_powder('xray', units='tth', energy_kev=8)

        :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']
        :param units: str : one of ['tth', 'dspace', 'q']
        :param peak_width: float : Peak with in units of inverse wavevector (Q)
        :param background: float : if >0, a normal background around this value will be added
        :param pixels: int : number of pixels to add to the resulting mesh
        :param powder_average: Bool : if True, intensities will be reduced for the powder average
        :param options: additional arguments to pass to intensity calculation
        :return xval: arrray : x-axis of powder scan (units)
        :return yval: array :  intensity values at each point in x-axis
        """
        if scattering_type is None:
            scattering_type = self._scattering_type
        if units is None:
            units = self._powder_units
        if peak_width is None:
            peak_width = self._powder_peak_width
        if background is None:
            background = self._powder_background
        if pixels is None:
            pixels = self._powder_pixels
        if powder_average is None:
            powder_average = self._powder_average
        if 'energy_kev' in options:
            energy_kev = options['energy_kev']
        else:
            energy_kev = self._energy_kev

        # Units
        min_twotheta = self._scattering_min_twotheta
        if min_twotheta <= 0: min_twotheta = 1.0
        max_twotheta = self._scattering_max_twotheta
        q_min = fc.calqmag(min_twotheta, energy_kev)
        q_max = fc.calqmag(max_twotheta, energy_kev)

        # Get reflections
        hmax, kmax, lmax = fc.maxHKL(q_max, self.xtl.Cell.UVstar())
        HKL = fc.genHKL([hmax, -hmax], [kmax, -kmax], [lmax, -lmax])
        HKL = self.xtl.Cell.sort_hkl(HKL)  # required for labels
        Qmag = self.xtl.Cell.Qmag(HKL)
        HKL = HKL[Qmag < q_max, :]
        Qmag = self.xtl.Cell.Qmag(HKL)

        # Calculate intensities
        I = self.intensity(HKL, scattering_type, **options)

        if powder_average:
            # Apply powder averging correction, I0/|Q|**2
            I = I/(Qmag+0.001)**2

        # create plotting mesh
        pixels = int(pixels * q_max)  # reduce this to make convolution faster
        pixel_size = q_max / (1.0 * pixels)
        peak_width_pixels = peak_width / (1.0 * pixel_size)
        mesh = np.zeros([pixels])
        mesh_q = np.linspace(q_min, q_max, pixels)

        # add reflections to background
        pixel_coord = Qmag / (1.0 * q_max)
        pixel_coord = (pixel_coord * (pixels - 1)).astype(int)

        for n in range(1, len(I)):
            mesh[pixel_coord[n]] = mesh[pixel_coord[n]] + I[n]

        # Convolve with a gaussian (if >0 or not None)
        if peak_width:
            gauss_x = np.arange(-3*peak_width_pixels, 3*peak_width_pixels + 1)  # gaussian width = 2*FWHM
            G = fg.gauss(gauss_x, None, height=1, cen=0, fwhm=peak_width_pixels, bkg=0)
            mesh = np.convolve(mesh, G, mode='same')

        # Add background (if >0 or not None)
        if background:
            bkg = np.random.normal(background, np.sqrt(background), [pixels])
            mesh = mesh + bkg

        # Change output units
        if units.lower() in ['tth', 'angle', 'two-theta', 'twotheta', 'theta']:
            xval = fc.cal2theta(mesh_q, energy_kev)
        elif units.lower() in ['d', 'dspace', 'd-spacing', 'dspacing']:
            xval = fc.q2dspace(mesh_q)
        else:
            xval = mesh_q
        return xval, mesh

    def x_ray(self, HKL):
        """
        Calculate the squared structure factor for the given HKL, using x-ray scattering factors
          Scattering.x_ray([1,0,0])
          Scattering.x_ray([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1, 3])
        Nref = len(HKL)
        
        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        Nat = len(uvw)
        
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get atomic form factors
        if self._use_waaskirf_scattering_factor:
            ff = fc.xray_scattering_factor_WaasKirf(atom_type, Qmag)
        else:
            ff = fc.xray_scattering_factor(atom_type, Qmag)
        
        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso, Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calcualtion not implemented yet')
        else:
            dw = 1
        
        # Calculate dot product
        dot_KR = np.dot(HKL, uvw.T)
        
        # Calculate structure factor
        # Broadcasting used on 2D ff
        SF =  np.sum(ff*dw*occ*np.exp(1j*2*np.pi*dot_KR),axis=1)
        #SF = np.zeros(Nref,dtype=np.complex)
        #for ref in range(Nref):
        #    for at in range(Nat): 
        #        SF[ref] += ff[ref,at]*dw[ref,at]*occ[at]*np.exp(1j*2*np.pi*dot_KR[ref,at])
        
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def x_ray_fast(self, HKL):
        """
        Calculate the squared structure factor for the given HKL, using atomic number as scattering length
          Scattering.x_ray_fast([1,0,0])
          Scattering.x_ray_fast([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        
        uvw,type,label,occ,uiso,mxmymz = self.xtl.Structure.get()
        
        # Get atomic form factors
        ff = fc.atom_properties(type, 'Z')
        
        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            Qmag = self.xtl.Cell.Qmag(HKL)
            dw = fc.debyewaller(uiso,Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        # Calculate structure factor
        SF = np.sum(ff * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)
        
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)

    def xray_dispersion(self, HKL, energy_kev):
        """
        Calculate the squared structure factor for the given HKL,
        using x-ray scattering factors with dispersion corrections
          Scattering.xray_dispersion([1,0,0], 9.71)
          Scattering.xray_dispersion([[1,0,0],[2,0,0],[3,0,0]], 2.838)
          Scattering.xray_dispersion([[1,0,0],[2,0,0],[3,0,0]], np.arange(2.83, 2.86, 0.001))
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        if energy_kev is an array, the returned array will have shape [len(HKL), len(energy_kev)]
        """

        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])

        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()

        Qmag = self.xtl.Cell.Qmag(HKL)

        # Get atomic form factors
        ff = fc.xray_scattering_factor_resonant(atom_type, Qmag, energy_kev)  # shape (len(HKL), len(type), len(en))

        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso, Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calcualtion not implemented yet')
        else:
            dw = 1

        # Calculate dot product
        dot_KR = np.dot(HKL, uvw.T)

        energy_kev = np.asarray(energy_kev).reshape(-1)
        sf = np.zeros([len(HKL), len(energy_kev)], dtype=np.complex)
        for n, en in enumerate(energy_kev):
            # Broadcasting used on 2D ff
            sf[:, n] = np.sum(ff[:, :, n] * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)
        if len(energy_kev) == 1:
            sf = sf[:, 0]
        elif len(HKL) == 1:
            sf = sf[0, :]

        sf = sf / self.xtl.scale

        if self._return_structure_factor: return sf

        # Calculate intensity
        I = sf * np.conj(sf)
        return np.real(I)
    
    def neutron(self, HKL):
        """
        Calculate the squared structure factor for the given HKL, using neutron scattering length
          Scattering.neutron([1,0,0])
          Scattering.neutron([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        """
        
        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])

        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        
        # Get atomic form factors
        ff = fc.atom_properties(atom_type, 'Coh_b')
        
        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            Qmag = self.xtl.Cell.Qmag(HKL)
            dw = fc.debyewaller(uiso,Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        # Calculate structure factor
        SF = np.sum(ff * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)
        
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def magnetic_neutron(self, HKL):
        """
        Calculate the magnetic component of the structure factor for the given HKL, using neutron rules and form factor
          Scattering.magnetic_neutron([1,0,0])
          Scattering.magnetic_neutron([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)

        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        Nat = len(uvw)
        
        Q = self.xtl.Cell.calculateQ(HKL)
        Qhat = fg.norm(Q).reshape([-1,3])
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get magnetic form factors
        if self._use_magnetic_form_factor:
            ff = fc.magnetic_form_factor(atom_type, Qmag)
        else:
            ff = np.ones([len(HKL), Nat])
        
        # Calculate moment
        momentmag = fg.mag(mxmymz).reshape([-1,1])
        momentxyz = self.xtl.Cell.calculateR(mxmymz)
        moment = momentmag*fg.norm(momentxyz) # broadcast n*1 x n*3 = n*3
        moment[np.isnan(moment)] = 0.
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)

        # Calculate structure factor
        SF = np.zeros(Nref,dtype=np.complex)
        for n,Qh in enumerate(Qhat):
            SFm = [0.,0.,0.]
            for m,mom in enumerate(moment):
                # Calculate Magnetic part
                QM = mom - np.dot(Qh,mom)*Qh

                # Calculate structure factor
                SFm = SFm + ff[n, m] * np.exp(1j * 2 * np.pi * dot_KR[n, m]) * QM
            
            # Calculate polarisation with incident neutron
            if self._polarised:
                SF[n] = np.dot(SFm, self._polarisation_vector_incident)
            else:
                #SF[n] = np.dot(SFm,SFm) # maximum possible
                SF[n] = (np.dot(SFm,[1,0,0]) + np.dot(SFm,[0,1,0]) + np.dot(SFm,[0,0,1]))/3 # average polarisation

        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def xray_magnetic(self, HKL):
        """
        Calculate the non-resonant magnetic component of the structure factor 
        for the given HKL, using x-ray rules and form factor
          Scattering.xray_magnetic([1,0,0])
          Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        
        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)
        Book: "X-ray Scattering and Absorption by Magnetic Materials" by <NAME> Collins. Ch 2. Eqn.2.21+1
        No orbital component assumed
        magnetic moments assumed to be in the same reference frame as the polarisation
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)
        
        uvw,type,label,occ,uiso,mxmymz = self.xtl.Structure.get()
        Nat = len(uvw)
        
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get magnetic form factors
        if self._use_magnetic_form_factor:
            ff = fc.magnetic_form_factor(type,Qmag)
        else:
            ff = np.ones([len(HKL),Nat])
        
        # Calculate moment
        momentmag = fg.mag(mxmymz).reshape([-1,1])
        momentxyz = self.xtl.Cell.calculateR(mxmymz) # moment direction in cartesian reference frame
        moment = momentmag*fg.norm(momentxyz) # broadcast n*1 x n*3 = n*3
        moment[np.isnan(moment)] = 0.
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        # Calculate structure factor
        SF = np.zeros(Nref,dtype=np.complex)
        for n in range(Nref):
            # Calculate vector structure factor
            SFm = [0.,0.,0.]
            for m,mom in enumerate(moment):
                SFm = SFm + ff[n,m]*np.exp(1j*2*np.pi*dot_KR[n,m])*mom
            
            # Calculate polarisation with incident x-ray
            # The reference frame of the x-ray and the crystal are assumed to be the same
            # i.e. pol=[1,0,0] || mom=[1,0,0] || (1,0,0)
            if self._polarised:
                SF[n] = np.dot(SFm,self._polarisation_vector_incident)
            else:
                #SF[n] = np.dot(SFm,SFm) # maximum possible
                SF[n] = (np.dot(SFm,[1,0,0]) + np.dot(SFm,[0,1,0]) + np.dot(SFm,[0,0,1]))/3 # average polarisation
        
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def xray_resonant(self, HKL, energy_kev=None, polarisation='sp', F0=1, F1=1, F2=1, azim_zero=[1,0,0], PSI=[0], disp=False):
        """
        Calculate structure factors using resonant scattering factors in the dipolar approximation
          I = Scattering.xray_resonant(HKL,energy_kev,polarisation,F0,F1,F2)
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
            energy_kev = x-ray energy in keV
            polarisation = x-ray polarisation: 'ss',{'sp'},'ps','pp'
            F0/1/2 = Resonance factor Flm
            azim_zero = [h,k,l] vector parallel to the origin of the azimuth
            psi = azimuthal angle defining the scattering plane
        
        Uses the E1E1 resonant x-ray scattering amplitude:
            fxr_n = (ef.ei)*F0 -i(ef X ei).z_n*F1 + (ef.z_n)(ei.z_n)F2
        
        Where ei and ef are the initial and final polarisation states, respectively,
        and z_n is a unit vector in the direction of the magnetic moment of the nth ion.
        The polarisation states are determined to be one of the natural synchrotron 
        states, where sigma (s) is perpendicular to the scattering plane and pi (p) is
        parallel to it.
                ( s-s  s-p )
                ( p-s  p-p )
        
        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (15)
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)
        
        PSI = np.asarray(PSI,dtype=np.float).reshape([-1])
        Npsi = len(PSI)

        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso,Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        SF = np.zeros([Nref,Npsi],dtype=np.complex)
        for psival in range(Npsi):
            # Get resonant form factor
            fxres = self.xray_resonant_scattering_factor(HKL,energy_kev,polarisation,F0,F1,F2,azim_zero,PSI[psival],disp=disp)
            
            # Calculate structure factor
            # Broadcasting used on 2D fxres
            SF[:,psival] = np.sum(fxres*dw*occ*np.exp(1j*2*np.pi*dot_KR),axis=1)
            
        SF = SF/self.xtl.scale
        
        if self._return_structure_factor: return SF
        
        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)
    
    def xray_resonant_scattering_factor(self, HKL, energy_kev=None, polarisation='sp', F0=1, F1=1, F2=1, azim_zero=[1,0,0], psi=0, disp=False):
        """
        Calcualte fxres, the resonant x-ray scattering factor
          fxres = Scattering.xray_resonant_scattering_factor(HKL,energy_kev,polarisation,F0,F1,F2,azim_zero,psi)
        energy_kev = x-ray energy in keV
            polarisation = x-ray polarisation: 'ss',{'sp'},'ps','pp'
            F0/1/2 = Resonance factor Flm
            azim_zero = [h,k,l] vector parallel to the origin of the azimuth {[1,0,0]}
            psi = azimuthal angle defining the scattering plane {0}
        
        Uses the E1E1 resonant x-ray scattering amplitude:
            fxr_n = (ef.ei)*F0 -i(ef X ei).z_n*F1 + (ef.z_n)(ei.z_n)F2
        
        Where ei and ef are the initial and final polarisation states, respectively,
        and z_n is a unit vector in the direction of the magnetic moment of the nth ion.
        The polarisation states are determined to be one of the natural synchrotron 
        states, where sigma (s) is perpendicular to the scattering plane and pi (p) is
        parallel to it.
                ( s-s  s-p )
                ( p-s  p-p )
        
        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (15)
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)
        Qmag = self.xtl.Cell.Qmag(HKL)
        tth = fc.cal2theta(Qmag, energy_kev)
        
        mxmymz = self.xtl.Structure.mxmymz()
        Nat = len(mxmymz)
        
        fxres = np.zeros([Nref,Nat],dtype=np.complex)
        for ref in range(Nref):
            # Resonant scattering factor
            # Electric Dipole transition at 3d L edge
            #F0,F1,F2 = 1,1,1 # Flm (form factor?)
            z1,z2,z3 = self.scatteringcomponents(mxmymz, HKL[ref], azim_zero, psi).T
            tthr = np.deg2rad(tth[ref])/2.0

            polarisation = polarisation.replace('-', '').replace(' ', '')
            if polarisation in ['sigmasigma', 'sigsig', 'ss']:    # Sigma-Sigma
                f0 = 1*np.ones(Nat)
                f1 = 0*np.ones(Nat)
                f2 = z2**2
            elif polarisation in ['sigmapi', 'sigpi', 'sp']:  # Sigma-Pi
                f0 = 0 * np.ones(Nat)
                f1 = z3 * np.sin(tthr) + z1 * np.cos(tthr)
                f2 = z2 * (z1 * np.sin(tthr) + z3 * np.cos(tthr))
            elif polarisation in ['pisigma', 'pisig', 'ps']:  # Pi-Sigma
                f0 = 0*np.ones(Nat)
                f1 = z1*np.cos(tthr) - z3*np.sin(tthr)
                f2 = -z2*(z1*np.sin(tthr)-z3*np.cos(tthr))
            elif polarisation in ['pipi', 'pp']:  # Pi-Pi
                f0 = np.cos(2*tthr)*np.ones(Nat)
                f1 = -z2*np.sin(2*tthr)
                f2 = -(np.cos(tthr)**2)*(z1**2*np.tan(tthr)**2 + z3**2)
            else:
                raise ValueError('Incorrect polarisation. pol should be e.g. ''ss'' or ''sp''')
            fxres[ref,:] = F0*f0 -1j*F1*f1 + F2*f2
            if disp:
                print('( h, k, l)   TTH  (    mx,    my,    mz)  (    z1,    z2,    z3)')
                fmt = '(%2.0f,%2.0f,%2.0f) %6.2f  (%6.3f,%6.3f,%6.3f)  (%6.3f,%6.3f,%6.3f)  f0=%8.4f  f1=%8.4f  f2=%8.4f fxres= (%8.4f + %8.4fi)'
                for at in range(Nat):
                    vals = (HKL[ref,0],HKL[ref,1],HKL[ref,2],tth[ref],
                            mxmymz[at,0],mxmymz[at,1],mxmymz[at,2],
                            z1[at],z2[at],z3[at],
                            f0[at],f1[at],f2[at],
                            fxres[ref,at].real,fxres[ref,at].imag)
                    print(fmt%vals)
        return fxres

    def xray_nonresonant_magnetic(self, HKL, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p', disp=False):
        """
        Calculate the non-resonant magnetic component of the structure factor
        for the given HKL, using x-ray rules and form factor
          Scattering.xray_magnetic([1,0,0])
          Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.

        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)
        Book: "X-ray Scattering and Absorption by Magnetic Materials" by Loevesy and Collins. Ch 2. Eqn.2.21+1
        No orbital component assumed
        magnetic moments assumed to be in the same reference frame as the polarisation
        """

        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])

        uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()

        kin, kout, ein, eout = self.scatteringvectors(HKL, energy_kev, azim_zero, psi, polarisation)

        Qmag = self.xtl.Cell.Qmag(HKL)

        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso, Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1

        # Get magnetic form factors
        if self._use_magnetic_form_factor:
            ff = fc.magnetic_form_factor(type, Qmag)
        else:
            ff = np.ones([len(HKL), len(uvw)])

        # Calculate moment
        momentmag = fg.mag(mxmymz).reshape([-1, 1])
        momentxyz = self.xtl.Cell.calculateR(mxmymz)  # moment direction in cartesian reference frame
        moment = momentmag * fg.norm(momentxyz)  # broadcast n*1 x n*3 = n*3
        moment[np.isnan(moment)] = 0.

        # Magnetic form factor
        # f_non-res_mag = i.r0.(hw/mc^2).fD.[.5.L.A + S.B] #equ 2 Hill+McMorrow 1996
        # ignore orbital moment L
        B = np.zeros([len(HKL), 3])
        for n in range(len(HKL)):
            #print(n,HKL[n],kin[n],kout[n],ein[n],eout[n])
            B[n, :] = np.cross(eout[n], ein[n]) + \
                np.cross(kout[n], eout[n]) * np.dot(kout[n], ein[n]) - \
                np.cross(kin[n], ein[n]) * np.dot(kin[n], eout[n]) - \
                np.cross(np.cross(kout[n], eout[n]), np.cross(kin[n], ein[n]))
        fspin = 1j * ff * np.dot(moment, B.T).T

        # Calculate dot product
        dot_KR = np.dot(HKL, uvw.T)

        # Calculate structure factor
        SF = np.sum(fspin * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)

        if disp:
            for n in range(len(HKL)):
                pr_h = '(%2.0f,%2.0f,%2.0f)'%tuple(HKL[n,:])
                pr_b = '(%6.2g,%6.2g,%6.2g)'%tuple(B[n,:])
                print('psi=%3.0f  %3d hkl=%10s B=%22s'%(psi, n, pr_h, pr_b))
                ctot = 0j
                for m in range(len(moment)):
                    if np.sum(moment[m,:]**2) < 0.01: continue
                    pr_m = '(%6.2g,%6.2g,%6.2g)'%tuple(moment[m,:])
                    dot = np.dot(moment[m],B[n,:])
                    pdot = '%6.3f'%dot
                    phase = np.exp(1j * 2 * np.pi * dot_KR[n, m])
                    prph = '%5.2f+i%5.2f'%(np.real(phase),np.imag(phase))
                    prod = dot*phase
                    pprd = '%5.2f+i%5.2f'%(np.real(prod),np.imag(prod))
                    ctot += prod
                    ptot = '%5.2f+i%5.2f'%(np.real(ctot),np.imag(ctot))
                    print('\t%3d mom=%22s dot(mom,B)=%6s   exp(ik.r)=%12s   sum=%12s   tot=%12s'%(m,pr_m,pdot,prph,pprd,ptot))

        SF = SF / self.xtl.scale

        if self._return_structure_factor: return SF

        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)

    def xray_resonant_magnetic(self, HKL, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p', F0=0, F1=1, F2=0, disp=True):
        """
        Calculate the non-resonant magnetic component of the structure factor
        for the given HKL, using x-ray rules and form factor
          Scattering.xray_magnetic([1,0,0])
          Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.

        From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)
        Book: "X-ray Scattering and Absorption by Magnetic Materials" by <NAME> Collins. Ch 2. Eqn.2.21+1
        No orbital component assumed
        magnetic moments assumed to be in the same reference frame as the polarisation
        """

        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])

        uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()

        kin, kout, ein, eout = self.scatteringvectors(HKL, energy_kev, azim_zero, psi, polarisation)

        Qmag = self.xtl.Cell.Qmag(HKL)

        # Get Debye-Waller factor
        if self._use_isotropic_thermal_factor:
            dw = fc.debyewaller(uiso, Qmag)
        elif self._use_anisotropic_thermal_factor:
            raise Exception('anisotropic thermal factor calculation not implemented yet')
        else:
            dw = 1

        # Calculate moment
        momentmag = fg.mag(mxmymz).reshape([-1, 1])
        momentxyz = self.xtl.Cell.calculateR(mxmymz)  # moment direction in cartesian reference frame
        moment = momentmag * fg.norm(momentxyz)  # broadcast n*1 x n*3 = n*3
        moment[np.isnan(moment)] = 0.
        fe1e1 = np.zeros([len(HKL), len(uvw)], dtype=np.complex)
        for ref in range(len(HKL)):
            # Magnetic form factor
            # f_res_mag = [(e'.e)F0 - i(e'xe).Z*F1 + (e'.Z)*(e.Z)*F2] #equ 7 Hill+McMorrow 1996
            f0 = np.dot(eout[ref], ein[ref])
            f1 = np.dot(np.cross(eout[ref], ein[ref]), moment.T)
            f2 = np.dot(eout[ref], moment.T) * np.dot(ein[ref], moment.T)
            fe1e1[ref, :] = f0*F0 - 1j*f1*F1 + f2*F2

        # Calculate dot product
        dot_KR = np.dot(HKL, uvw.T)

        if disp:
            for ref in range(len(HKL)):
                pr_h = '(%2.0f,%2.0f,%2.0f)'%tuple(HKL[ref,:])
                print('psi=%3.0f  %3d hkl=%10s'%(psi, ref, pr_h))
                ctot = 0j
                for m in range(len(moment)):
                    if np.sum(moment[m,:]**2) < 0.01: continue
                    pr_m = '(%6.2g,%6.2g,%6.2g)'%tuple(moment[m,:])
                    f0 = np.dot(eout[ref], ein[ref])
                    f1 = np.dot(np.cross(eout[ref], ein[ref]), moment[m,:])
                    f2 = np.dot(eout[ref], moment[m,:]) * np.dot(ein[ref], moment[m,:])
                    dot = f0*F0 - 1j*f1*F1 + f2*F2
                    phase = np.exp(1j * 2 * np.pi * dot_KR[ref, m])
                    prph = '%5.2f+i%5.2f'%(np.real(phase),np.imag(phase))
                    prod = dot*phase
                    pprd = '%5.2f+i%5.2f'%(np.real(prod),np.imag(prod))
                    ctot += prod
                    ptot = '%5.2f+i%5.2f'%(np.real(ctot),np.imag(ctot))
                    print('\t%3d mom=%22s f0=%5.2f f1=%5.2f f2=%5.2f   exp(ik.r)=%12s   sum=%12s   tot=%12s'%(m,pr_m,f0,f1,f2,prph,pprd,ptot))

        # Calculate structure factor
        SF = np.sum(fe1e1 * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)

        SF = SF / self.xtl.scale

        if self._return_structure_factor: return SF

        # Calculate intensity
        I = SF * np.conj(SF)
        return np.real(I)

    def scatteringvectors(self, hkl, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p'):
        """
        Determine the scattering and polarisation vectors of a reflection based on energy, azimuth and polarisation.
        :param xtl: Crystal Class
        :param hkl: [n,3] array of reflections
        :param energy_kev: x-ray scattering energy in keV
        :param azim_zero: [1,3] direction along which the azimuthal zero angle is determind
        :param psi: float angle in degrees about the azimuth
        :param polarisation: polarisation with respect to the scattering plane, options:
                    'ss' : sigma-sigma polarisation
                    'sp' : sigma-pi polarisation
                    'ps' : pi-sigma polarisation
                    'pp' : pi-pi polarisation
                or: polarisation: float polarisation angle of scattered vector in degrees
        :return: kin, kout, ein, eout
        Returned values are [n,3] arrays
            kin : [n,3] array of incident wavevectors
            kout: [n,3] array of scattered wavevectors
            ein : [n,3] array of incident polarisation
            eout: [n,3] array of scattered polarisation

        The basis is chosen such that Q defines the scattering plane, sigma and pi directions are normal to this plane.
        Q is defined as Q = kout - kin, with kout +ve along the projection of azim_zero
        """

        hkl = np.asarray(hkl).reshape([-1, 3])

        if energy_kev is None:
            energy_kev = self._energy_kev

        out_kin = np.zeros([len(hkl), 3])
        out_kout = np.zeros([len(hkl), 3])
        out_ein = np.zeros([len(hkl), 3])
        out_eout = np.zeros([len(hkl), 3])
        for n in range(len(hkl)):
            # Define coordinate system I,J,Q (U1,U2,U3)
            Ihat, Jhat, Qhat = self.scatteringbasis(hkl[n, :], azim_zero, psi)
            Ihat = Ihat.reshape([-1, 3])
            Jhat = Jhat.reshape([-1, 3])
            Qhat = Qhat.reshape([-1, 3])

            # Determine wavevectors
            bragg = self.xtl.Cell.tth(hkl[n, :], energy_kev) / 2.
            rb = np.deg2rad(bragg).reshape([-1, 1])
            kin = np.cos(rb) * Ihat - np.sin(rb) * Qhat
            kout = np.cos(rb) * Ihat + np.sin(rb) * Qhat
            esig = Jhat  # sigma polarisation (in or out)
            piin = np.cross(kin, esig)  # pi polarisation in
            piout = np.cross(kout, esig)  # pi polarisation out

            # Polarisation
            try:
                # polarisation = 'ss' or 's-s'
                polarisation = polarisation.replace('-', '').replace(' ', '')
            except AttributeError:
                # polarisation = angle in deg from sigma' to pi'
                ein = 1.0 * esig
                pol = np.deg2rad(polarisation)
                eout = np.cos(pol)*esig + np.sin(pol)*piout
            if polarisation in ['sigmasigma', 'sigsig', 'ss']:
                ein = 1.0 * esig
                eout = 1.0 * esig
            elif polarisation in ['sigmapi', 'sigpi', 'sp']:
                ein = 1.0 * esig
                eout = 1.0 * piout
            elif polarisation in ['pisigma', 'pisig', 'ps']:
                ein = 1.0 * piin
                eout = 1.0 * esig
            elif polarisation in ['pipi', 'pp']:
                ein = 1.0 * piin
                eout = 1.0 * piout
            out_kin[n, :] = kin
            out_kout[n, :] = kout
            out_ein[n, :] = ein
            out_eout[n, :] = eout

        return out_kin, out_kout, out_ein, out_eout

    def scatteringcomponents(self, mxmymz, hkl, azim_zero=[1,0,0], psi=0):
        """
        Transform magnetic vector into components within the scattering plane
            ***warning - may not be correct for non-cubic systems***
        """

        # Define coordinate system I,J,Q (U1,U2,U3)
        U = self.scatteringbasis(hkl, azim_zero, psi)
        
        # Determine components of the magnetic vector
        z1z2z3 = np.dot(mxmymz, U.T) # [mxmymz.I, mxmymz.J, mxmymz.Q]
        return fg.norm(z1z2z3)

    def scatteringbasis(self, hkl, azim_zero=[1, 0, 0], psi=0):
        """
        Determine the scattering and polarisation vectors of a reflection based on energy, azimuth and polarisation.
        :param hkl: [n,3] array of reflections
        :param azim_zero: [1,3] direction along which the azimuthal zero angle is determind
        :param psi: float azimuthal angle about U3 in degrees
        :return: U1, U2, U3
        The basis is chosen such that Q defines the scattering plane, the sigma direction is normal to this plane,
        the pi direction is always within this plane.
        The azimuthal angle defines a rotation about the Q axis in a clockwise mannor, matching I16.
        At an azimuth of 0degrees, U1 is perpendicular to Q, along the direction of azim_zero.
        """

        # Define coordinate system I,J,Q (U1,U2,U3)
        # See FDMNES User's Guide p20 'II-11) Anomalous or resonant diffraction'
        # U1 || projection of azim_zero
        # U2 _|_ U1,U3
        # U3 || Q = kf-ki
        azim_zero = fg.norm(self.xtl.Cell.calculateQ(azim_zero)) # put in orthogonal basis
        Qhat = fg.norm(self.xtl.Cell.calculateQ(hkl)).reshape([-1,3])  # || Q
        AxQ = fg.norm(np.cross(azim_zero, Qhat))
        Ihat = fg.norm(np.cross(Qhat, AxQ)).reshape([-1,3])  # || to projection of azim_zero
        Jhat = fg.norm(np.cross(Qhat, Ihat)).reshape([-1,3])  # _|_ to I and Q

        # Rotate psi about Qhat
        rpsi = np.deg2rad(psi)
        # -ve sin makes clockwise rotation
        # This was checked on 21/1/19 vs CRO paper + sergio's calculations and seems to agree with experiment,
        # however we never did an azimuthal scan of the (103) which would have distinguished this completely.
        Ihat_psi = fg.norm(np.cos(rpsi) * Ihat - np.sin(rpsi) * Jhat)
        Jhat_psi = fg.norm(np.cross(Qhat, Ihat_psi))
        return np.vstack([Ihat_psi, Jhat_psi, Qhat])

    def print_scattering_coordinates(self,hkl,azim_zero=[1,0,0],psi=0):
        """
        Transform magnetic vector into components within the scattering plane
            ***warning - may not be correct for non-cubic systems***
        """
        
        # Define coordinate system I,J,Q (U1,U2,U3)
        Qhat = fg.norm(self.xtl.Cell.calculateQ(hkl)) # || Q
        AxQ = fg.norm(np.cross(azim_zero,Qhat))
        Ihat = fg.norm(np.cross(Qhat,AxQ)) # || to azim_zero
        Jhat = fg.norm(np.cross(Qhat,Ihat)) # -| to I and Q
        
        # Rotate coordinate system by azimuth
        Ihat_psi = fg.norm(np.cos(np.deg2rad(psi))*Ihat + np.sin(np.deg2rad(psi))*Jhat)
        Jhat_psi = fg.norm(np.cross(Qhat,Ihat_psi))
        
        # Determine components of the magnetic vector
        U=np.vstack([Ihat_psi,Jhat_psi,Qhat])
        print('U1 = (%5.2f,%5.2f,%5.2f)'%(U[0,0],U[0,1],U[0,2]))
        print('U2 = (%5.2f,%5.2f,%5.2f)'%(U[1,0],U[1,1],U[1,2]))
        print('U3 = (%5.2f,%5.2f,%5.2f)'%(U[2,0],U[2,1],U[2,2]))
    
    def print_intensity(self, HKL):
        """
        Print intensities calcualted in different ways
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Qmag =  self.xtl.Cell.Qmag(HKL)
        srt = np.argsort(Qmag)
        HKL = HKL[srt,:]
        
        IN=self.neutron(HKL)
        IX=self.x_ray(HKL)
        INM=self.magnetic_neutron(HKL)*1e4
        IXM=self.xray_magnetic(HKL)*1e4
        IXRss=self.xray_resonant(HKL, None, 'ss')
        IXRsp=self.xray_resonant(HKL, None, 'sp')
        IXRps=self.xray_resonant(HKL, None, 'ps')
        IXRpp=self.xray_resonant(HKL, None, 'pp')
        
        fmt = '(%2.0f,%2.0f,%2.0f)  %8.1f  %8.1f  %8.2f  %8.2f  ss=%8.2f  sp=%8.2f  ps=%8.2f  pp=%8.2f'
        print('( h, k, l)   Neutron      xray   Magn. N  Magn. XR   sig-sig    sig-pi    pi-sig     pi-pi')
        for n in range(len(HKL)):
            vals=(HKL[n][0],HKL[n][1],HKL[n][2],IN[n],IX[n],INM[n],IXM[n],IXRss[n],IXRsp[n],IXRps[n],IXRpp[n])
            print(fmt%vals)
    
    def old_intensity(self, HKL, scattering_type=None):
        """
        Calculate the squared structure factor for the given HKL
          Crystal.intensity([1,0,0])
          Crystal.intensity([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the real intensity at each reflection.
        
        Notes:
        - Uses x-ray atomic form factors, calculated from approximated tables in the ITC
        - This may be a little slow for large numbers of reflections, as it is not currently
         possible to use accelerated calculation methods in Jython.
        - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs
        - Crystal.scale is used to scale the complex structure factor, so the intensity is  
         reduced by (Crystal.scale)^2
        - Testing against structure factors calculated by Vesta.exe is very close, though there
          are some discrepancies, probably due to the method of calculation of the form factor.
        """

        if scattering_type is None:
            scattering_type = self._scattering_type
        scattering_type = scattering_type.lower()

        # Break up long lists of HKLs
        n_arrays = np.ceil(len(HKL)*len(self.xtl.Structure.u)/10000.)
        hkl_array = np.array_split(HKL, n_arrays)

        intensity = []
        for _hkl in hkl_array:
            if scattering_type in ['xray','x','x-ray','thomson','charge']:
                intensity += self.x_ray(_hkl).tolist()
            elif scattering_type in ['neutron','n','nuclear']:
                intensity += self.neutron(_hkl).tolist()
            elif scattering_type in ['xray magnetic','magnetic xray','spin xray','xray spin']:
                intensity += list(self.xray_magnetic(_hkl)*1e4)
            elif scattering_type in ['neutron magnetic','magnetic neutron','magnetic']:
                intensity += list(self.magnetic_neutron(_hkl)*1e4)
            elif scattering_type in ['xray dispersion']:
                intensity += self.xray_dispersion(_hkl, self._energy_kev).tolist()
            elif scattering_type in ['xray resonant','resonant','resonant xray','rxs']:
                intensity += self.xray_resonant(_hkl).tolist()
            elif scattering_type in ['xray resonant magnetic', 'xray magnetic resonant',
                                     'resonant magnetic', 'magnetic resonant']:
                intensity += self.xray_resonant_magnetic(
                    _hkl,
                    self._energy_kev,
                    self._azimuthal_reference,
                    self._azimuthal_angle,
                    self._polarisation,
                    F0=0, F1=1, F2=0).tolist()
            elif scattering_type in ['xray nonresonant magnetic', 'xray magnetic nonresonant',
                                     'nonresonant magnetic', 'magnetic nonresonant',
                                     'xray non-resonant magnetic', 'xray magnetic non-resonant',
                                     'non-resonant magnetic', 'magnetic non-resonant']:
                intensity += self.xray_resonant_magnetic(
                    _hkl,
                    self._energy_kev,
                    self._azimuthal_reference,
                    self._azimuthal_angle,
                    self._polarisation).tolist()
            else:
                print('Scattering type not defined')
        return np.array(intensity)

    def old_structure_factor(self, HKL, scattering_type=None):
        """
        Calculate the complex structure factor for the given HKL
          Crystal.structure_factor([1,0,0])
          Crystal.structure_factor([[1,0,0],[2,0,0],[3,0,0])
        Returns an array with the same length as HKL, giving the complex structure factor at each reflection.

        Notes:
        - Uses x-ray atomic form factors, calculated from approximated tables in the ITC
        - This may be a little slow for large numbers of reflections, as it is not currently
         possible to use accelerated calculation methods in Jython.
        - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs
        - Crystal.scale is used to scale the complex structure factor, so the intensity is
         reduced by (Crystal.scale)^2
        - Testing against structure factors calculated by Vesta.exe is very close, though there
          are some discrepancies, probably due to the method of calculation of the form factor.
        """
        prev_sf_setting = self._return_structure_factor
        self._return_structure_factor = True
        sf = self.intensity(HKL, scattering_type)
        self._return_structure_factor = prev_sf_setting
        return sf
    
    def hkl(self, HKL, energy_kev=None):
        """ Calculate the two-theta and intensity of the given HKL, display the result"""
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        tth = self.xtl.Cell.tth(HKL,energy_kev)
        inten = self.intensity(HKL)
        
        print('Energy = %6.3f keV' % energy_kev)
        print('( h, k, l) TwoTheta  Intensity')
        for n in range(len(tth)):
            print('(%2.0f,%2.0f,%2.0f) %8.2f  %9.2f' % (HKL[n,0],HKL[n,1],HKL[n,2],tth[n],inten[n]))
    
    def hkl_reflection(self, HKL, energy_kev=None):
        """
        Calculate the theta, two-theta and intensity of the given HKL in reflection geometry, display the result
        Uses sample orientation set up in setup_scatter
        :param HKL: [h,k,l] or list of hkl
        :param energy_kev: None or float
        :return: str
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        tth = self.xtl.Cell.tth(HKL,energy_kev)
        theta = self.xtl.Cell.theta_reflection(HKL, energy_kev, self._scattering_specular_direction, self._scattering_theta_offset)
        inten = self.intensity(HKL)
        
        print('Energy = %6.3f keV' % energy_kev)
        print('Specular Direction = (%1.0g,%1.0g,%1.0g)' %
              (self._scattering_specular_direction[0],
               self._scattering_specular_direction[1],
               self._scattering_specular_direction[2]))
        print('( h, k, l)    Theta TwoTheta  Intensity')
        for n in range(len(tth)):
            print('(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f' %
                  (HKL[n, 0], HKL[n, 1], HKL[n, 2], theta[n], tth[n], inten[n]))
    
    def hkl_transmission(self,HKL,energy_kev=None):
        " Calculate the theta, two-theta and intensity of the given HKL in transmission geometry, display the result"
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        tth = self.xtl.Cell.tth(HKL,energy_kev)
        theta = self.xtl.Cell.theta_transmission(HKL, energy_kev, self._scattering_specular_direction,self._scattering_theta_offset)
        inten = self.intensity(HKL)
        
        print('Energy = %6.3f keV' % energy_kev)
        print('Direction parallel to beam  = (%1.0g,%1.0g,%1.0g)' %(self._scattering_parallel_direction[0],self._scattering_parallel_direction[1],self._scattering_parallel_direction[2]))
        print('( h, k, l)    Theta TwoTheta  Intensity')
        for n in range(len(tth)):
            print('(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f' % (HKL[n,0],HKL[n,1],HKL[n,2],theta[n],tth[n],inten[n]))

    def generate_powder(self, q_max=8, peak_width=0.01, background=0, powder_average=True):
        """
        Generates array of intensities along a spaced grid, equivalent to a powder pattern.
          Q,I = generate_powder(energy_kev=8.0,peak_width=0.05,background=0)
            q_max = maximum Q, in A-1
            peak_width = width of convolution, in A-1
            background = average of normal background
            powder_average = True*/False, apply the powder averaging correction
          Returns:
            Q = [1000x1] array of wave-vector values
            I = [1000x1] array of intensity values

        Note: To get two-theta values use:
            tth = fc.cal2theta(Q, energy_kev)
        Note: To get d-spacing values use:
            dspace = fc.caldspace(Q)
        """

        # Get reflections
        hmax, kmax, lmax = fc.maxHKL(q_max, self.xtl.Cell.UVstar())
        HKL = fc.genHKL([hmax, -hmax], [kmax, -kmax], [lmax, -lmax])
        HKL = self.xtl.Cell.sort_hkl(HKL)  # required for labels
        Qmag = self.xtl.Cell.Qmag(HKL)
        HKL = HKL[Qmag < q_max, :]
        Qmag = self.xtl.Cell.Qmag(HKL)
        # Qmag = Qmag[Qmag<q_max]

        # Calculate intensities
        I = self.intensity(HKL)

        if powder_average:
            # Apply powder averging correction, I0/|Q|**2
            I = I/(Qmag+0.001)**2

        # create plotting mesh
        pixels = int(2000 * q_max)  # reduce this to make convolution faster
        pixel_size = q_max / (1.0 * pixels)
        peak_width_pixels = peak_width / (1.0 * pixel_size)
        mesh = np.zeros([pixels])
        mesh_q = np.linspace(0, q_max, pixels)

        # add reflections to background
        pixel_coord = Qmag / (1.0 * q_max)
        pixel_coord = (pixel_coord * (pixels - 1)).astype(int)

        for n in range(1, len(I)):
            mesh[pixel_coord[n]] = mesh[pixel_coord[n]] + I[n]

        # Convolve with a gaussian (if >0 or not None)
        if peak_width:
            gauss_x = np.arange(-3*peak_width_pixels, 3*peak_width_pixels + 1)  # gaussian width = 2*FWHM
            G = fg.gauss(gauss_x, None, height=1, cen=0, fwhm=peak_width_pixels, bkg=0)
            mesh = np.convolve(mesh, G, mode='same')

            # Add background (if >0 or not None)
        if background:
            bkg = np.random.normal(background, np.sqrt(background), [pixels])
            mesh = mesh + bkg
        return mesh_q, mesh

    def powder_correction(self, HKL, intensities, symmetric_multiplyer=True, powder_average=True):
        """
        Averages symmetric reflections and applies symmetry multipliers and 1/q^2 correction
        Ic = I0*C
        :param HKL: [nx3] array of [h,k,l] reflections
        :param intensities: [nx1] array of reflection intensities
        :return: [mx3], [mx1] arrays of averaged, corrected reflections + intensity
        """
        # Average symmetric reflections
        rhkl, rinten = self.xtl.Symmetry.average_symmetric_intensity(HKL, intensities)

        if symmetric_multiplyer:
            multiplyer = self.xtl.Symmetry.reflection_multiplyer(rhkl)
            rinten = rinten*multiplyer

        if powder_average:
            q = self.xtl.Cell.Qmag(rhkl)
            rinten = rinten/(q+0.001)**2

        return rhkl, rinten

    def print_all_reflections(self, energy_kev=None, print_symmetric=False,
                              min_intensity=0.01, max_intensity=None, units=None):
        """
        Prints a list of all allowed reflections at this energy
            energy_kev = energy in keV
            print_symmetric = False*/True : omits reflections with the same intensity at the same angle
            min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)
            max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)
            units = None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        if min_intensity is None: min_intensity = -1
        if max_intensity is None: max_intensity = np.inf
        
        hkl = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        if not print_symmetric:
            hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)
        hkl = self.xtl.Cell.sort_hkl(hkl)

        tth = self.xtl.Cell.tth(hkl, energy_kev)
        inrange = np.all([tth < self._scattering_max_twotheta, tth > self._scattering_min_twotheta], axis=0)
        hkl = hkl[inrange, :]
        tth = tth[inrange]
        #inten = np.sqrt(self.intensity(hkl)) # structure factor
        inten = self.intensity(hkl)

        if units is None:
            units = self._powder_units
        units = units.lower()
        if units in ['d', 'dspc', 'dspace', 'd space', 'd-space', 'dspacing', 'd spacing', 'd-spacing']:
            unit_str = 'd-spacing'
            unit = fc.caldspace(tth, energy_kev)
        elif units in ['q', 'wavevector']:
            unit_str = 'Q'
            unit = fc.calqmag(tth, energy_kev)
        else:
            unit_str = 'TwoTheta'
            unit = tth
        
        fmt = '(%3.0f,%3.0f,%3.0f) %10.2f  %9.2f\n'
        outstr = ''
        
        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Radiation: %s\n' % self._scattering_type
        outstr += '( h, k, l)    %10s  Intensity\n' % unit_str
        #outstr+= fmt % (hkl[0,0],hkl[0,1],hkl[0,2],tth[0],inten[0]) # hkl(0,0,0)
        count = 0
        for n in range(1, len(tth)):
            if inten[n] < min_intensity: continue
            if inten[n] > max_intensity: continue
            count += 1
            outstr += fmt % (hkl[n,0], hkl[n,1], hkl[n,2],unit[n],inten[n])
        outstr += 'Reflections: %1.0f\n' % count
        return outstr
    
    def print_ref_reflections(self, energy_kev=None, min_intensity=0.01, max_intensity=None):
        """
        Prints a list of all allowed reflections at this energy in the reflection geometry
            energy = energy in keV
            min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)
            max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)
                       |
                    // |\
                    \\ |
                     \\|___/____
                      \\   \
                       \\__
        
        Note, to change min/max theta values or specular direciton, change the following attributres of 
        the crystal object:
            self._scattering_max_two_theta   :  maximum detector (two-theta) angle
            self._scattering_min_theta       :  minimum sample angle = -opening angle
            self._scattering_max_theta       :  maximum sample angle = opening angle
            self._scattering_theta_offset    :  sample offset angle
            self._scattering_specular_direction : [h,k,l] : reflections normal to sample surface
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev

        if min_intensity is None: min_intensity = -1
        if max_intensity is None: max_intensity = np.inf

        HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        HKL = self.xtl.Cell.sort_hkl(HKL)
        tth = self.xtl.Cell.tth(HKL, energy_kev)
        HKL = HKL[tth > self._scattering_min_twotheta, :]
        tth = tth[tth > self._scattering_min_twotheta]
        theta = self.xtl.Cell.theta_reflection(HKL, energy_kev, self._scattering_specular_direction,
                                               self._scattering_theta_offset)
        # inten = np.sqrt(self.intensity(HKL)) # structure factor
        inten = self.intensity(HKL)

        p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)
        p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))
        pos_theta = p1 * p2
        
        fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'
        outstr = ''

        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Radiation: %s\n' % self._scattering_type
        outstr += 'Specular Direction = (%1.0g,%1.0g,%1.0g)\n' % (
            self._scattering_specular_direction[0],
            self._scattering_specular_direction[1],
            self._scattering_specular_direction[2]
        )
        outstr += '( h, k, l) TwoTheta    Theta  Intensity\n'
        # outstr+= fmt % (HKL[0,0],HKL[0,1],HKL[0,2],tth[0],theta[0],inten[0]) # hkl(0,0,0)
        count = 0
        for n in range(1, len(tth)):
            if inten[n] < min_intensity: continue
            if inten[n] > max_intensity: continue
            if not pos_theta[n]: continue
            # if not print_symmetric and np.abs(tth[n]-tth[n-1]) < 0.01: continue # only works if sorted
            count += 1
            outstr += fmt % (HKL[n, 0], HKL[n, 1], HKL[n, 2], tth[n], theta[n], inten[n])
        outstr += 'Reflections: %1.0f\n' % count
        return outstr
    
    def print_tran_reflections(self, energy_kev=None, min_intensity=0.01, max_intensity=None):
        """
        Prints a list of all allowed reflections at this energy in the transmission geometry
            energy = energy in keV
            min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)
            max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)
                   \ /      
             --<-- || --<-- 
                  / \       
        
        Note, to change min/max theta values or specular direciton, change the following attributres of 
        the crystal object:
            self._scattering_max_two_theta   :  maximum detector (two-theta) angle
            self._scattering_min_theta       :  minimum sample angle = -opening angle
            self._scattering_max_theta       :  maximum sample angle = opening angle
            self._scattering_theta_offset    :  sample offset angle
            self._scattering_parallel_direction : [h,k,l] : reflections normal to sample surface, parallel to beam at theta = 0
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev

        if min_intensity is None: min_intensity = -1
        if max_intensity is None: max_intensity = np.inf

        HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        HKL = self.xtl.Cell.sort_hkl(HKL)
        tth = self.xtl.Cell.tth(HKL, energy_kev)
        HKL = HKL[tth > self._scattering_min_twotheta, :]
        tth = tth[tth > self._scattering_min_twotheta]
        theta = self.xtl.Cell.theta_transmission(HKL, energy_kev, self._scattering_parallel_direction)
        # inten = np.sqrt(self.intensity(HKL)) # structure factor
        inten = self.intensity(HKL)

        p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)
        p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))
        pos_theta = p1 * p2

        fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'
        outstr = ''

        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Radiation: %s\n' % self._scattering_type
        outstr += 'Direction parallel to beam  = (%1.0g,%1.0g,%1.0g)\n' % (
            self._scattering_parallel_direction[0],
            self._scattering_parallel_direction[1],
            self._scattering_parallel_direction[2]
        )
        outstr += '( h, k, l) TwoTheta    Theta  Intensity\n'
        # outstr+= fmt % (HKL[0,0],HKL[0,1],HKL[0,2],tth[0],theta[0],inten[0]) # hkl(0,0,0)
        count = 0
        for n in range(1, len(tth)):
            if inten[n] < min_intensity: continue
            if inten[n] > max_intensity: continue
            if not pos_theta[n]: continue
            # if not print_symmetric and np.abs(tth[n]-tth[n-1]) < 0.01: continue # only works if sorted
            count += 1
            outstr += fmt % (HKL[n, 0], HKL[n, 1], HKL[n, 2], tth[n], theta[n], inten[n])
        outstr += ('Reflections: %1.0f\n' % count)
        return outstr
    
    def print_symmetric_reflections(self,HKL):
        "Prints equivalent reflections"
        
        symHKL = self.xtl.Symmetry.symmetric_reflections(HKL)
        Nsyms = len(symHKL)
        outstr = ''
        outstr+= 'Equivalent reflections: %d\n' % Nsyms
        for n in range(Nsyms):
            outstr+= '(%5.3g,%5.3g,%5.3g)\n' % (symHKL[n,0],symHKL[n,1],symHKL[n,2])
        return outstr
    
    def print_atomic_contributions(self,HKL):
        """
        Prints the atomic contributions to the structure factor
        """
        
        HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])
        Nref = len(HKL)
        
        # Calculate the full intensity
        I = self.intensity(HKL)
        
        # Calculate the structure factors of the symmetric atomic sites
        base_label = self.xtl.Atoms.label
        uvw,type,label,occ,uiso,mxmymz = self.xtl.Structure.get()
        
        Qmag = self.xtl.Cell.Qmag(HKL)
        
        # Get atomic form factors
        ff = fc.xray_scattering_factor(type,Qmag)
        
        # Get Debye-Waller factor
        dw = fc.debyewaller(uiso,Qmag)
        
        # Calculate dot product
        dot_KR = np.dot(HKL,uvw.T)
        
        # Calculate structure factor
        SF =  ff*dw*occ*np.exp(1j*2*np.pi*dot_KR)
        
        # Sum structure factors of each base label in atoms
        SFbase = np.zeros([len(HKL),len(base_label)],dtype=np.complex128)
        for n in range(len(base_label)):
            label_idx = label == base_label[n]
            SFbase[:,n] = np.sum(SF[:,label_idx],axis=1)
        
        # Get the real part of the structure factor
        #SFtot = np.sqrt(np.real(SF * np.conj(SF)))
        SFrel = np.real(SFbase)
        SFimg = np.imag(SFbase)
        
        # Generate the results
        outstr = ''
        outstr+= '( h, k, l) Intensity' + ' '.join(['%12s    '%x for x in base_label])+'\n'
        for n in range(Nref):
            ss = ' '.join(['%6.1f + i%-6.1f' % (x,y) for x,y in zip(SFrel[n],SFimg[n])])
            outstr+= '(%2.0f,%2.0f,%2.0f) %9.2f    %s\n' % (HKL[n,0],HKL[n,1],HKL[n,2],I[n],ss)
        return outstr

    def print_symmetry_contributions(self,HKL):
        """
        Prints the symmetry contributions to the structure factor for each atomic site
        """

        HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])
        Nref = len(HKL)

        # Calculate the full intensity
        I = self.intensity(HKL)

        # Calculate the structure factors of the symmetric atomic sites
        buvw, btype, base_label, bocc, buiso, bmxmymz = self.xtl.Atoms.get()
        operations = np.hstack([self.xtl.Symmetry.symmetric_coordinate_operations(buvw[n]) for n in range(len(buvw))])
        rotations = np.hstack(
            [self.xtl.Symmetry.symmetric_coordinate_operations(buvw[n], True)[1] for n in range(len(buvw))])

        # Calculate the structure factors
        uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
        Qmag = self.xtl.Cell.Qmag(HKL)
        ff = fc.xray_scattering_factor(type, Qmag)
        dw = fc.debyewaller(uiso, Qmag)
        dot_KR = np.dot(HKL, uvw.T)
        phase = np.exp(1j * 2 * np.pi * dot_KR)
        sf = ff * dw * occ * phase

        # Generate the results
        outstr = ''
        for n in range(Nref):
            ss = '\n'
            all_phase = 0j
            all_sf = 0j
            for lab in base_label:
                label_idx = np.argwhere(label == lab)
                ss += '  %s\n' % lab
                tot_phase = 0j
                tot_sf = 0j
                for a in label_idx:
                    uvwstr = '(%-7.3g,%-7.3g,%-7.3g)' % (uvw[a, 0], uvw[a, 1], uvw[a, 2])
                    phstr = fg.complex2str(phase[n, a])
                    sfstr = fg.complex2str(sf[n, a])
                    val = (a, uvwstr, operations[a[0]], rotations[a[0]], phstr, sfstr)
                    ss += '    %3d %s %25s %20s  %s  %s\n' % val
                    tot_phase += phase[n, a]
                    tot_sf += sf[n, a]
                ss += '%74sTotal:  %s  %s\n' % (' ', fg.complex2str(tot_phase), fg.complex2str(tot_sf))
                all_phase += tot_phase
                all_sf += tot_sf
            ss += '%62s Reflection Total:  %s  %s\n' % (' ', fg.complex2str(all_phase), fg.complex2str(all_sf))
            outstr += '(%2.0f,%2.0f,%2.0f) I = %9.2f    %s\n' % (HKL[n, 0], HKL[n, 1], HKL[n, 2], I[n], ss)
        return outstr
    
    def find_close_reflections(self,HKL,energy_kev=None,max_twotheta=2,max_angle=10):
        """
        Find and print list of reflections close to the given one
        """
        
        if energy_kev is None:
            energy_kev = self._energy_kev
        
        all_HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        all_HKL = self.xtl.Cell.sort_hkl(all_HKL)
        
        all_tth = self.xtl.Cell.tth(all_HKL,energy_kev)
        tth1 = self.xtl.Cell.tth(HKL, energy_kev)
        tth_dif = np.abs(all_tth-tth1)
        
        all_Q = self.xtl.Cell.calculateQ(all_HKL)
        Q1 = self.xtl.Cell.calculateQ(HKL)
        all_angles = np.abs([fg.ang(Q1,Q2,'deg') for Q2 in all_Q])
        
        selected = (tth_dif < max_twotheta)*(all_angles < max_angle)
        sel_HKL = all_HKL[selected,:]
        sel_tth = all_tth[selected]
        sel_angles = all_angles[selected]
        sel_intensity = self.intensity(sel_HKL)
        
        # Generate Results
        fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'
        outstr = '' 
        
        outstr+= 'Energy = %6.3f keV\n' % energy_kev
        outstr+= 'Close to Reflection (%2.0f,%2.0f,%2.0f)\n' %(HKL[0],HKL[1],HKL[2])
        outstr+= '( h, k, l) TwoTheta    Angle  Intensity\n'
        count = 0
        for n in range(0,len(sel_HKL)):
            count += 1
            outstr+= fmt % (sel_HKL[n,0],sel_HKL[n,1],sel_HKL[n,2],sel_tth[n],sel_angles[n],sel_intensity[n])
        outstr+= 'Reflections: %1.0f\n'%count
        return outstr

    def multiple_scattering(self, hkl, azir=[0, 0, 1], pv=[1, 0], energy_range=[7.8, 8.2], numsteps=60,
                            full=False, pv1=False, pv2=False, sfonly=True, pv1xsf1=False):
        """
        Run multiple scattering code, plot result.

        mslist = xtl.Scatter.multiple_scattering([h,k,l], energy_range=[7.8, 8.2])

        See multiple_scattering.py for more details. Code created by Dr <NAME>, DLS
        :param hkl: [h,k,l] principle reflection
        :param azir: [h,k,l] reference of azimuthal 0 angle
        :param pv: [s,p] polarisation vector
        :param energy_range: [min, max] energy range in keV
        :param numsteps: int: number of calculation steps from energy min to max
        :param full: True/False: calculation type: full
        :param pv1: True/False: calculation type: pv1
        :param pv2: True/False: calculation type: pv2
        :param sfonly: True/False: calculation type: sfonly *default
        :param pv1xsf1: True/False: calculation type: pv1xsf1?
        :return: array
        """
        return ms.run_calcms(self.xtl, hkl, azir, pv, energy_range, numsteps,
                             full=full, pv1=pv1, pv2=pv2, sfonly=sfonly, pv1xsf1=pv1xsf1)

    def ms_azimuth(self, hkl, energy_kev, azir=[0, 0, 1], pv=[1, 0], numsteps=3, peak_width=0.1,
                   full=False, pv1=False, pv2=False, sfonly=True, pv1xsf1=False):
        """
        Returns an azimuthal dependence at a particular energy

        :param xtl: Crystal structure from Dans_Diffraction
        :param hkl: [h,k,l] principle reflection
        :param energy_kev: calculation energy
        :param azir: [h,k,l] reference of azimuthal 0 angle
        :param pv: [s,p] polarisation vector
        :param numsteps: int: number of calculation steps from energy min to max
        :param peak_width: convolution width
        :param full: True/False: calculation type: full
        :param pv1: True/False: calculation type: pv1
        :param pv2: True/False: calculation type: pv2
        :param sfonly: True/False: calculation type: sfonly *default
        :param pv1xsf1: True/False: calculation type: pv1xsf1?
        :return: None
        """

        energy_range = [energy_kev-0.001, energy_kev+0.001]
        mslist = self.xtl.Scatter.multiple_scattering(hkl, azir, pv, energy_range, numsteps,
                                                      full=full, pv1=pv1, pv2=pv2, sfonly=sfonly, pv1xsf1=pv1xsf1)

        if pv1 + pv2 + sfonly + full + pv1xsf1 != 0:
            azimuth = np.concatenate(mslist[:, [3, 4]])
            intensity = np.concatenate(mslist[:, [-1, -1]])
        else:
            azimuth = np.concatenate(mslist[:, [3, 4]])
            intensity = np.ones(azimuth)

        # create plotting mesh
        peak_width_pixels = 10
        pixels = int(360*(peak_width_pixels/ peak_width))
        mesh_azi = np.linspace(-180, 180, pixels)
        mesh = np.zeros(mesh_azi.shape)

        # add reflections to background
        pixel_coord = (azimuth-180) / 360.
        pixel_coord = (pixel_coord * pixels).astype(int)

        for n in range(len(intensity)):
            mesh[pixel_coord[n]] = mesh[pixel_coord[n]] + intensity[n]

        # Convolve with a gaussian (if >0 or not None)
        if peak_width:
            gauss_x = np.arange(-3 * peak_width_pixels, 3 * peak_width_pixels + 1)  # gaussian width = 2*FWHM
            G = fg.gauss(gauss_x, None, height=1, cen=0, fwhm=peak_width_pixels, bkg=0)
            mesh = np.convolve(mesh, G, mode='same')
        return mesh_azi, mesh

    '''  Removed tensor scattering 26/05/20 V1.7
    def tensor_scattering(self, atom_label, hkl, energy_kev=None, azir=[0, 0, 1], psideg=0, process='E1E1',
                          rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None):
        """
        Return tensor scattering intensities
          ss, sp, ps, pp = tensor_scattering('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param hkl: list/array, [h,k,l] reflection to calculate
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param psideg: float/array, azimuthal angle/ range
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :return: ss, sp, ps, pp intensity values
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        # Calculate tensor scattering
        ss, sp, ps, pp = ts.CalculateIntensityInPolarizationChannels(
            process, B, sitevec, sglist, lam, hkl, azir, psideg,
            K=rank, Time=time, Parity=parity, mk=mk, lk=lk, sk=sk
        )
        return ss, sp, ps, pp

    def tensor_scattering_stokes(self, atom_label, hkl, energy_kev=None, azir=[0, 0, 1], psideg=0, stokes=0,
                                 pol_theta=45, process='E1E1', rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None):
        """
        Return tensor scattering intensities for non-standard polarisation
          pol = tensor_scattering_stokes('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90, stokes=45)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param hkl: list/array, [h,k,l] reflection to calculate
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param psideg: float, azimuthal angle
        :param stokes: float/array, rotation of polarisation analyser (0=sigma), degrees
        :param pol_theta: float, scattering angle of polarisation analyser, degrees
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :return: array of intensity values
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        # Calculate tensor scattering
        stokesvec_swl = [0, 0, 1]
        pol = ts.CalculateIntensityFromPolarizationAnalyser(
            process, B, sitevec, sglist, lam, hkl, azir, psideg, stokes, pol_theta,
            stokesvec_swl, rank, time, parity, mk, lk, sk
        )
        return pol

    def print_tensor_scattering(self, atom_label, hkl, energy_kev=None, azir=[0, 0, 1], psideg=0, process='E1E1',
                                rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None):
        """
        Return tensor scattering intensities
          ss, sp, ps, pp = tensor_scattering('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param hkl: list/array, [h,k,l] reflection to calculate
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param psideg: float, azimuthal angle
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :return: str
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        # Calculate tensor scattering
        ss, sp, ps, pp = ts.CalculateIntensityInPolarizationChannels(
            process, B, sitevec, sglist, lam, hkl, azir, psideg,
            K=rank, Time=time, Parity=parity, mk=mk, lk=lk, sk=sk
        )
        outstr1 = ts.tensorproperties(sitevec, sglist, hkl, Parity=parity, Time=time)
        outstr2 = ts.print_tensors(B, sitevec, sglist, hkl, K=rank, Parity=parity, Time=time)
        outstr3 = "\nScattering Tensor:\n\n    [ss, sp] = [%5.2f, %5.2f]\n    [ps, pp]   [%5.2f, %5.2f]"
        outstr3 = outstr3 % (ss, sp, ps, pp)
        return outstr1 + outstr2 + outstr3

    def print_tensor_scattering_refs(self, atom_label, energy_kev=None, azir=[0, 0, 1], psideg=0, process='E1E1',
                                     rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None,
                                     print_symmetric=False, units=None):
        """
        Return tensor scattering intensities for all reflections at given azimuth and energy
          ss, sp, ps, pp = tensor_scattering('Ru1', 2.838, [0,1,0], psideg=90)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param psideg: float, azimuthal angle
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :param print_symmetric: False*/True : omits reflections with the same intensity at the same angle
        :param units: None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display
        :return: str
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        hkl = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        if not print_symmetric:
            hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)
        hkl = self.xtl.Cell.sort_hkl(hkl)

        tth = self.xtl.Cell.tth(hkl, energy_kev)
        inrange = np.all([tth < self._scattering_max_twotheta, tth > self._scattering_min_twotheta], axis=0)
        hkl = hkl[inrange, :]
        tth = tth[inrange]

        if units is None:
            units = self._powder_units
        units = units.lower()
        if units in ['d', 'dspc', 'dspace', 'd space', 'd-space', 'dspacing', 'd spacing', 'd-spacing']:
            unit_str = 'd-spacing'
            unit = fc.caldspace(tth, energy_kev)
        elif units in ['q', 'wavevector']:
            unit_str = 'Q'
            unit = fc.calqmag(tth, energy_kev)
        else:
            unit_str = 'TwoTheta'
            unit = tth

        fmt = '(%3.0f,%3.0f,%3.0f) %10.2f  %11.2f %11.2f %11.2f %11.2f\n'
        outstr = 'Tensor Scattering %s\n' % self.xtl.name
        outstr += 'Process: %s, site: %s\n' % (process, atom_label)
        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Psi_0 = (%.2g,%.2g,%.2g)  Psi = %3.3g\n' % (azir[0], azir[1], azir[2], psideg)
        outstr += '( h, k, l)    %10s  Sigma-Sigma    Sigma-Pi    Pi-Sigma       Pi-Pi\n' % unit_str

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        for n in range(1, len(tth)):
            ss, sp, ps, pp = ts.CalculateIntensityInPolarizationChannels(
                process, B, sitevec, sglist, lam, hkl[n, :], azir, psideg,
                K=rank, Time=time, Parity=parity, mk=mk, lk=lk, sk=sk
            )
            outstr += fmt % (hkl[n, 0], hkl[n, 1], hkl[n, 2], unit[n], ss, sp, ps, pp)
        outstr += 'Reflections: %1.0f\n' % len(tth)
        return outstr

    def print_tensor_scattering_refs_max(self, atom_label, energy_kev=None, azir=[0, 0, 1], process='E1E1',
                                     rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None,
                                     print_symmetric=False, units=None):
        """
        Return tensor scattering intensities for all reflections at given energy at maximum intensity psi
          ss, sp, ps, pp = tensor_scattering('Ru1', 2.838, [0,1,0], psideg=90)
        Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.
        :param atom_label: str atom site label, e.g. Ru1
        :param energy_kev: float
        :param azir: list/array, [h,k,l] azimuthal reference
        :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'
        :param rank: int, 1,2,3: tensor rank. Only required
        :param time: +/-1 time symmetry
        :param parity: +/-1 parity
        :param mk:
        :param lk:
        :param sk:
        :param print_symmetric: False*/True : omits reflections with the same intensity at the same angle
        :param units: None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display
        :return: str
        """

        if energy_kev is None:
            energy_kev = self._energy_kev

        hkl = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)
        if not print_symmetric:
            hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)
        hkl = self.xtl.Cell.sort_hkl(hkl)

        tth = self.xtl.Cell.tth(hkl, energy_kev)
        inrange = np.all([tth < self._scattering_max_twotheta, tth > self._scattering_min_twotheta], axis=0)
        hkl = hkl[inrange, :]
        tth = tth[inrange]
        # Caluclate structure factor **2
        inten = self.intensity(hkl)

        if units is None:
            units = self._powder_units
        units = units.lower()
        if units in ['d', 'dspc', 'dspace', 'd space', 'd-space', 'dspacing', 'd spacing', 'd-spacing']:
            unit_str = 'd-spacing'
            unit = fc.caldspace(tth, energy_kev)
        elif units in ['q', 'wavevector']:
            unit_str = 'Q'
            unit = fc.calqmag(tth, energy_kev)
        else:
            unit_str = 'TwoTheta'
            unit = tth

        fmt = '(%3.0f,%3.0f,%3.0f) %10.2f %10.2f  %5.2f (%3.0f) %5.2f (%3.0f) %5.2f (%3.0f) %5.2f (%3.0f)\n'
        outstr = 'Tensor Scattering %s\n' % self.xtl.name
        outstr += 'Process: %s, site: %s\n' % (process, atom_label)
        outstr += 'Energy = %6.3f keV\n' % energy_kev
        outstr += 'Psi_0 = (%.2g,%.2g,%.2g)\n' % (azir[0], azir[1], azir[2])
        outstr += '( h, k, l)    %10s         I0  Sigma-Sigma    Sigma-Pi    Pi-Sigma       Pi-Pi\n' % unit_str

        # Structure data for input
        B = ts.latt2b(self.xtl.Cell.lp())
        sitevec = self.xtl.Atoms.uvw()[self.xtl.Atoms.label.index(atom_label)]
        sglist = ts.spacegroup_list_from_genpos_list(self.xtl.Symmetry.symmetry_operations)
        lam = fc.energy2wave(energy_kev)

        psideg = np.arange(0, 361, 5)
        for n in range(1, len(tth)):
            ss, sp, ps, pp = ts.CalculateIntensityInPolarizationChannels(
                process, B, sitevec, sglist, lam, hkl[n, :], azir, psideg,
                K=rank, Time=time, Parity=parity, mk=mk, lk=lk, sk=sk
            )
            iss = np.argmax(ss)
            isp = np.argmax(sp)
            ips = np.argmax(ps)
            ipp = np.argmax(pp)
            ss = ss[iss]
            sp = sp[isp]
            ps = ps[ips]
            pp = pp[ipp]
            psi_ss = psideg[iss]
            psi_sp = psideg[isp]
            psi_ps = psideg[ips]
            psi_pp = psideg[ipp]
            outstr += fmt % (hkl[n, 0], hkl[n, 1], hkl[n, 2], unit[n], inten[n],
                             ss, psi_ss, sp, psi_sp, ps, psi_ps, pp, psi_pp)
        outstr += 'Reflections: %1.0f\n' % len(tth)
        return outstr
    '''


class ScatteringTypes:
    """
    Container for available scattering types
    """
    def __init__(self, parent, typedict):
        self.parent = parent
        self.typedict = typedict

        typenames = typedict.keys()

        for name in typenames:
            attrname = name.replace(' ','_')
            setattr(self, attrname, self.ScatteringType(parent, name))

    def __call__(self, type):
        """
        Calling this container with a scattering type will set the scattering type in the parent class
        :param type: str scattering type
        :return: None
        """
        for typename in self.typedict:
            if type.lower() in self.typedict[typename]:
                self.parent._scattering_type = typename
                print('Changed scattering type to: %s'%typename)

    def i16(self):
        """
        Set max/min angles consistent with reflection geometry on beamline I16
        :return: None
        """
        self.parent.setup_scatter(type='xray', energy_kev=8.0, min_theta=-20, max_theta=140, min_twotheta=0, max_twotheta=140)

    def wish(self):
        """
        Set max/min angles consistent with reflection geometry on beamline I16
        :return: None
        """
        self.parent.setup_scatter(type='neutron', wavelength_a=1.5, min_theta=-180, max_theta=180, min_twotheta=-180, max_twotheta=180)

    class ScatteringType:
        """
        Container for scattering type switcher
        """
        def __init__(self, parent, typename):
            self.parent = parent
            self.typename = typename

        def __call__(self):
            self.parent._scattering_type = self.typename

# -*- coding: utf-8 -*- Scattering Class "classes_scattering.py"     Scattering functions for Crystal class.  By <NAME>, PhD Diamond 2017  Version 2.0.1 Last updated: 28/09/21  Version History: 10/09/17 0.1    Program created 30/10/17 1.0    Main functions finshed, some testing done 06/01/18 1.1    Renamed classes_scattering.py 31/10/18 1.2    Added print_symmetric_contributions 21/01/19 1.3    Added non-resonant diffraction, corrected resonant diffraction 16/12/19 1.4    Added multiple_scattering code, print_all_reflections updated with units 18/02/20 1.5    Added tensor_scattering code 20/03/20 1.6    Increased powder gauss width from 2fwhm to 6fwhm, added powder averaging 14/04/20 1.6    Added powder_correction 26/05/20 1.7    Removed tensor_scattering 16/06/20 1.7.1  Added output option of setup_scatter 04/01/21 1.7.2  Added structure_factor function 21/01/21 1.8    Added xray_dispersion scattering function 10/06/21 1.9    Added x_ray calculation using Waasmaier and Kirfel scattering factors. 09/07/21 1.9    Added new scattering factors as option on normal scattering functions 20/08/21 2.0    Switched over to new scattering module, added self.powder() 28/09/21 2.0.1  Added __repr__  @author: DGPorter # from . import tensor_scattering as ts  # Removed V1.7 Simulate diffraction from Crystal     Useage:         xtl = Crystal()         xtl.Scatter.setup_scatter(type='x-ray',energy_keV=8.0)         xtl.Scatter.intensity([h,k,l]) # Returns intensity         print(xtl.Scatter.print_all_refelctions()) # Returns formated string of all allowed reflections                  Allowed radiation types:             'xray','neutron','xray magnetic','neutron magnetic','xray resonant' #------Options------- # Standard Options # added so recalculation not required # 'xray','neutron','xray magnetic','neutron magnetic','xray resonant' # reflection # transmission # powder options # tth (two theta), Q, d # in Deg # Complex Structure factor # Uses the coefficients for analytical approximation to the scattering factors from: #        "<NAME>, Acta Cryst. (1995) A51, 416-431" # Thermal Factors # Magnetic Options # Polarisation Options # Radiation energy # Resonant X-ray Options # Initialise the scattering type container #out += '  Anisotropic ADPs: %s\n' % self._use_anisotropic_thermal_factor #out += '  use e2e2 approx.: %s\n' % self._resonant_approximation_e2e2 #out += '  use e1e2 approx.: %s\n' % self._resonant_approximation_e1e2 #out += '  use m1m1 approx.: %s\n' % self._resonant_approximation_m1m1 Simple way to set scattering parameters, each parameter is internal to xtl (self)          scattering_type: self._scattering type            :  'xray','neutron','xray magnetic','neutron magnetic','xray resonant', 'xray dispersion'         energy_kev  : self._energy_kev                    :  radiation energy in keV         wavelength_a: self._wavelength_a                  :  radiation wavelength in Angstrom         powder_units: self._powder_units                  :  units to use when displaying/ plotting ['twotheta', 'd',' 'q']         min_twotheta: self._scattering_min_two_theta      :  minimum detector (two-theta) angle         max_twotheta: self._scattering_max_two_theta      :  maximum detector (two-theta) angle         min_theta   : self._scattering_min_theta          :  minimum sample angle = -opening angle         max_theta   : self._scattering_max_theta          :  maximum sample angle = opening angle         theta_offset: self._scattering_theta_offset       :  sample offset angle         specular    : self._scattering_specular_direction : [h,k,l] : reflections normal to sample surface         parallel    : self._scattering_parallel_direction : [h,k,l] : reflections normal to sample surface         scattering_factors: self._use_waaskirf_scattering_factor : xray scattering factor ['waaskirf', 'itc']         magnetic_formfactor: self._use_magnetic_form_factor: True/False magnetic form factor for magnetic SF         polarisation: self._polarisation                  : beam polarisation setting ['ss', 'sp'*, 'sp', 'pp']         polarisation_vector: _polarisation_vector_incident: [x,y,z] incident polarisation vector         azimuthal_reference: self._azimuthal_reference    : [h,k,l] direction of azimuthal zero angle         azimuth    : self._azimuthal_angle                : azimuthal angle in deg         flm        : self._resonant_flm                   : Resonant settings (flm1, flm2, flm3) # print('Scattering Options:') # print('                            Type : %s' % self._scattering_type) # print('                  Default Energy : %6.3f keV' % self._energy_kev) # print('                    Powder Units : %s' % self._powder_units) # print('  Specular Direction (reflection): (%2.0f,%2.0f,%2.0f)' % ( #     self._scattering_specular_direction[0], self._scattering_specular_direction[1], #     self._scattering_specular_direction[2])) # print('Parallel Direction (transmission): (%2.0f,%2.0f,%2.0f)' % ( #     self._scattering_parallel_direction[0], self._scattering_parallel_direction[1], #     self._scattering_parallel_direction[2])) # print('                   Sample Offset : %5.2f' % self._scattering_theta_offset) # print('             Minimum Theta angle : %5.2f' % self._scattering_min_theta) # print('             Maximum Theta angle : %5.2f' % self._scattering_max_theta) # print('         Minimum Two-Theta angle : %5.2f' % self._scattering_min_twotheta) # print('         Maximum Two-Theta angle : %5.2f' % self._scattering_max_twotheta) Return stored hkl or generate         :param regenerate: if True, hkl list will be regenerated, if False - previous list will be returned         :param remove_symmetric: generate only non-symmetric hkl values         :param reflection: generate only reflections possible in reflection geometry         :param transmission: generate only reflections possible in transmission geometry         :param kwargs: additional options to pass to setup_scatter()         :return: array # remove (0, 0, 0) Toggle on debug mode Calculate the structure factor at reflection indexes (h,k,l)          Notes:         - Uses x-ray atomic form factors, calculated from approximated tables in the ITC         - This may be a little slow for large numbers of reflections, as it is not currently          possible to use accelerated calculation methods in Jython.         - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs         - Crystal.scale is used to scale the complex structure factor, so the intensity is          reduced by (Crystal.scale)^2         - Testing against structure factors calculated by Vesta.exe is exactly the same when using Waasmaier structure factors.          :param hkl: array[n,3] : reflection indexes (h, k, l)         :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']         :param kwargs: additional options to pass to scattering function         :return: complex array[n] : structure factors # print('Scattering function: %s' % scattering_fun.__name__) # Break up long lists of HKLs # Scattering factors # Get Debye-Waller factor # Get magnetic form factors # shape(nref) # shape(nref, nenergy) # shape(nref, nspi) :param hkl:         :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']         :return: Generates array of intensities along a spaced grid, equivalent to a powder pattern.           tth, iten = generate_powder('xray', units='tth', energy_kev=8)          :param scattering_type: str : one of ['xray','neutron','xray magnetic','neutron magnetic','xray resonant']         :param units: str : one of ['tth', 'dspace', 'q']         :param peak_width: float : Peak with in units of inverse wavevector (Q)         :param background: float : if >0, a normal background around this value will be added         :param pixels: int : number of pixels to add to the resulting mesh         :param powder_average: Bool : if True, intensities will be reduced for the powder average         :param options: additional arguments to pass to intensity calculation         :return xval: arrray : x-axis of powder scan (units)         :return yval: array :  intensity values at each point in x-axis # Units # Get reflections # required for labels # Calculate intensities # Apply powder averging correction, I0/|Q|**2 # create plotting mesh # reduce this to make convolution faster # add reflections to background # Convolve with a gaussian (if >0 or not None) # gaussian width = 2*FWHM # Add background (if >0 or not None) # Change output units Calculate the squared structure factor for the given HKL, using x-ray scattering factors           Scattering.x_ray([1,0,0])           Scattering.x_ray([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the real intensity at each reflection. # Get atomic form factors # Get Debye-Waller factor # Calculate dot product # Calculate structure factor # Broadcasting used on 2D ff #SF = np.zeros(Nref,dtype=np.complex) #for ref in range(Nref): #    for at in range(Nat): #        SF[ref] += ff[ref,at]*dw[ref,at]*occ[at]*np.exp(1j*2*np.pi*dot_KR[ref,at]) # Calculate intensity Calculate the squared structure factor for the given HKL, using atomic number as scattering length           Scattering.x_ray_fast([1,0,0])           Scattering.x_ray_fast([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the real intensity at each reflection. # Get atomic form factors # Get Debye-Waller factor # Calculate dot product # Calculate structure factor # Calculate intensity Calculate the squared structure factor for the given HKL,         using x-ray scattering factors with dispersion corrections           Scattering.xray_dispersion([1,0,0], 9.71)           Scattering.xray_dispersion([[1,0,0],[2,0,0],[3,0,0]], 2.838)           Scattering.xray_dispersion([[1,0,0],[2,0,0],[3,0,0]], np.arange(2.83, 2.86, 0.001))         Returns an array with the same length as HKL, giving the real intensity at each reflection.         if energy_kev is an array, the returned array will have shape [len(HKL), len(energy_kev)] # Get atomic form factors # shape (len(HKL), len(type), len(en)) # Get Debye-Waller factor # Calculate dot product # Broadcasting used on 2D ff # Calculate intensity Calculate the squared structure factor for the given HKL, using neutron scattering length           Scattering.neutron([1,0,0])           Scattering.neutron([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the real intensity at each reflection. # Get atomic form factors # Get Debye-Waller factor # Calculate dot product # Calculate structure factor # Calculate intensity Calculate the magnetic component of the structure factor for the given HKL, using neutron rules and form factor           Scattering.magnetic_neutron([1,0,0])           Scattering.magnetic_neutron([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the real intensity at each reflection. # Get magnetic form factors # Calculate moment # broadcast n*1 x n*3 = n*3 # Calculate dot product # Calculate structure factor # Calculate Magnetic part # Calculate structure factor # Calculate polarisation with incident neutron #SF[n] = np.dot(SFm,SFm) # maximum possible # average polarisation # Calculate intensity Calculate the non-resonant magnetic component of the structure factor          for the given HKL, using x-ray rules and form factor           Scattering.xray_magnetic([1,0,0])           Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the real intensity at each reflection.                  From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)         Book: "X-ray Scattering and Absorption by Magnetic Materials" by <NAME> Collins. Ch 2. Eqn.2.21+1         No orbital component assumed         magnetic moments assumed to be in the same reference frame as the polarisation # Get magnetic form factors # Calculate moment # moment direction in cartesian reference frame # broadcast n*1 x n*3 = n*3 # Calculate dot product # Calculate structure factor # Calculate vector structure factor # Calculate polarisation with incident x-ray # The reference frame of the x-ray and the crystal are assumed to be the same # i.e. pol=[1,0,0] || mom=[1,0,0] || (1,0,0) #SF[n] = np.dot(SFm,SFm) # maximum possible # average polarisation # Calculate intensity Calculate structure factors using resonant scattering factors in the dipolar approximation           I = Scattering.xray_resonant(HKL,energy_kev,polarisation,F0,F1,F2)         Returns an array with the same length as HKL, giving the real intensity at each reflection.             energy_kev = x-ray energy in keV             polarisation = x-ray polarisation: 'ss',{'sp'},'ps','pp'             F0/1/2 = Resonance factor Flm             azim_zero = [h,k,l] vector parallel to the origin of the azimuth             psi = azimuthal angle defining the scattering plane                  Uses the E1E1 resonant x-ray scattering amplitude:             fxr_n = (ef.ei)*F0 -i(ef X ei).z_n*F1 + (ef.z_n)(ei.z_n)F2                  Where ei and ef are the initial and final polarisation states, respectively,         and z_n is a unit vector in the direction of the magnetic moment of the nth ion.         The polarisation states are determined to be one of the natural synchrotron          states, where sigma (s) is perpendicular to the scattering plane and pi (p) is         parallel to it.                 ( s-s  s-p )                 ( p-s  p-p )                  From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (15) # Get Debye-Waller factor # Calculate dot product # Get resonant form factor # Calculate structure factor # Broadcasting used on 2D fxres # Calculate intensity Calcualte fxres, the resonant x-ray scattering factor           fxres = Scattering.xray_resonant_scattering_factor(HKL,energy_kev,polarisation,F0,F1,F2,azim_zero,psi)         energy_kev = x-ray energy in keV             polarisation = x-ray polarisation: 'ss',{'sp'},'ps','pp'             F0/1/2 = Resonance factor Flm             azim_zero = [h,k,l] vector parallel to the origin of the azimuth {[1,0,0]}             psi = azimuthal angle defining the scattering plane {0}                  Uses the E1E1 resonant x-ray scattering amplitude:             fxr_n = (ef.ei)*F0 -i(ef X ei).z_n*F1 + (ef.z_n)(ei.z_n)F2                  Where ei and ef are the initial and final polarisation states, respectively,         and z_n is a unit vector in the direction of the magnetic moment of the nth ion.         The polarisation states are determined to be one of the natural synchrotron          states, where sigma (s) is perpendicular to the scattering plane and pi (p) is         parallel to it.                 ( s-s  s-p )                 ( p-s  p-p )                  From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (15) # Resonant scattering factor # Electric Dipole transition at 3d L edge #F0,F1,F2 = 1,1,1 # Flm (form factor?) # Sigma-Sigma # Sigma-Pi # Pi-Sigma # Pi-Pi )             fxres[ref,:] = F0*f0 -1j*F1*f1 + F2*f2             if disp:                 print('( h, k, l)   TTH  (    mx,    my,    mz)  (    z1,    z2,    z3)')                 fmt = '(%2.0f,%2.0f,%2.0f) %6.2f  (%6.3f,%6.3f,%6.3f)  (%6.3f,%6.3f,%6.3f)  f0=%8.4f  f1=%8.4f  f2=%8.4f fxres= (%8.4f + %8.4fi)'                 for at in range(Nat):                     vals = (HKL[ref,0],HKL[ref,1],HKL[ref,2],tth[ref],                             mxmymz[at,0],mxmymz[at,1],mxmymz[at,2],                             z1[at],z2[at],z3[at],                             f0[at],f1[at],f2[at],                             fxres[ref,at].real,fxres[ref,at].imag)                     print(fmt%vals)         return fxres      def xray_nonresonant_magnetic(self, HKL, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p', disp=False):         """         Calculate the non-resonant magnetic component of the structure factor         for the given HKL, using x-ray rules and form factor           Scattering.xray_magnetic([1,0,0])           Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the real intensity at each reflection.          From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)         Book: "X-ray Scattering and Absorption by Magnetic Materials" by Loevesy and Collins. Ch 2. Eqn.2.21+1         No orbital component assumed         magnetic moments assumed to be in the same reference frame as the polarisation         """          HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])          uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()          kin, kout, ein, eout = self.scatteringvectors(HKL, energy_kev, azim_zero, psi, polarisation)          Qmag = self.xtl.Cell.Qmag(HKL)          # Get Debye-Waller factor         if self._use_isotropic_thermal_factor:             dw = fc.debyewaller(uiso, Qmag)         elif self._use_anisotropic_thermal_factor:             raise Exception('anisotropic thermal factor calculation not implemented yet')         else:             dw = 1          # Get magnetic form factors         if self._use_magnetic_form_factor:             ff = fc.magnetic_form_factor(type, Qmag)         else:             ff = np.ones([len(HKL), len(uvw)])          # Calculate moment         momentmag = fg.mag(mxmymz).reshape([-1, 1])         momentxyz = self.xtl.Cell.calculateR(mxmymz)  # moment direction in cartesian reference frame         moment = momentmag * fg.norm(momentxyz)  # broadcast n*1 x n*3 = n*3         moment[np.isnan(moment)] = 0.          # Magnetic form factor         # f_non-res_mag = i.r0.(hw/mc^2).fD.[.5.L.A + S.B] #equ 2 Hill+McMorrow 1996         # ignore orbital moment L         B = np.zeros([len(HKL), 3])         for n in range(len(HKL)):             #print(n,HKL[n],kin[n],kout[n],ein[n],eout[n])             B[n, :] = np.cross(eout[n], ein[n]) + \                 np.cross(kout[n], eout[n]) * np.dot(kout[n], ein[n]) - \                 np.cross(kin[n], ein[n]) * np.dot(kin[n], eout[n]) - \                 np.cross(np.cross(kout[n], eout[n]), np.cross(kin[n], ein[n]))         fspin = 1j * ff * np.dot(moment, B.T).T          # Calculate dot product         dot_KR = np.dot(HKL, uvw.T)          # Calculate structure factor         SF = np.sum(fspin * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)          if disp:             for n in range(len(HKL)):                 pr_h = '(%2.0f,%2.0f,%2.0f)'%tuple(HKL[n,:])                 pr_b = '(%6.2g,%6.2g,%6.2g)'%tuple(B[n,:])                 print('psi=%3.0f  %3d hkl=%10s B=%22s'%(psi, n, pr_h, pr_b))                 ctot = 0j                 for m in range(len(moment)):                     if np.sum(moment[m,:]**2) < 0.01: continue                     pr_m = '(%6.2g,%6.2g,%6.2g)'%tuple(moment[m,:])                     dot = np.dot(moment[m],B[n,:])                     pdot = '%6.3f'%dot                     phase = np.exp(1j * 2 * np.pi * dot_KR[n, m])                     prph = '%5.2f+i%5.2f'%(np.real(phase),np.imag(phase))                     prod = dot*phase                     pprd = '%5.2f+i%5.2f'%(np.real(prod),np.imag(prod))                     ctot += prod                     ptot = '%5.2f+i%5.2f'%(np.real(ctot),np.imag(ctot))                     print('\t%3d mom=%22s dot(mom,B)=%6s   exp(ik.r)=%12s   sum=%12s   tot=%12s'%(m,pr_m,pdot,prph,pprd,ptot))          SF = SF / self.xtl.scale          if self._return_structure_factor: return SF          # Calculate intensity         I = SF * np.conj(SF)         return np.real(I)      def xray_resonant_magnetic(self, HKL, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p', F0=0, F1=1, F2=0, disp=True):         """         Calculate the non-resonant magnetic component of the structure factor         for the given HKL, using x-ray rules and form factor           Scattering.xray_magnetic([1,0,0])           Scattering.xray_magnetic([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the real intensity at each reflection.          From Hill+McMorrow Acta Cryst. 1996 A52, 236-244 Equ. (2)         Book: "X-ray Scattering and Absorption by Magnetic Materials" by <NAME> Collins. Ch 2. Eqn.2.21+1         No orbital component assumed         magnetic moments assumed to be in the same reference frame as the polarisation         """          HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])          uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()          kin, kout, ein, eout = self.scatteringvectors(HKL, energy_kev, azim_zero, psi, polarisation)          Qmag = self.xtl.Cell.Qmag(HKL)          # Get Debye-Waller factor         if self._use_isotropic_thermal_factor:             dw = fc.debyewaller(uiso, Qmag)         elif self._use_anisotropic_thermal_factor:             raise Exception('anisotropic thermal factor calculation not implemented yet')         else:             dw = 1          # Calculate moment         momentmag = fg.mag(mxmymz).reshape([-1, 1])         momentxyz = self.xtl.Cell.calculateR(mxmymz)  # moment direction in cartesian reference frame         moment = momentmag * fg.norm(momentxyz)  # broadcast n*1 x n*3 = n*3         moment[np.isnan(moment)] = 0.         fe1e1 = np.zeros([len(HKL), len(uvw)], dtype=np.complex)         for ref in range(len(HKL)):             # Magnetic form factor             # f_res_mag = [(e'.e)F0 - i(e'xe).Z*F1 + (e'.Z)*(e.Z)*F2] #equ 7 Hill+McMorrow 1996             f0 = np.dot(eout[ref], ein[ref])             f1 = np.dot(np.cross(eout[ref], ein[ref]), moment.T)             f2 = np.dot(eout[ref], moment.T) * np.dot(ein[ref], moment.T)             fe1e1[ref, :] = f0*F0 - 1j*f1*F1 + f2*F2          # Calculate dot product         dot_KR = np.dot(HKL, uvw.T)          if disp:             for ref in range(len(HKL)):                 pr_h = '(%2.0f,%2.0f,%2.0f)'%tuple(HKL[ref,:])                 print('psi=%3.0f  %3d hkl=%10s'%(psi, ref, pr_h))                 ctot = 0j                 for m in range(len(moment)):                     if np.sum(moment[m,:]**2) < 0.01: continue                     pr_m = '(%6.2g,%6.2g,%6.2g)'%tuple(moment[m,:])                     f0 = np.dot(eout[ref], ein[ref])                     f1 = np.dot(np.cross(eout[ref], ein[ref]), moment[m,:])                     f2 = np.dot(eout[ref], moment[m,:]) * np.dot(ein[ref], moment[m,:])                     dot = f0*F0 - 1j*f1*F1 + f2*F2                     phase = np.exp(1j * 2 * np.pi * dot_KR[ref, m])                     prph = '%5.2f+i%5.2f'%(np.real(phase),np.imag(phase))                     prod = dot*phase                     pprd = '%5.2f+i%5.2f'%(np.real(prod),np.imag(prod))                     ctot += prod                     ptot = '%5.2f+i%5.2f'%(np.real(ctot),np.imag(ctot))                     print('\t%3d mom=%22s f0=%5.2f f1=%5.2f f2=%5.2f   exp(ik.r)=%12s   sum=%12s   tot=%12s'%(m,pr_m,f0,f1,f2,prph,pprd,ptot))          # Calculate structure factor         SF = np.sum(fe1e1 * dw * occ * np.exp(1j * 2 * np.pi * dot_KR), axis=1)          SF = SF / self.xtl.scale          if self._return_structure_factor: return SF          # Calculate intensity         I = SF * np.conj(SF)         return np.real(I)      def scatteringvectors(self, hkl, energy_kev=None, azim_zero=[1, 0, 0], psi=0, polarisation='s-p'):         """         Determine the scattering and polarisation vectors of a reflection based on energy, azimuth and polarisation.         :param xtl: Crystal Class         :param hkl: [n,3] array of reflections         :param energy_kev: x-ray scattering energy in keV         :param azim_zero: [1,3] direction along which the azimuthal zero angle is determind         :param psi: float angle in degrees about the azimuth         :param polarisation: polarisation with respect to the scattering plane, options:                     'ss' : sigma-sigma polarisation                     'sp' : sigma-pi polarisation                     'ps' : pi-sigma polarisation                     'pp' : pi-pi polarisation                 or: polarisation: float polarisation angle of scattered vector in degrees         :return: kin, kout, ein, eout         Returned values are [n,3] arrays             kin : [n,3] array of incident wavevectors             kout: [n,3] array of scattered wavevectors             ein : [n,3] array of incident polarisation             eout: [n,3] array of scattered polarisation          The basis is chosen such that Q defines the scattering plane, sigma and pi directions are normal to this plane.         Q is defined as Q = kout - kin, with kout +ve along the projection of azim_zero         """          hkl = np.asarray(hkl).reshape([-1, 3])          if energy_kev is None:             energy_kev = self._energy_kev          out_kin = np.zeros([len(hkl), 3])         out_kout = np.zeros([len(hkl), 3])         out_ein = np.zeros([len(hkl), 3])         out_eout = np.zeros([len(hkl), 3])         for n in range(len(hkl)):             # Define coordinate system I,J,Q (U1,U2,U3)             Ihat, Jhat, Qhat = self.scatteringbasis(hkl[n, :], azim_zero, psi)             Ihat = Ihat.reshape([-1, 3])             Jhat = Jhat.reshape([-1, 3])             Qhat = Qhat.reshape([-1, 3])              # Determine wavevectors             bragg = self.xtl.Cell.tth(hkl[n, :], energy_kev) / 2.             rb = np.deg2rad(bragg).reshape([-1, 1])             kin = np.cos(rb) * Ihat - np.sin(rb) * Qhat             kout = np.cos(rb) * Ihat + np.sin(rb) * Qhat             esig = Jhat  # sigma polarisation (in or out)             piin = np.cross(kin, esig)  # pi polarisation in             piout = np.cross(kout, esig)  # pi polarisation out              # Polarisation             try:                 # polarisation = 'ss' or 's-s'                 polarisation = polarisation.replace('-', '').replace(' ', '')             except AttributeError:                 # polarisation = angle in deg from sigma' to pi'                 ein = 1.0 * esig                 pol = np.deg2rad(polarisation)                 eout = np.cos(pol)*esig + np.sin(pol)*piout             if polarisation in ['sigmasigma', 'sigsig', 'ss']:                 ein = 1.0 * esig                 eout = 1.0 * esig             elif polarisation in ['sigmapi', 'sigpi', 'sp']:                 ein = 1.0 * esig                 eout = 1.0 * piout             elif polarisation in ['pisigma', 'pisig', 'ps']:                 ein = 1.0 * piin                 eout = 1.0 * esig             elif polarisation in ['pipi', 'pp']:                 ein = 1.0 * piin                 eout = 1.0 * piout             out_kin[n, :] = kin             out_kout[n, :] = kout             out_ein[n, :] = ein             out_eout[n, :] = eout          return out_kin, out_kout, out_ein, out_eout      def scatteringcomponents(self, mxmymz, hkl, azim_zero=[1,0,0], psi=0):         """         Transform magnetic vector into components within the scattering plane             ***warning - may not be correct for non-cubic systems***         """          # Define coordinate system I,J,Q (U1,U2,U3)         U = self.scatteringbasis(hkl, azim_zero, psi)                  # Determine components of the magnetic vector         z1z2z3 = np.dot(mxmymz, U.T) # [mxmymz.I, mxmymz.J, mxmymz.Q]         return fg.norm(z1z2z3)      def scatteringbasis(self, hkl, azim_zero=[1, 0, 0], psi=0):         """         Determine the scattering and polarisation vectors of a reflection based on energy, azimuth and polarisation.         :param hkl: [n,3] array of reflections         :param azim_zero: [1,3] direction along which the azimuthal zero angle is determind         :param psi: float azimuthal angle about U3 in degrees         :return: U1, U2, U3         The basis is chosen such that Q defines the scattering plane, the sigma direction is normal to this plane,         the pi direction is always within this plane.         The azimuthal angle defines a rotation about the Q axis in a clockwise mannor, matching I16.         At an azimuth of 0degrees, U1 is perpendicular to Q, along the direction of azim_zero.         """          # Define coordinate system I,J,Q (U1,U2,U3)         # See FDMNES User's Guide p20 'II-11) Anomalous or resonant diffraction'         # U1 || projection of azim_zero         # U2 _|_ U1,U3         # U3 || Q = kf-ki         azim_zero = fg.norm(self.xtl.Cell.calculateQ(azim_zero)) # put in orthogonal basis         Qhat = fg.norm(self.xtl.Cell.calculateQ(hkl)).reshape([-1,3])  # || Q         AxQ = fg.norm(np.cross(azim_zero, Qhat))         Ihat = fg.norm(np.cross(Qhat, AxQ)).reshape([-1,3])  # || to projection of azim_zero         Jhat = fg.norm(np.cross(Qhat, Ihat)).reshape([-1,3])  # _|_ to I and Q          # Rotate psi about Qhat         rpsi = np.deg2rad(psi)         # -ve sin makes clockwise rotation         # This was checked on 21/1/19 vs CRO paper + sergio's calculations and seems to agree with experiment,         # however we never did an azimuthal scan of the (103) which would have distinguished this completely.         Ihat_psi = fg.norm(np.cos(rpsi) * Ihat - np.sin(rpsi) * Jhat)         Jhat_psi = fg.norm(np.cross(Qhat, Ihat_psi))         return np.vstack([Ihat_psi, Jhat_psi, Qhat])      def print_scattering_coordinates(self,hkl,azim_zero=[1,0,0],psi=0):         """         Transform magnetic vector into components within the scattering plane             ***warning - may not be correct for non-cubic systems***         """                  # Define coordinate system I,J,Q (U1,U2,U3)         Qhat = fg.norm(self.xtl.Cell.calculateQ(hkl)) # || Q         AxQ = fg.norm(np.cross(azim_zero,Qhat))         Ihat = fg.norm(np.cross(Qhat,AxQ)) # || to azim_zero         Jhat = fg.norm(np.cross(Qhat,Ihat)) # -| to I and Q                  # Rotate coordinate system by azimuth         Ihat_psi = fg.norm(np.cos(np.deg2rad(psi))*Ihat + np.sin(np.deg2rad(psi))*Jhat)         Jhat_psi = fg.norm(np.cross(Qhat,Ihat_psi))                  # Determine components of the magnetic vector         U=np.vstack([Ihat_psi,Jhat_psi,Qhat])         print('U1 = (%5.2f,%5.2f,%5.2f)'%(U[0,0],U[0,1],U[0,2]))         print('U2 = (%5.2f,%5.2f,%5.2f)'%(U[1,0],U[1,1],U[1,2]))         print('U3 = (%5.2f,%5.2f,%5.2f)'%(U[2,0],U[2,1],U[2,2]))          def print_intensity(self, HKL):         """         Print intensities calcualted in different ways         """                  HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])         Qmag =  self.xtl.Cell.Qmag(HKL)         srt = np.argsort(Qmag)         HKL = HKL[srt,:]                  IN=self.neutron(HKL)         IX=self.x_ray(HKL)         INM=self.magnetic_neutron(HKL)*1e4         IXM=self.xray_magnetic(HKL)*1e4         IXRss=self.xray_resonant(HKL, None, 'ss')         IXRsp=self.xray_resonant(HKL, None, 'sp')         IXRps=self.xray_resonant(HKL, None, 'ps')         IXRpp=self.xray_resonant(HKL, None, 'pp')                  fmt = '(%2.0f,%2.0f,%2.0f)  %8.1f  %8.1f  %8.2f  %8.2f  ss=%8.2f  sp=%8.2f  ps=%8.2f  pp=%8.2f'         print('( h, k, l)   Neutron      xray   Magn. N  Magn. XR   sig-sig    sig-pi    pi-sig     pi-pi')         for n in range(len(HKL)):             vals=(HKL[n][0],HKL[n][1],HKL[n][2],IN[n],IX[n],INM[n],IXM[n],IXRss[n],IXRsp[n],IXRps[n],IXRpp[n])             print(fmt%vals)          def old_intensity(self, HKL, scattering_type=None):         """         Calculate the squared structure factor for the given HKL           Crystal.intensity([1,0,0])           Crystal.intensity([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the real intensity at each reflection.                  Notes:         - Uses x-ray atomic form factors, calculated from approximated tables in the ITC         - This may be a little slow for large numbers of reflections, as it is not currently          possible to use accelerated calculation methods in Jython.         - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs         - Crystal.scale is used to scale the complex structure factor, so the intensity is            reduced by (Crystal.scale)^2         - Testing against structure factors calculated by Vesta.exe is very close, though there           are some discrepancies, probably due to the method of calculation of the form factor.         """          if scattering_type is None:             scattering_type = self._scattering_type         scattering_type = scattering_type.lower()          # Break up long lists of HKLs         n_arrays = np.ceil(len(HKL)*len(self.xtl.Structure.u)/10000.)         hkl_array = np.array_split(HKL, n_arrays)          intensity = []         for _hkl in hkl_array:             if scattering_type in ['xray','x','x-ray','thomson','charge']:                 intensity += self.x_ray(_hkl).tolist()             elif scattering_type in ['neutron','n','nuclear']:                 intensity += self.neutron(_hkl).tolist()             elif scattering_type in ['xray magnetic','magnetic xray','spin xray','xray spin']:                 intensity += list(self.xray_magnetic(_hkl)*1e4)             elif scattering_type in ['neutron magnetic','magnetic neutron','magnetic']:                 intensity += list(self.magnetic_neutron(_hkl)*1e4)             elif scattering_type in ['xray dispersion']:                 intensity += self.xray_dispersion(_hkl, self._energy_kev).tolist()             elif scattering_type in ['xray resonant','resonant','resonant xray','rxs']:                 intensity += self.xray_resonant(_hkl).tolist()             elif scattering_type in ['xray resonant magnetic', 'xray magnetic resonant',                                      'resonant magnetic', 'magnetic resonant']:                 intensity += self.xray_resonant_magnetic(                     _hkl,                     self._energy_kev,                     self._azimuthal_reference,                     self._azimuthal_angle,                     self._polarisation,                     F0=0, F1=1, F2=0).tolist()             elif scattering_type in ['xray nonresonant magnetic', 'xray magnetic nonresonant',                                      'nonresonant magnetic', 'magnetic nonresonant',                                      'xray non-resonant magnetic', 'xray magnetic non-resonant',                                      'non-resonant magnetic', 'magnetic non-resonant']:                 intensity += self.xray_resonant_magnetic(                     _hkl,                     self._energy_kev,                     self._azimuthal_reference,                     self._azimuthal_angle,                     self._polarisation).tolist()             else:                 print('Scattering type not defined')         return np.array(intensity)      def old_structure_factor(self, HKL, scattering_type=None):         """         Calculate the complex structure factor for the given HKL           Crystal.structure_factor([1,0,0])           Crystal.structure_factor([[1,0,0],[2,0,0],[3,0,0])         Returns an array with the same length as HKL, giving the complex structure factor at each reflection.          Notes:         - Uses x-ray atomic form factors, calculated from approximated tables in the ITC         - This may be a little slow for large numbers of reflections, as it is not currently          possible to use accelerated calculation methods in Jython.         - Debye-Waller factor (atomic displacement) is applied for isotropic ADPs         - Crystal.scale is used to scale the complex structure factor, so the intensity is          reduced by (Crystal.scale)^2         - Testing against structure factors calculated by Vesta.exe is very close, though there           are some discrepancies, probably due to the method of calculation of the form factor.         """         prev_sf_setting = self._return_structure_factor         self._return_structure_factor = True         sf = self.intensity(HKL, scattering_type)         self._return_structure_factor = prev_sf_setting         return sf          def hkl(self, HKL, energy_kev=None):         """ Calculate the two-theta and intensity of the given HKL, display the result"""                  if energy_kev is None:             energy_kev = self._energy_kev                  HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])         tth = self.xtl.Cell.tth(HKL,energy_kev)         inten = self.intensity(HKL)                  print('Energy = %6.3f keV' % energy_kev)         print('( h, k, l) TwoTheta  Intensity')         for n in range(len(tth)):             print('(%2.0f,%2.0f,%2.0f) %8.2f  %9.2f' % (HKL[n,0],HKL[n,1],HKL[n,2],tth[n],inten[n]))          def hkl_reflection(self, HKL, energy_kev=None):         """         Calculate the theta, two-theta and intensity of the given HKL in reflection geometry, display the result         Uses sample orientation set up in setup_scatter         :param HKL: [h,k,l] or list of hkl         :param energy_kev: None or float         :return: str         """                  if energy_kev is None:             energy_kev = self._energy_kev                  HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])         tth = self.xtl.Cell.tth(HKL,energy_kev)         theta = self.xtl.Cell.theta_reflection(HKL, energy_kev, self._scattering_specular_direction, self._scattering_theta_offset)         inten = self.intensity(HKL)                  print('Energy = %6.3f keV' % energy_kev)         print('Specular Direction = (%1.0g,%1.0g,%1.0g)' %               (self._scattering_specular_direction[0],                self._scattering_specular_direction[1],                self._scattering_specular_direction[2]))         print('( h, k, l)    Theta TwoTheta  Intensity')         for n in range(len(tth)):             print('(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f' %                   (HKL[n, 0], HKL[n, 1], HKL[n, 2], theta[n], tth[n], inten[n]))          def hkl_transmission(self,HKL,energy_kev=None):         " Calculate the theta, two-theta and intensity of the given HKL in transmission geometry, display the result"                  if energy_kev is None:             energy_kev = self._energy_kev                  HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])         tth = self.xtl.Cell.tth(HKL,energy_kev)         theta = self.xtl.Cell.theta_transmission(HKL, energy_kev, self._scattering_specular_direction,self._scattering_theta_offset)         inten = self.intensity(HKL)                  print('Energy = %6.3f keV' % energy_kev)         print('Direction parallel to beam  = (%1.0g,%1.0g,%1.0g)' %(self._scattering_parallel_direction[0],self._scattering_parallel_direction[1],self._scattering_parallel_direction[2]))         print('( h, k, l)    Theta TwoTheta  Intensity')         for n in range(len(tth)):             print('(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f' % (HKL[n,0],HKL[n,1],HKL[n,2],theta[n],tth[n],inten[n]))      def generate_powder(self, q_max=8, peak_width=0.01, background=0, powder_average=True):         """         Generates array of intensities along a spaced grid, equivalent to a powder pattern.           Q,I = generate_powder(energy_kev=8.0,peak_width=0.05,background=0)             q_max = maximum Q, in A-1             peak_width = width of convolution, in A-1             background = average of normal background             powder_average = True*/False, apply the powder averaging correction           Returns:             Q = [1000x1] array of wave-vector values             I = [1000x1] array of intensity values          Note: To get two-theta values use:             tth = fc.cal2theta(Q, energy_kev)         Note: To get d-spacing values use:             dspace = fc.caldspace(Q)         """          # Get reflections         hmax, kmax, lmax = fc.maxHKL(q_max, self.xtl.Cell.UVstar())         HKL = fc.genHKL([hmax, -hmax], [kmax, -kmax], [lmax, -lmax])         HKL = self.xtl.Cell.sort_hkl(HKL)  # required for labels         Qmag = self.xtl.Cell.Qmag(HKL)         HKL = HKL[Qmag < q_max, :]         Qmag = self.xtl.Cell.Qmag(HKL)         # Qmag = Qmag[Qmag<q_max]          # Calculate intensities         I = self.intensity(HKL)          if powder_average:             # Apply powder averging correction, I0/|Q|**2             I = I/(Qmag+0.001)**2          # create plotting mesh         pixels = int(2000 * q_max)  # reduce this to make convolution faster         pixel_size = q_max / (1.0 * pixels)         peak_width_pixels = peak_width / (1.0 * pixel_size)         mesh = np.zeros([pixels])         mesh_q = np.linspace(0, q_max, pixels)          # add reflections to background         pixel_coord = Qmag / (1.0 * q_max)         pixel_coord = (pixel_coord * (pixels - 1)).astype(int)          for n in range(1, len(I)):             mesh[pixel_coord[n]] = mesh[pixel_coord[n]] + I[n]          # Convolve with a gaussian (if >0 or not None)         if peak_width:             gauss_x = np.arange(-3*peak_width_pixels, 3*peak_width_pixels + 1)  # gaussian width = 2*FWHM             G = fg.gauss(gauss_x, None, height=1, cen=0, fwhm=peak_width_pixels, bkg=0)             mesh = np.convolve(mesh, G, mode='same')              # Add background (if >0 or not None)         if background:             bkg = np.random.normal(background, np.sqrt(background), [pixels])             mesh = mesh + bkg         return mesh_q, mesh      def powder_correction(self, HKL, intensities, symmetric_multiplyer=True, powder_average=True):         """         Averages symmetric reflections and applies symmetry multipliers and 1/q^2 correction         Ic = I0*C         :param HKL: [nx3] array of [h,k,l] reflections         :param intensities: [nx1] array of reflection intensities         :return: [mx3], [mx1] arrays of averaged, corrected reflections + intensity         """         # Average symmetric reflections         rhkl, rinten = self.xtl.Symmetry.average_symmetric_intensity(HKL, intensities)          if symmetric_multiplyer:             multiplyer = self.xtl.Symmetry.reflection_multiplyer(rhkl)             rinten = rinten*multiplyer          if powder_average:             q = self.xtl.Cell.Qmag(rhkl)             rinten = rinten/(q+0.001)**2          return rhkl, rinten      def print_all_reflections(self, energy_kev=None, print_symmetric=False,                               min_intensity=0.01, max_intensity=None, units=None):         """         Prints a list of all allowed reflections at this energy             energy_kev = energy in keV             print_symmetric = False*/True : omits reflections with the same intensity at the same angle             min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)             max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)             units = None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display         """                  if energy_kev is None:             energy_kev = self._energy_kev                  if min_intensity is None: min_intensity = -1         if max_intensity is None: max_intensity = np.inf                  hkl = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)         if not print_symmetric:             hkl = self.xtl.Symmetry.remove_symmetric_reflections(hkl)         hkl = self.xtl.Cell.sort_hkl(hkl)          tth = self.xtl.Cell.tth(hkl, energy_kev)         inrange = np.all([tth < self._scattering_max_twotheta, tth > self._scattering_min_twotheta], axis=0)         hkl = hkl[inrange, :]         tth = tth[inrange]         #inten = np.sqrt(self.intensity(hkl)) # structure factor         inten = self.intensity(hkl)          if units is None:             units = self._powder_units         units = units.lower()         if units in ['d', 'dspc', 'dspace', 'd space', 'd-space', 'dspacing', 'd spacing', 'd-spacing']:             unit_str = 'd-spacing'             unit = fc.caldspace(tth, energy_kev)         elif units in ['q', 'wavevector']:             unit_str = 'Q'             unit = fc.calqmag(tth, energy_kev)         else:             unit_str = 'TwoTheta'             unit = tth                  fmt = '(%3.0f,%3.0f,%3.0f) %10.2f  %9.2f\n'         outstr = ''                  outstr += 'Energy = %6.3f keV\n' % energy_kev         outstr += 'Radiation: %s\n' % self._scattering_type         outstr += '( h, k, l)    %10s  Intensity\n' % unit_str         #outstr+= fmt % (hkl[0,0],hkl[0,1],hkl[0,2],tth[0],inten[0]) # hkl(0,0,0)         count = 0         for n in range(1, len(tth)):             if inten[n] < min_intensity: continue             if inten[n] > max_intensity: continue             count += 1             outstr += fmt % (hkl[n,0], hkl[n,1], hkl[n,2],unit[n],inten[n])         outstr += 'Reflections: %1.0f\n' % count         return outstr          def print_ref_reflections(self, energy_kev=None, min_intensity=0.01, max_intensity=None):         """         Prints a list of all allowed reflections at this energy in the reflection geometry             energy = energy in keV             min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)             max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)                        |                     // |\                     \\ |                      \\|___/____                       \\   \                        \\__                  Note, to change min/max theta values or specular direciton, change the following attributres of          the crystal object:             self._scattering_max_two_theta   :  maximum detector (two-theta) angle             self._scattering_min_theta       :  minimum sample angle = -opening angle             self._scattering_max_theta       :  maximum sample angle = opening angle             self._scattering_theta_offset    :  sample offset angle             self._scattering_specular_direction : [h,k,l] : reflections normal to sample surface         """                  if energy_kev is None:             energy_kev = self._energy_kev          if min_intensity is None: min_intensity = -1         if max_intensity is None: max_intensity = np.inf          HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)         HKL = self.xtl.Cell.sort_hkl(HKL)         tth = self.xtl.Cell.tth(HKL, energy_kev)         HKL = HKL[tth > self._scattering_min_twotheta, :]         tth = tth[tth > self._scattering_min_twotheta]         theta = self.xtl.Cell.theta_reflection(HKL, energy_kev, self._scattering_specular_direction,                                                self._scattering_theta_offset)         # inten = np.sqrt(self.intensity(HKL)) # structure factor         inten = self.intensity(HKL)          p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)         p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))         pos_theta = p1 * p2                  fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'         outstr = ''          outstr += 'Energy = %6.3f keV\n' % energy_kev         outstr += 'Radiation: %s\n' % self._scattering_type         outstr += 'Specular Direction = (%1.0g,%1.0g,%1.0g)\n' % (             self._scattering_specular_direction[0],             self._scattering_specular_direction[1],             self._scattering_specular_direction[2]         )         outstr += '( h, k, l) TwoTheta    Theta  Intensity\n'         # outstr+= fmt % (HKL[0,0],HKL[0,1],HKL[0,2],tth[0],theta[0],inten[0]) # hkl(0,0,0)         count = 0         for n in range(1, len(tth)):             if inten[n] < min_intensity: continue             if inten[n] > max_intensity: continue             if not pos_theta[n]: continue             # if not print_symmetric and np.abs(tth[n]-tth[n-1]) < 0.01: continue # only works if sorted             count += 1             outstr += fmt % (HKL[n, 0], HKL[n, 1], HKL[n, 2], tth[n], theta[n], inten[n])         outstr += 'Reflections: %1.0f\n' % count         return outstr          def print_tran_reflections(self, energy_kev=None, min_intensity=0.01, max_intensity=None):         """         Prints a list of all allowed reflections at this energy in the transmission geometry             energy = energy in keV             min_intensity = None/ 0.01 : omits reflections less than this (remove extinctions)             max_intensity = None/ 0.01 : omits reflections greater than this (show extinctions only)                    \ /                    --<-- || --<--                    / \                         Note, to change min/max theta values or specular direciton, change the following attributres of          the crystal object:             self._scattering_max_two_theta   :  maximum detector (two-theta) angle             self._scattering_min_theta       :  minimum sample angle = -opening angle             self._scattering_max_theta       :  maximum sample angle = opening angle             self._scattering_theta_offset    :  sample offset angle             self._scattering_parallel_direction : [h,k,l] : reflections normal to sample surface, parallel to beam at theta = 0         """                  if energy_kev is None:             energy_kev = self._energy_kev          if min_intensity is None: min_intensity = -1         if max_intensity is None: max_intensity = np.inf          HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)         HKL = self.xtl.Cell.sort_hkl(HKL)         tth = self.xtl.Cell.tth(HKL, energy_kev)         HKL = HKL[tth > self._scattering_min_twotheta, :]         tth = tth[tth > self._scattering_min_twotheta]         theta = self.xtl.Cell.theta_transmission(HKL, energy_kev, self._scattering_parallel_direction)         # inten = np.sqrt(self.intensity(HKL)) # structure factor         inten = self.intensity(HKL)          p1 = (theta > self._scattering_min_theta) * (theta < self._scattering_max_theta)         p2 = (tth > (theta + self._scattering_min_theta)) * (tth < (theta + self._scattering_max_theta))         pos_theta = p1 * p2          fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'         outstr = ''          outstr += 'Energy = %6.3f keV\n' % energy_kev         outstr += 'Radiation: %s\n' % self._scattering_type         outstr += 'Direction parallel to beam  = (%1.0g,%1.0g,%1.0g)\n' % (             self._scattering_parallel_direction[0],             self._scattering_parallel_direction[1],             self._scattering_parallel_direction[2]         )         outstr += '( h, k, l) TwoTheta    Theta  Intensity\n'         # outstr+= fmt % (HKL[0,0],HKL[0,1],HKL[0,2],tth[0],theta[0],inten[0]) # hkl(0,0,0)         count = 0         for n in range(1, len(tth)):             if inten[n] < min_intensity: continue             if inten[n] > max_intensity: continue             if not pos_theta[n]: continue             # if not print_symmetric and np.abs(tth[n]-tth[n-1]) < 0.01: continue # only works if sorted             count += 1             outstr += fmt % (HKL[n, 0], HKL[n, 1], HKL[n, 2], tth[n], theta[n], inten[n])         outstr += ('Reflections: %1.0f\n' % count)         return outstr          def print_symmetric_reflections(self,HKL):         "Prints equivalent reflections"                  symHKL = self.xtl.Symmetry.symmetric_reflections(HKL)         Nsyms = len(symHKL)         outstr = ''         outstr+= 'Equivalent reflections: %d\n' % Nsyms         for n in range(Nsyms):             outstr+= '(%5.3g,%5.3g,%5.3g)\n' % (symHKL[n,0],symHKL[n,1],symHKL[n,2])         return outstr          def print_atomic_contributions(self,HKL):         """         Prints the atomic contributions to the structure factor         """                  HKL = np.asarray(np.rint(HKL),dtype=np.float).reshape([-1,3])         Nref = len(HKL)                  # Calculate the full intensity         I = self.intensity(HKL)                  # Calculate the structure factors of the symmetric atomic sites         base_label = self.xtl.Atoms.label         uvw,type,label,occ,uiso,mxmymz = self.xtl.Structure.get()                  Qmag = self.xtl.Cell.Qmag(HKL)                  # Get atomic form factors         ff = fc.xray_scattering_factor(type,Qmag)                  # Get Debye-Waller factor         dw = fc.debyewaller(uiso,Qmag)                  # Calculate dot product         dot_KR = np.dot(HKL,uvw.T)                  # Calculate structure factor         SF =  ff*dw*occ*np.exp(1j*2*np.pi*dot_KR)                  # Sum structure factors of each base label in atoms         SFbase = np.zeros([len(HKL),len(base_label)],dtype=np.complex128)         for n in range(len(base_label)):             label_idx = label == base_label[n]             SFbase[:,n] = np.sum(SF[:,label_idx],axis=1)                  # Get the real part of the structure factor         #SFtot = np.sqrt(np.real(SF * np.conj(SF)))         SFrel = np.real(SFbase)         SFimg = np.imag(SFbase)                  # Generate the results         outstr = ''         outstr+= '( h, k, l) Intensity' + ' '.join(['%12s    '%x for x in base_label])+'\n'         for n in range(Nref):             ss = ' '.join(['%6.1f + i%-6.1f' % (x,y) for x,y in zip(SFrel[n],SFimg[n])])             outstr+= '(%2.0f,%2.0f,%2.0f) %9.2f    %s\n' % (HKL[n,0],HKL[n,1],HKL[n,2],I[n],ss)         return outstr      def print_symmetry_contributions(self,HKL):         """         Prints the symmetry contributions to the structure factor for each atomic site         """          HKL = np.asarray(np.rint(HKL), dtype=np.float).reshape([-1, 3])         Nref = len(HKL)          # Calculate the full intensity         I = self.intensity(HKL)          # Calculate the structure factors of the symmetric atomic sites         buvw, btype, base_label, bocc, buiso, bmxmymz = self.xtl.Atoms.get()         operations = np.hstack([self.xtl.Symmetry.symmetric_coordinate_operations(buvw[n]) for n in range(len(buvw))])         rotations = np.hstack(             [self.xtl.Symmetry.symmetric_coordinate_operations(buvw[n], True)[1] for n in range(len(buvw))])          # Calculate the structure factors         uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()         Qmag = self.xtl.Cell.Qmag(HKL)         ff = fc.xray_scattering_factor(type, Qmag)         dw = fc.debyewaller(uiso, Qmag)         dot_KR = np.dot(HKL, uvw.T)         phase = np.exp(1j * 2 * np.pi * dot_KR)         sf = ff * dw * occ * phase          # Generate the results         outstr = ''         for n in range(Nref):             ss = '\n'             all_phase = 0j             all_sf = 0j             for lab in base_label:                 label_idx = np.argwhere(label == lab)                 ss += '  %s\n' % lab                 tot_phase = 0j                 tot_sf = 0j                 for a in label_idx:                     uvwstr = '(%-7.3g,%-7.3g,%-7.3g)' % (uvw[a, 0], uvw[a, 1], uvw[a, 2])                     phstr = fg.complex2str(phase[n, a])                     sfstr = fg.complex2str(sf[n, a])                     val = (a, uvwstr, operations[a[0]], rotations[a[0]], phstr, sfstr)                     ss += '    %3d %s %25s %20s  %s  %s\n' % val                     tot_phase += phase[n, a]                     tot_sf += sf[n, a]                 ss += '%74sTotal:  %s  %s\n' % (' ', fg.complex2str(tot_phase), fg.complex2str(tot_sf))                 all_phase += tot_phase                 all_sf += tot_sf             ss += '%62s Reflection Total:  %s  %s\n' % (' ', fg.complex2str(all_phase), fg.complex2str(all_sf))             outstr += '(%2.0f,%2.0f,%2.0f) I = %9.2f    %s\n' % (HKL[n, 0], HKL[n, 1], HKL[n, 2], I[n], ss)         return outstr          def find_close_reflections(self,HKL,energy_kev=None,max_twotheta=2,max_angle=10):         """         Find and print list of reflections close to the given one         """                  if energy_kev is None:             energy_kev = self._energy_kev                  all_HKL = self.xtl.Cell.all_hkl(energy_kev, self._scattering_max_twotheta)         all_HKL = self.xtl.Cell.sort_hkl(all_HKL)                  all_tth = self.xtl.Cell.tth(all_HKL,energy_kev)         tth1 = self.xtl.Cell.tth(HKL, energy_kev)         tth_dif = np.abs(all_tth-tth1)                  all_Q = self.xtl.Cell.calculateQ(all_HKL)         Q1 = self.xtl.Cell.calculateQ(HKL)         all_angles = np.abs([fg.ang(Q1,Q2,'deg') for Q2 in all_Q])                  selected = (tth_dif < max_twotheta)*(all_angles < max_angle)         sel_HKL = all_HKL[selected,:]         sel_tth = all_tth[selected]         sel_angles = all_angles[selected]         sel_intensity = self.intensity(sel_HKL)                  # Generate Results         fmt = '(%2.0f,%2.0f,%2.0f) %8.2f %8.2f  %9.2f\n'         outstr = ''                   outstr+= 'Energy = %6.3f keV\n' % energy_kev         outstr+= 'Close to Reflection (%2.0f,%2.0f,%2.0f)\n' %(HKL[0],HKL[1],HKL[2])         outstr+= '( h, k, l) TwoTheta    Angle  Intensity\n'         count = 0         for n in range(0,len(sel_HKL)):             count += 1             outstr+= fmt % (sel_HKL[n,0],sel_HKL[n,1],sel_HKL[n,2],sel_tth[n],sel_angles[n],sel_intensity[n])         outstr+= 'Reflections: %1.0f\n'%count         return outstr      def multiple_scattering(self, hkl, azir=[0, 0, 1], pv=[1, 0], energy_range=[7.8, 8.2], numsteps=60,                             full=False, pv1=False, pv2=False, sfonly=True, pv1xsf1=False):         """         Run multiple scattering code, plot result.          mslist = xtl.Scatter.multiple_scattering([h,k,l], energy_range=[7.8, 8.2])          See multiple_scattering.py for more details. Code created by Dr <NAME>, DLS         :param hkl: [h,k,l] principle reflection         :param azir: [h,k,l] reference of azimuthal 0 angle         :param pv: [s,p] polarisation vector         :param energy_range: [min, max] energy range in keV         :param numsteps: int: number of calculation steps from energy min to max         :param full: True/False: calculation type: full         :param pv1: True/False: calculation type: pv1         :param pv2: True/False: calculation type: pv2         :param sfonly: True/False: calculation type: sfonly *default         :param pv1xsf1: True/False: calculation type: pv1xsf1?         :return: array         """         return ms.run_calcms(self.xtl, hkl, azir, pv, energy_range, numsteps,                              full=full, pv1=pv1, pv2=pv2, sfonly=sfonly, pv1xsf1=pv1xsf1)      def ms_azimuth(self, hkl, energy_kev, azir=[0, 0, 1], pv=[1, 0], numsteps=3, peak_width=0.1,                    full=False, pv1=False, pv2=False, sfonly=True, pv1xsf1=False):         """         Returns an azimuthal dependence at a particular energy          :param xtl: Crystal structure from Dans_Diffraction         :param hkl: [h,k,l] principle reflection         :param energy_kev: calculation energy         :param azir: [h,k,l] reference of azimuthal 0 angle         :param pv: [s,p] polarisation vector         :param numsteps: int: number of calculation steps from energy min to max         :param peak_width: convolution width         :param full: True/False: calculation type: full         :param pv1: True/False: calculation type: pv1         :param pv2: True/False: calculation type: pv2         :param sfonly: True/False: calculation type: sfonly *default         :param pv1xsf1: True/False: calculation type: pv1xsf1?         :return: None         """          energy_range = [energy_kev-0.001, energy_kev+0.001]         mslist = self.xtl.Scatter.multiple_scattering(hkl, azir, pv, energy_range, numsteps,                                                       full=full, pv1=pv1, pv2=pv2, sfonly=sfonly, pv1xsf1=pv1xsf1)          if pv1 + pv2 + sfonly + full + pv1xsf1 != 0:             azimuth = np.concatenate(mslist[:, [3, 4]])             intensity = np.concatenate(mslist[:, [-1, -1]])         else:             azimuth = np.concatenate(mslist[:, [3, 4]])             intensity = np.ones(azimuth)          # create plotting mesh         peak_width_pixels = 10         pixels = int(360*(peak_width_pixels/ peak_width))         mesh_azi = np.linspace(-180, 180, pixels)         mesh = np.zeros(mesh_azi.shape)          # add reflections to background         pixel_coord = (azimuth-180) / 360.         pixel_coord = (pixel_coord * pixels).astype(int)          for n in range(len(intensity)):             mesh[pixel_coord[n]] = mesh[pixel_coord[n]] + intensity[n]          # Convolve with a gaussian (if >0 or not None)         if peak_width:             gauss_x = np.arange(-3 * peak_width_pixels, 3 * peak_width_pixels + 1)  # gaussian width = 2*FWHM             G = fg.gauss(gauss_x, None, height=1, cen=0, fwhm=peak_width_pixels, bkg=0)             mesh = np.convolve(mesh, G, mode='same')         return mesh_azi, mesh Return tensor scattering intensities           ss, sp, ps, pp = tensor_scattering('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90)         Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.         :param atom_label: str atom site label, e.g. Ru1         :param hkl: list/array, [h,k,l] reflection to calculate         :param energy_kev: float         :param azir: list/array, [h,k,l] azimuthal reference         :param psideg: float/array, azimuthal angle/ range         :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'         :param rank: int, 1,2,3: tensor rank. Only required         :param time: +/-1 time symmetry         :param parity: +/-1 parity         :param mk:         :param lk:         :param sk:         :return: ss, sp, ps, pp intensity values # Structure data for input # Calculate tensor scattering Return tensor scattering intensities for non-standard polarisation           pol = tensor_scattering_stokes('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90, stokes=45)         Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.         :param atom_label: str atom site label, e.g. Ru1         :param hkl: list/array, [h,k,l] reflection to calculate         :param energy_kev: float         :param azir: list/array, [h,k,l] azimuthal reference         :param psideg: float, azimuthal angle         :param stokes: float/array, rotation of polarisation analyser (0=sigma), degrees         :param pol_theta: float, scattering angle of polarisation analyser, degrees         :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'         :param rank: int, 1,2,3: tensor rank. Only required         :param time: +/-1 time symmetry         :param parity: +/-1 parity         :param mk:         :param lk:         :param sk:         :return: array of intensity values # Structure data for input # Calculate tensor scattering Return tensor scattering intensities           ss, sp, ps, pp = tensor_scattering('Ru1', [0,0,3], 2.838, [0,1,0], psideg=90)         Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.         :param atom_label: str atom site label, e.g. Ru1         :param hkl: list/array, [h,k,l] reflection to calculate         :param energy_kev: float         :param azir: list/array, [h,k,l] azimuthal reference         :param psideg: float, azimuthal angle         :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'         :param rank: int, 1,2,3: tensor rank. Only required         :param time: +/-1 time symmetry         :param parity: +/-1 parity         :param mk:         :param lk:         :param sk:         :return: str # Structure data for input # Calculate tensor scattering Return tensor scattering intensities for all reflections at given azimuth and energy           ss, sp, ps, pp = tensor_scattering('Ru1', 2.838, [0,1,0], psideg=90)         Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.         :param atom_label: str atom site label, e.g. Ru1         :param energy_kev: float         :param azir: list/array, [h,k,l] azimuthal reference         :param psideg: float, azimuthal angle         :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'         :param rank: int, 1,2,3: tensor rank. Only required         :param time: +/-1 time symmetry         :param parity: +/-1 parity         :param mk:         :param lk:         :param sk:         :param print_symmetric: False*/True : omits reflections with the same intensity at the same angle         :param units: None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display         :return: str # Structure data for input Return tensor scattering intensities for all reflections at given energy at maximum intensity psi           ss, sp, ps, pp = tensor_scattering('Ru1', 2.838, [0,1,0], psideg=90)         Uses TensorScatteringClass.py by Prof. <NAME>, Dimaond Light Source Ltd.         :param atom_label: str atom site label, e.g. Ru1         :param energy_kev: float         :param azir: list/array, [h,k,l] azimuthal reference         :param process: str: 'Scalar', 'E1E1', 'E1E2', 'E2E2', 'E1E1mag', 'NonResMag'         :param rank: int, 1,2,3: tensor rank. Only required         :param time: +/-1 time symmetry         :param parity: +/-1 parity         :param mk:         :param lk:         :param sk:         :param print_symmetric: False*/True : omits reflections with the same intensity at the same angle         :param units: None/ 'twotheta'/ 'q'/ 'dspace' : choose scattering angle units to display         :return: str # Caluclate structure factor **2 # Structure data for input Container for available scattering types Calling this container with a scattering type will set the scattering type in the parent class         :param type: str scattering type         :return: None Set max/min angles consistent with reflection geometry on beamline I16         :return: None Set max/min angles consistent with reflection geometry on beamline I16         :return: None Container for scattering type switcher