# Python S(q,w) Models

A S(q,w) Python module has to define the two functions "TakinInit" and "TakinSqw".

"TakinInit" is called after one or several parameters have changed (for example after each minimisation step in the convolution fitter). It can be used to check if e.g. pre-calculated tables, variables, etc. need to be recalculated.

The "TakinSqw" function receives four floating-point parameters h,k,l, and E and returns a single floating-point value S, the dynamical structure factor. The function is called for every Monte-Carlo point.

All global variables that are defined in an S(q,w) Python module and that are prefixed with "g_" (for "global") are made available as settable parameters in the convolution dialog and as fit parameters for the convolution fitter.

The minimal interface to Takin is defined as follows: ``

    def TakinInit():
    	# reinitialise variables here
    	pass

    def TakinSqw(h, k, l, E):
    	S = 0.
    	# calculate S here
    	return S

A full example is given below (it can also be found in the subdirectory "examples/sqw_py" of the source code distribution): ``

    #
    # Sample Python S(q,w) module for (anti-)ferromagnetic dispersions
    #
    # @author Tobias Weber
    # @license GPLv2
    # @date jun-2016
    #

    import math as m

    import numpy as np
    import numpy.linalg as la
    from numpy import array	# in global namespace so that Takin can access it

    import scipy as sp
    import scipy.constants as const

    # -----------------------------------------------------------------------------
    # dispersion
    # -----------------------------------------------------------------------------

    # kB in meV/K
    kB = const.k / const.e * 1e3

    # dispersion relations
    def disp_ferro(q, D, offs):
    	return D*q**2. + offs

    def disp_antiferro(q, D, offs):
    	return D*q + offs

    # Gaussian peak
    def gauss(x, x0, sig, amp):
    	norm = (np.sqrt(2.*m.pi) * sig)
    	return amp * np.exp(-0.5*((x-x0)/sig)**2.) / norm

    # Bose factor
    def bose(E, T):
    	n = 1./(m.exp(abs(E)/(kB*T)) - 1.)
    	if E >= 0.:
    		n += 1.
    	return n

    # Bose factor which is cut off below Ecut
    def bose_cutoff(E, T, Ecut=0.02):
    	Ecut = abs(Ecut)

    	if abs(E) < Ecut:
    		b = bose(np.sign(E)*Ecut, T)
    	else:
    		b = bose(E, T)

    	return b

    # -----------------------------------------------------------------------------

    # -----------------------------------------------------------------------------
    # Takin interface
    # -----------------------------------------------------------------------------

    # global variables which can be accessed / changed by Takin
    g_G = np.array([1., 1., 0.])	# Bragg peak

    g_D = 20.		# magnon stiffness
    g_offs = 0.		# energy gap
    g_sig = 0.02		# linewidth
    g_S0 = 1.		# intensity

    g_inc_sig = 0.02	# incoherent width
    g_inc_amp = 1.		# incoherent intensity

    g_T = 300.		# temperature

    g_bose_cut = 0.02	# cutoff energy for Bose factor

    g_disp = 0		# which dispersion (0: ferro./quadr., 1: antiferro./linear)?

    #
    # the init function is called after Takin has changed a global variable (optional)
    #
    def TakinInit():
    	print("Init: G=" + repr(g_G) + ", T=" + repr(g_T))

    #
    # dispersion E(Q) and weight factor (optional)
    #
    def TakinDisp(h, k, l):
    	E_peak = 0.		# energy
    	w_peak = 1.		# weight

    	try:
    		Q = np.array([h,k,l])
    		q = la.norm(Q - g_G)

    		E_peak = 0.
    		if g_disp == 0:
    			E_peak = disp_ferro(q, g_D, g_offs)
    		elif g_disp == 1:
    			E_peak = disp_antiferro(q, g_D, g_offs)
    		else:
    			return [0., 0.]
    	except ZeroDivisionError:
    		return [0., 0.]

    	return [[E_peak, -E_peak], [w_peak, w_peak]]

    #
    # S(Q,E) function, called for every Monte-Carlo point
    #
    def TakinSqw(h, k, l, E):
    	try:
    #		print("h={0}, k={1}, l={2}, E={3}".format(h,k,l,E))
    		[Ep_peak, Em_peak], [wp_peak, wm_peak] = TakinDisp(h,k,l)

    		S_p = gauss(E, Ep_peak, g_sig, g_S0*wp_peak)
    		S_m = gauss(E, Em_peak, g_sig, g_S0*wm_peak)
    		incoh = gauss(E, 0., g_inc_sig, g_inc_amp)

    		S = (S_p + S_m)*bose_cutoff(E, g_T, g_bose_cut) + incoh
    #		print("S={0}".format(S))
    		return S
    	except ZeroDivisionError:
    		return 0.

    # -----------------------------------------------------------------------------

    import os
    print("Script working directory: " + os.getcwd())

    # test
    #print(TakinSqw(1.1, 0.9, 0., 0.4))