Ticket #576: hexagonal_cylinder.py

# cylinder model
# Note: model title and parameter table are inserted automatically
r"""

...
"""

import sys
sys.path.append('/usr/local/lib/python2.7/site-packages/')

import numpy as np  # type: ignore
from numpy import pi, inf  # type: ignore
import mpmath
import scipy

name = "hexagonal_cylinder"
title = "Circular cylinder in a hexagonal lattice"
description = """
     
     I(q) = scale*P(q)*S(q) + background
     where S(q) is the structure factor for a hexagonal lattice.
"""
category = "shape-independent"

#             [ "name", "units", default, [lower, upper], "type", "description"],
parameters = [["length", "Ang", 5000, [0, inf], "","Cylinder length"],
              ["radius", "Ang", 150, [0, inf], "","Cylinder radius"],              
              ["a", "Ang", 500, [0, inf], "","Lattice parameter"],
              ["D", "Ang", 5000, [0, inf], "","Size of the ordered domain"],
              ["Rsigma", "", 0.1, [0, inf], "","Polydispersity of radii"],
              ["Asigma", "", 0.002, [0, inf], "","Disorder parameter"],
              ["c", "", 6, [0, inf], "","Correction for Porod invariant"],
             ]

#source = ["hexagonal_cylinder.c"]

def Reflections(hmax):    
    nReflections = (hmax+1)*(hmax+2)/2-1;
    reflections = np.empty([nReflections,4])
    N = -1
    for h in range(1,hmax+1):
        for k in range(h+1):
            N += 1
            reflections[N,0] = h
            reflections[N,1] = k
            reflections[N,3] = 1 # f_hkl
            
            ### m_hkl
            if abs(h)==abs(k):
                reflections[N,2] = 6 
            elif h==0:
                reflections[N,2] = 6
            elif k==0:
                reflections[N,2] = 6
            else:
                reflections[N,2] = 12
    
    return reflections

def LatticeFactor(q,a,D,c):
    d = 2
    n = 1
    delta = 2*np.pi/D
    vd = np.sqrt(3)/2*a**2
    Omega_d = 2*np.pi
    hmax = 7
    
    reflections = Reflections(hmax)
    h_seq = reflections[:,0]
    k_seq = reflections[:,1]
    m_hkl = reflections[:,2]
    f_hkl = reflections[:,3]
    
    Z0 = 0
    for ii in range(h_seq.shape[0]):
        h = h_seq[ii]
        k = k_seq[ii]
        q_hkl = (4*np.pi)/(a*np.sqrt(3))*np.sqrt(h**2 + h*k + k**2)
        x = q - q_hkl
        
        L_hkl = 1/np.pi * (delta/2)/(x**2 + (delta/2)**2) # Lorentzian peak shape
        #L_hkl = 2/(np.pi*delta)*np.exp(-(4*x**2)/(np.pi*delta**2)) # Gaussian peak shape

        dZ =  m_hkl[ii]*(f_hkl[ii]**2)*L_hkl
        Z0 += dZ
    Z0prefactor = (c*(2*np.pi)**(d-1))/(n*vd*Omega_d*q**(d-1))
    Z0 = Z0prefactor*Z0
    
    return Z0

def DebyeWaller(q,a,Asigma):
    abar = a    
    G = np.exp(-(Asigma**2)*(abar**2)*(q**2))
    #G = np.exp(-(q**2))
    
    return G

def betafactor(q,R,Rsigma):
    #beta = np.exp(-(Rsigma**2)*(R**2)*(q**2))
    
    z = 1/(Rsigma**2) - 1
    Rbar = R
    Rz = Rbar/(1+z)
    
    a1 = 3.0/2.0
    a2 = (z+1)/2.0
    a3 = (z+2)/2.0
    b1 = 2.0
    b2 = 3.0
    x1 = -4.0*(q**2)*(Rz**2)
    avgFd2perp = mpmath.hyp3f2(a1,a2,a3,b1,b2,x1)

    x2 = -(q**2)*(Rz**2)
    avgFdperp = mpmath.hyp2f1(a2,a3,b1,x2)
    
    beta = (avgFdperp**2)/avgFd2perp
        
    return beta, avgFd2perp

def Iq(q, length, radius, a, D, Rsigma, Asigma, c):
    """
    :param q:           Input q-value
    :param length:   Intrecept in linear model
    :return:            Calculated Intensity
    """
    
    # Structure factor
    Z0 = LatticeFactor(q,a,D,c)
    G = DebyeWaller(q,a,Asigma)
    beta,avgFd2perp = betafactor(q,radius,Rsigma)

    S = 1 + beta*(Z0-1)*G

    ## Form factor
    si,ci = scipy.special.sici(q*length)
    F2dparallel_or = (2/q*length)*si - (np.sin(q*length/2.0)/(q*length/2))**2
    P = avgFd2perp*F2dparallel_or

    ## Intensity
    I = S*P
    
    return I

def Iqxy(qx, qy, *args):
    """
    :param qx:   Input q_x-value
    :param qy:   Input q_y-value
    :param args: Remaining arguments
    :return:     2D-Intensity
    """
    
    return Iq(qx, *args)*Iq(qy, *args)

    
# ER defaults to 1.0

# VR defaults to 1.0

demo = dict(scale=1, background=0, length=50)