def e2pz(w1,w2,th):
	"""
	Calculates the momentum scale and the relativistic Compton cross section 
	correction according to P. Holm, PRA 37, 3706 (1988).
	from KH 29.05.96
	input:
	w1 = incident energy  [keV]
	w2 = scattered energy [keV]
	th = scattering angle [deg]
	returns:
	pz = momentum scale   [a.u.]
	cf = cross section correction factor such that:
	J(pz) = cf * d^2(sigma)/d(w2)*d(Omega) [barn/atom/keV/srad]
	"""
	w1  = np.array(w1)    # make sure arrays are used
	w2  = np.array(w2)           
	m   = constants.value('electron mass energy equivalent in MeV')*1e3 #511.003      # Mass in natural units
	th  = math.radians(th) # th/180.0*np.pi  # Angle to radians
	alp = constants.value('fine-structure constant') #1.0/137.036  # Fine structure constant
	r0  = constants.value('classical electron radius') #2.8179e-15   # Electron radius
	q   = np.sqrt(w1**2 + w2**2-2*w1*w2*np.cos(th))                        # Momentum transfer    
	pz  = q/2.0 - (w1-w2) * np.sqrt(1/4.0 + m**2/(2*w1*w2*(1-np.cos(th)))) # In natural units
	E   = np.sqrt(m**2+pz**2)
	A   = ((w1-w2)*E-w1*w2*(1-np.cos(th)))/q
	D   = (w1-w2*np.cos(th))*A/q
	R   = w1*(E-D)
	R2  = R-w1*w2*(1-np.cos(th))
	chi = R/R2 + R2/R + 2.0*m**2 * (1/R-1/R2) + m**4 * (1/R-1/R2)**2
	cf  = 2.0*w1*q*E/(m**2*r0**2*w2*chi)
	cf  = cf*(1e-28*(m*alp))			                             # Cross section now in barns/atom/keV/srad
	pz  = pz/(m*alp)                                                 # pz to atomic units (a.u.)
	return pz, cf

def electrostatic_units(energy):
    """ coefficient to convert to appropriate electrostatic units """
    area = (1 / ((2 * sqrt(2)) ** 2)) * 2 * sqrt(3)
    factor = (1j * ((2 *
              constants.value("Rydberg constant times hc in eV")) ** 5) *
              1e-5 * 2.08e-15 *
            ((constants.value("lattice parameter of silicon") / 1e-10) ** 3))\
              / (area * ((energy) ** 3))
    return factor

    def loadCIF(self,event):
    
        wildcards = 'XRD cifile (*.cif)|*.cif|All files (*.*)|*.*'
        dlg = wx.FileDialog(self, message='Choose CIF',
                           defaultDir=os.getcwd(),
                           wildcard=wildcards, style=wx.FD_OPEN)

        path, read = None, False
        if dlg.ShowModal() == wx.ID_OK:
            read = True
            path = dlg.GetPath().replace('\\', '/')
        dlg.Destroy()
        
        if read:
            cifile = os.path.split(path)[-1]

            try:
                cif = xu.materials.Crystal.fromCIF(path)
            except:
                print('incorrect file format: %s' % os.path.split(path)[-1])
                return

            ## generate hkl list
            hkllist = []
            maxhkl = 8
            for i in range(-maxhkl,maxhkl+1):
                for j in range(-maxhkl,maxhkl+1):
                    for k in range(-maxhkl,maxhkl+1):
                        if i+j+k > 0: # as long as h,k,l all positive, eliminates 0,0,0
                            hkllist.append([i,j,k])
            
            hc = constants.value(u'Planck constant in eV s') * \
                       constants.value(u'speed of light in vacuum') * 1e-3 ## units: keV-m

            if self.wavelength is not None:
                qlist = cif.Q(hkllist)
                Flist = cif.StructureFactorForQ(qlist,(hc/(self.wavelength*(1e-10))*1e3))
            
                Fall = []
                qall = []
                hklall = []
                for i,hkl in enumerate(hkllist):
                    if np.abs(Flist[i]) > 0.01:
                        Fadd = np.abs(Flist[i])
                        qadd = np.linalg.norm(qlist[i])
                        if qadd not in qall and qadd < 6:
                            Fall.extend((0,Fadd,0))
                            qall.extend((qadd,qadd,qadd))
                if np.shape(Fall)[0] > 0:
                    Fall = np.array(Fall)
                    qall = np.array(qall)
                    self.add1Ddata(qall,Fall,name=os.path.split(path)[-1],cif=True)
                else:
                    print('Could not calculate real structure factors.')
            else:
                print('Wavelength/energy must be specified for structure factor calculations.')

 def __init__(self, frep_MHz = None, n = None):
     if frep_MHz is not None:
         self._frep_MHz = frep_MHz
         if frep_MHz > 1.0e6:
             warnings.warn("frep should be specified in MHz; large value given.")
     if n is not None:
         self.set_NPTS(n)            
     # Constants, moved here so that module runs through Sphynx autodoc when
     # scipy is Mocked out.
     self._c_nmps = constants.value('speed of light in vacuum')*1e9/1e12 # c in nm/ps
     self._c_mks  = constants.value('speed of light in vacuum') # m/s        

    def __init__(self):
    
        ## Constructor
        dialog = wx.Dialog.__init__(self, None, title='Define wavelength/energy')#,size=(500, 440))
        ## remember: size=(width,height)
        
        panel = wx.Panel(self)
        
        main = wx.BoxSizer(wx.VERTICAL)
        hmain1 = wx.BoxSizer(wx.HORIZONTAL)
        
        #####
        ## Energy or Wavelength
        hmain1 = wx.BoxSizer(wx.HORIZONTAL)
        self.ch_EorL = wx.Choice(panel,choices=['Energy (keV)','Wavelength (A)'])
        self.entr_EorL = wx.TextCtrl(panel)#, size=(110, -1))
 
        self.ch_EorL.Bind(wx.EVT_CHOICE, self.onEorLSel)
 
        hmain1.Add(self.ch_EorL,   flag=wx.RIGHT,  border=8)
        hmain1.Add(self.entr_EorL, flag=wx.RIGHT,  border=8)
        
        #####
        ## OKAY!
        hmain2 = wx.BoxSizer(wx.HORIZONTAL)
        #hlpBtn = wx.Button(panel, wx.ID_HELP    )
        okBtn  = wx.Button(panel, wx.ID_OK      )
        canBtn = wx.Button(panel, wx.ID_CANCEL  )

        #hmain2.Add(hlpBtn,flag=wx.RIGHT,  border=8)
        hmain2.Add(canBtn, flag=wx.RIGHT, border=8) 
        hmain2.Add(okBtn,  flag=wx.RIGHT,  border=8)

        main.Add(hmain1, flag=wx.ALL, border=10) 
        main.Add(hmain2, flag=wx.ALL|wx.ALIGN_RIGHT, border=10) 

        panel.SetSizer(main)

        self.Show()
        ix,iy = panel.GetBestSize()
        self.SetSize((ix+20, iy+20))
        
        ## set default
        self.hc = constants.value(u'Planck constant in eV s') * \
                       constants.value(u'speed of light in vacuum') * 1e-3 ## units: keV-m
        self.energy = 19.0        
        self.wavelength = self.hc/(self.energy)*1e10 ## units: A
        self.ch_EorL.SetSelection(0)
        self.entr_EorL.SetValue(str(self.energy))

def convert():
    """ converts xyz file from angstroms to bohrs """
    conversion_factor = constants.angstrom / constants.value("Bohr radius")
    xx_coord, yy_coord, zz_coord = genfromtxt(XYZ_CONVERT, unpack="True")
    xyz = column_stack((xx_coord, yy_coord, zz_coord))
    new_xyz = xyz * 1/conversion_factor
    savetxt(XYZ_CONVERT_OUT, new_xyz, fmt=('%3.15f'), delimiter=' ')

 def __init__(self,Ts,qE0T,phis,lamb,length=22.e-3):
     self.Ts       = Ts
     self.Tsf      = Ts     # final kin. energy
     self.qE0T     = qE0T
     self.phis     = phis   # rad
     self.lamb     = lamb
     self.length   = length
     self.nbslices = 10     # nbof slices
     dl            = length/self.nbslices
     m0c2          = C.value('proton mass energy equivalent in MeV')  # MeV
     g             = 1.+Ts/m0c2   # Lorentz gamma
     bg            = sqrt(g**2-1.)
     beta          = bg/g         # Lorentz beta
     # phases,energies & maps of slices
     phases    = []
     dtks      = []
     self.maps = []
     for i in range(self.nbslices+1):
         phase = self.phis+2*pi/(beta*lamb)*(i*dl-self.length/2.)
         phases.append(phase)
         dtk = qE0T*dl*cos(phase)
         dtks.append(dtk)
     for i in range(self.nbslices):
         dgdp,f,pot,ham = Hamiltonian(self.Tsf,qE0T,phases[i],lamb)
         self.maps.append(Order3Map(+dl,dgdp,f))
         self.Tsf += dtks[i]
     self.Tsf -= dtks[i]     # one too much

def Hamiltonian(Ts,qE0T,phis,lamb):       
    """ 
        kanonische Koordinaten: p (aka w), x (aka phi)
        unabhaegige Variable: t (aka s)
        Hamiltonian H(p,x,t) = g(p) + V(x,t)
        g = A*p**2/2
        V = B*(sin(x)-x*cos(phis))
        f = -dV/dx = -B*(cos(x)-cos(phis))
        dgdp = dH/dp = A*p
    """
    # closure
    m0c2 = C.value('proton mass energy equivalent in MeV')  # MeV
    g    = 1+Ts/m0c2
    bgs  = sqrt(g**2-1.)
    A    = 2.*pi/bgs**3/lamb
    B    = qE0T/m0c2
    # potential
    def V(x,B=B,phis=phis):
        return B*(sin(x)-x*cos(phis))
    # hamiltonian
    def H(p,x,B=B,phis=phis):
        return A*p**2+V(x)
    # dH/dp
    def dgdp(p,A=A):
        return A*p
    # -dV/dx
    def f(x,t,B=B,phis=phis):
        return -B*(cos(x)-cos(phis))

    return dgdp,f,V,H

def freq_units(units):
    """
    Returns conversion factor from THz to the requred units and the label in the form of a namedtuple
    Accepted values: thz, ev, mev, ha, cm-1, cm^-1
    """

    d = {"thz": FreqUnits(1, "THz"),
         "ev": FreqUnits(const.value("hertz-electron volt relationship") * const.tera, "eV"),
         "mev": FreqUnits(const.value("hertz-electron volt relationship") * const.tera / const.milli, "meV"),
         "ha": FreqUnits(const.value("hertz-hartree relationship") * const.tera, "Ha"),
         "cm-1": FreqUnits(const.value("hertz-inverse meter relationship") * const.tera * const.centi, "cm^{-1}"),
         'cm^-1': FreqUnits(const.value("hertz-inverse meter relationship") * const.tera * const.centi, "cm^{-1}")
         }
    try:
        return d[units.lower().strip()]
    except KeyError:
        raise KeyError('Value for units `{}` unknown\nPossible values are:\n {}'.format(units, list(d.keys())))

 def rotational_constants(self, units='ghz'):
     """Compute the rotational constants in 1/cm or GHz."""
     choices = ('invcm', 'ghz')
     units = units.lower()
     if units not in choices:
         raise ValueError("Invalid units, pick one of {}".format(choices))
     principal_moments = self.principal_moments_of_inertia()[0]
     _check_scipy(_found_scipy)
     bohr2ang = spc.value('atomic unit of length') / spc.angstrom
     xfamu = 1 / spc.value('electron mass in u')
     xthz = spc.value('hartree-hertz relationship')
     rotghz = xthz * (bohr2ang ** 2) / (2 * xfamu * spc.giga)
     if units == 'ghz':
         conv = rotghz
     if units == 'invcm':
         ghz2invcm = spc.giga * spc.centi / spc.c
         conv = rotghz * ghz2invcm
     return conv / principal_moments

def water_term_pol(B0, T):
    gamma = constants.value('proton gyromag. ratio')
    hbar = constants.hbar
    k = constants.k

    larmor_w = gamma * B0
    pol = hbar * larmor_w
    pol /= (2 * k * T)
    return pol

    def _read_funcfl_file(self):
        rphi = np.sqrt(27.2*0.529)
        eV2J = spc.value('electron volt')

        print('Reading EAM potential in single-element (funcfl) format')
        f = open(self.filename, 'r')

        line = f.readline()  # comment line
        self.atomic_number = int(f.readline().split()[0])

        line = f.readline().split()
        self.nrho    =   int(line[0]) 
        self.drho    = float(line[1])
        self.nr      =   int(line[2])
        self.dr      = float(line[3])
        self.cutoff  = float(line[4])

        self.fembed  = np.zeros((self.nrho,2))
        self.effchg  = np.zeros((self.nr,2))
        self.fdens   = np.zeros((self.nr,2))

        # read embedding functional
        line    = f.readline().split()
        n_col   = len(line)
        n_lines = int(self.nrho/self.ncol)-1

        i_rows = 0
        while i_rows < self.nrho:
            for val in line:
                self.fembed[i_rows,0] = i_rows*self.drho
                self.fembed[i_rows,1] = float(val)
                i_rows += 1
                if i_rows >= self.nrho: 
                    break
            line=f.readline().split()
        
        # read charge function
        i_rows = 0
        while i_rows < self.nr:
            for val in line:
                self.effchg[i_rows,0]=i_rows*self.dr
                self.effchg[i_rows,1]=float(val)*rphi/(1.e-30+self.effchg[self.i_rows,0])
                i_rows += 1
                if i_rows >= self.nr: 
                    break
            line=f.readline().split()

        # read density function
        i_rows = 0
        while i_rows < self.nr:
            for val in line:
                self.fdens[i_rows,0]=i_rows*self.dr
                self.fdens[i_rows,1]=float(val)
                i_rows += 1
                if i_rows >= self.nr:
                    break
            line=f.readline().split()

    def __init__(self, N=DEFAULT_N):
        """
        Parameters
        ----------
        N : int
            the number of grid points, it should be a power of 2

        Examples
        --------
        >>> from generic_potential import GenericPotential
        >>> device = GenericPotential(1024)
        """

        # default grid size
        self.N = N
        self.hN = int(self.N/2) # half grid
        self.points_before = int(self.N/2)
        self.points_after = int(self.N/2)

        # useful atomic unities
        self.au_l    = cte.value('atomic unit of length')
        self.au_t    = cte.value('atomic unit of time')
        self.au_e    = cte.value('atomic unit of energy')
        self.au_v    = cte.value('atomic unit of electric potential')
        self.hbar_au = 1.0
        self.me_au   = 1.0

        # other useful constants
        self.ev = cte.value('electron volt')
        self.c  = cte.value('speed of light in vacuum') # m/s
        self.me = cte.value('electron mass')
        self.q  = cte.value('elementary charge')

        # relations of interest
        self.au2ev  = self.au_e / self.ev
        self.au2ang = self.au_l / 1e-10

        # specific for default quantum harmonic oscillator
        self.l     = 0.0000081 # m
        self.f     = self.c / self.l # Hz
        self.w     = 2.0 * np.pi * self.f
        self.z_m   = np.linspace(-5e-9,5e-9, self.N)
        self.v_j   = 0.5 * self.me * self.z_m**2 * self.w**2
        self.z_nm  = self.z_m / 1e-9
        self.v_ev  = self.v_j / self.ev
        self.m_eff = np.ones(self.z_nm.size)

        # set default time increment
        self._set_dt()

        # adjust grid
        self.normalize_device()

    def onLoadCIF(self, event=None):

        wildcards = 'CIF file (*.cif)|*.cif|All files (*.*)|*.*'
        dlg = wx.FileDialog(self, message='Choose CIF file',
                           defaultDir=os.getcwd(),
                           wildcard=wildcards, style=wx.FD_OPEN)

        path, read = None, False
        if dlg.ShowModal() == wx.ID_OK:
            read = True
            path = dlg.GetPath().replace('\\', '/')
        dlg.Destroy()

        if read:
            cry_strc = struc_from_cif(path)

            if cry_strc:
                hc = constants.value(u'Planck constant in eV s') * \
                         constants.value(u'speed of light in vacuum') * 1e-3 ## units: keV-m
                energy = hc/(self.xrd.wavelength) ## units: keV
                q,F = calc_all_F(cry_strc,energy,maxhkl=10,qmax=5)

                #self.plot1Dxrd(xrddata, show_xrd2=True)
                print('Values are calculated; plotting not yet implemented.')

 def principal_moments_of_inertia(self, units='amu_bohr_2'):
     """Return the principal moments of inertia in 3 kinds of units:
     1. [amu][bohr]^2
     2. [amu][angstrom]^2
     3. [g][cm]^2
     and the principal axes.
     """
     choices = ('amu_bohr_2', 'amu_angstrom_2', 'g_cm_2')
     units = units.lower()
     if units not in choices:
         raise ValueError("Invalid units, pick one of {}".format(choices))
     moi_tensor = self.moment_of_inertia_tensor()
     principal_moments, principal_axes = np.linalg.eigh(moi_tensor)
     if units == 'amu_bohr_2':
         conv = 1
     if units == 'amu_angstrom_2':
         _check_scipy(_found_scipy)
         bohr2ang = spc.value('atomic unit of length') / spc.angstrom
         conv = bohr2ang ** 2
     if units == 'g_cm_2':
         _check_scipy(_found_scipy)
         amu2g = spc.value('unified atomic mass unit') * spc.kilo
         conv = amu2g * (spc.value('atomic unit of length') * spc.centi) ** 2
     return conv * principal_moments, principal_axes

def test1(Ts,qE0T,phis,lamb):
    print("----------------------Test 1--------------")
    # prep plot
    m0c2   = C.value('proton mass energy equivalent in MeV')  # MeV
    phimin = radians(-60.)
    phimax = -phimin
    # dphi = (phimax-phimin)/20.
    wmin = -0.05/m0c2
    wmax = +0.05/m0c2
    # dw = (wmax-wmin)/20.
    pou    = []
    pin    = []
    nbprtcls  = 5000
    nbgaps    = 40
    # loop particles
    for j in range(nbprtcls):
        # x0 = phi0 = rn.uniform(phimin,phimax)
        # p0 = w0   = rn.uniform(wmin,wmax)
        x0 = phi0 = rn.gauss(phis,2*phis)
        p0 = w0   = rn.gauss(0,2*wmax)
        t0 = 0
        v = (p0,x0,t0)
        pin.append(v)
        gap  = Gap(Ts,qE0T,phis,lamb)
        # loop gaps
        for i in range(nbgaps):
            v = gap.map(v)
            # increase particle impulse for next gap
            gap(gap.Tsf)
        pou.append(v)
    print("Ts {}".format(gap.Tsf))
 
    # plot results
    win   = [x[0]*m0c2     for x in pin]
    phin  = [degrees(x[1]) for x in pin]
    wf    = [x[0]*m0c2     for x in pou]
    phif  = [degrees(x[1]) for x in pou]
 
    ax1   = plt.subplot(121)
    plt.xlim([-100,200])
    plt.ylim([-0.4,1.2])
    ax1.autoscale(enable=False)
    ax1.scatter(phin,win,s=1,c='r')

    ax2   = plt.subplot(122,sharex=ax1,sharey=ax1)
    ax2.autoscale(enable=False)
    ax2.scatter(phif,wf,s=1,c='r')
    plt.show()

    def axis_conversion(self, data, unit, name):
        # Converts the x-axis units to a format Mantid accepts,
        # frequency to wavenumber, time to time of flight

        logger.debug('Axis for conversion: name={0}, unit={1}'.format(name, unit))
        if name == 'frequency' and unit == 'THz':
            logger.information('Axis {0} will be converted to wave number in cm^-1'.format(name))
            unit = 'Energy_inWavenumber'
            data *= sc.tera
            data *= sc.value('Planck constant in eV s')
            data *= 8065.54
        elif name == 'time' and unit == 'ps':
            logger.information('Axis {0} will be converted to time in microsecond'.format(name))
            unit = 'TOF'
            data *= sc.micro
        else:
            unit = 'Empty'
        return (data, unit, name)

def test0(Ts,qE0T,phis,lamb):
    print("----------------------Test 0--------------")
    m0c2 = C.value('proton mass energy equivalent in MeV')  # MeV
    dgdp,dfdx,pot,ham  = Hamiltonian(Ts,qE0T,phis,lamb)
    dt = 1.e-2
    map3p = Order3Map(+dt,dgdp,dfdx)
    map3m = Order3Map(-dt,dgdp,dfdx)

    ax = plt.subplot()
    phimin = 1.5*phis
    phimax = -phis
    dphi = (phimax-phimin)/20.
    phi0 = phimin

    wmin = -0.05/m0c2
    wmax = +0.05/m0c2
    dw = (wmax-wmin)/20.
    w0 = wmin
    # while phi0 <= phimax:
    for i in range(50):
        # phi0 = rn.gauss(degrees(phis),1.*degrees(abs(phis)))
        phi0 = rn.uniform(phimin,phimax)
        w0   = rn.uniform(wmin,wmax)

        t0 = 0.
        p0 = w0
        x0 = phi0
        vp = (p0,x0,t0)
        vm = (p0,x0,t0)
        pou = [vp]
        for step in range(300):
            vp = map3p.map(vp)
            vm = map3m.map(vm)
            pou.append(vp)
            pou.append(vm)
        # phi0 += dphi
        w   = [x[0]*m0c2     for x in pou]
        phi = [degrees(x[1]) for x in pou]
        plt.xlim([-60,60])
        plt.ylim([-0.1,0.1])
        plt.autoscale(enable=False,axis='both')
        ax.scatter(phi,w,s=1,c='r')
    plt.show()

    def _axis_conversion(self, data, unit, name):
        """
        Converts an axis to a Mantid axis type (possibly performing a unit
        conversion).

        @param data The axis data as Numpy array
        @param unit The axis unit as read from the file
        @param name The axis name as read from the file
        @return Tuple containing updated axis details
        """
        logger.debug('Axis for conversion: name={0}, unit={1}'.format(name, unit))

        # Q (nm**-1) to Q (Angstrom**-1)
        if name.lower() == 'q' and unit.lower() == 'inv_nm':
            logger.information('Axis {0} will be converted to Q in Angstrom**-1'.format(name))
            unit = 'MomentumTransfer'
            data /= sc.nano # nm to m
            data *= sc.angstrom # m to Angstrom

        # Frequency (THz) to Energy (meV)
        elif name.lower() == 'frequency' and unit.lower() == 'thz':
            logger.information('Axis {0} will be converted to energy in meV'.format(name))
            unit = 'Energy'
            data *= sc.tera # THz to Hz
            data *= sc.value('Planck constant in eV s') # Hz to eV
            data /= sc.milli # eV to meV

        # Time (ps) to TOF (s)
        elif name.lower() == 'time' and unit.lower() == 'ps':
            logger.information('Axis {0} will be converted to time in microsecond'.format(name))
            unit = 'TOF'
            data *= sc.micro # ps to us

        # No conversion
        else:
            unit = 'Empty'

        return (data, unit, name)

def disorder():
    """ extracts selected rows from input, converts to matrix,
    disorders, and returns matrix which is written to file """
    count = 1
    chosen = [x - 1 for x in ATOMS_DISORDERED]
    blbohr = BOND_LENGTH * constants.angstrom / constants.value("Bohr radius")
    xyz = loadtxt(XYZ_DISORDER)
    data = xyz[chosen]
    new_xyz = xyz.copy()
    final = zip(data, DISORDER_AMOUNT)
    while count <= REPEAT:
        for atom in range(0, len(final)):
            polar = random() * constants.pi
            azimuthal = 2 * polar
            new_xyz[chosen[atom], 0] = final[atom][0][0] + (final[atom][1] *
                    (blbohr / 2) * sin(polar) * cos(azimuthal))
            new_xyz[chosen[atom], 1] = final[atom][0][1] + (final[atom][1] *
                    (blbohr / 2) * sin(polar) * sin(azimuthal))
            new_xyz[chosen[atom], 2] = final[atom][0][2] + (final[atom][1] *
                    (blbohr / 2) * cos(polar))
            out = "data/disorder/si_h_14_" + str(count).zfill(3) + ".xyz"
            savetxt(out, new_xyz, fmt=('%5.14e'), delimiter='\t')
        count += 1