Changeset 2755 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Mar 20, 2017 9:53:00 AM (5 years ago)
Author:
vondreele
Message:

extend Export HKLs for powder data to include most all columns in GSASII.py
narrow HorizontalLine? display
add a ReloadSubstances? option in case user changed wavelength in SASD & REFD substances
Add a default unit scatter to substances - only in new projects
add modeling of REFD patterns for x-rays (not yet tested - no data) & CW neutrons (tested)
modify XScattDen to give clearly fo, f' & f" components
add NCScattDen for CW neutrons - uses wave to generate real/imaginary components
fix printing of scattering density units - now correct symbols
remove unused imports from GSASIIsasd.py

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r2754 r2755  
    6060npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)    #=2pi*d*
    6161ateln2 = 8.0*math.log(2.0)
     62nxs = np.newaxis
    6263
    6364################################################################################
     
    19091910   
    19101911################################################################################
     1912# Reflectometry calculations
     1913################################################################################
     1914
     1915def REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data):
     1916   
     1917    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
     1918    Qmin = Limits[1][0]
     1919    Qmax = Limits[1][1]
     1920    iBeg = np.searchsorted(Q,Qmin)
     1921    iFin = np.searchsorted(Q,Qmax)+1    #include last point
     1922    Ib[:] = data['FltBack'][0]
     1923    Ic[:] = 0
     1924    Scale = data['Scale'][0]
     1925    Nlayers = len(data['Layers'])
     1926    depth = np.zeros(Nlayers)
     1927    rho = np.zeros(Nlayers)
     1928    irho = np.zeros(Nlayers)
     1929    sigma = np.zeros(Nlayers)
     1930    for ilayer,layer in enumerate(data['Layers']):
     1931        name = layer['Name']
     1932        if 'Thick' in layer:    #skips first & last layers
     1933            depth[ilayer] = layer['Thick'][0]
     1934        if 'Rough' in layer:    #skips first layer
     1935            sigma[ilayer] = layer['Rough'][0]
     1936        rho[ilayer] = Substances[name]['Scatt density']*layer['DenMul'][0]
     1937        irho[ilayer] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
     1938        A,B = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
     1939    Ic[iBeg:iFin] = (A**2+B**2)*Scale+Ib[iBeg:iFin]
     1940
     1941def abeles(kz, depth, rho, irho=0, sigma=0):
     1942    """
     1943    Optical matrix form of the reflectivity calculation.
     1944    O.S. Heavens, Optical Properties of Thin Solid Films
     1945   
     1946    Reflectometry as a function of kz for a set of slabs.
     1947
     1948    :Parameters:
     1949
     1950    *kz* : float[n] | |1/Ang|
     1951        Scattering vector $2\pi\sin(\theta)/\lambda$. This is $\tfrac12 Q_z$.
     1952    *depth* :  float[m] | |Ang|
     1953        thickness of each layer.  The thickness of the incident medium
     1954        and substrate are ignored.
     1955    *rho*, *irho* :  float[n,k] | |1e-6/Ang^2|
     1956        real and imaginary scattering length density for each layer for each kz
     1957        Note: absorption cross section mu = 2 irho/lambda for neutrons
     1958    *sigma* : float[m-1] | |Ang|
     1959        interfacial roughness.  This is the roughness between a layer
     1960        and the previous layer. The sigma array should have m-1 entries.
     1961
     1962    Slabs are ordered with the surface SLD at index 0 and substrate at
     1963    index -1, or reversed if kz < 0.
     1964    """
     1965    def calc(kz, depth, rho, irho, sigma):
     1966        if len(kz) == 0: return kz
     1967   
     1968        # Complex index of refraction is relative to the incident medium.
     1969        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
     1970        # in place of kz^2, and ignoring rho_o
     1971        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
     1972        k = kz
     1973   
     1974        # According to Heavens, the initial matrix should be [ 1 F; F 1],
     1975        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
     1976        # multiply versus some coding convenience.
     1977        B11 = 1
     1978        B22 = 1
     1979        B21 = 0
     1980        B12 = 0
     1981        for i in range(0, len(depth)-1):
     1982            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
     1983            F = (k - k_next) / (k + k_next)
     1984            F *= np.exp(-2*k*k_next*sigma[i]**2)
     1985            #print "==== layer",i
     1986            #print "kz:", kz
     1987            #print "k:", k
     1988            #print "k_next:",k_next
     1989            #print "F:",F
     1990            #print "rho:",rho[:,i+1]
     1991            #print "irho:",irho[:,i+1]
     1992            #print "d:",depth[i],"sigma:",sigma[i]
     1993            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
     1994            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
     1995            M21 = F*M11
     1996            M12 = F*M22
     1997            C1 = B11*M11 + B21*M12
     1998            C2 = B11*M21 + B21*M22
     1999            B11 = C1
     2000            B21 = C2
     2001            C1 = B12*M11 + B22*M12
     2002            C2 = B12*M21 + B22*M22
     2003            B12 = C1
     2004            B22 = C2
     2005            k = k_next
     2006   
     2007        r = B12/B11
     2008        return np.real(r),np.imag(r)
     2009
     2010    if np.isscalar(kz): kz = np.asarray([kz], 'd')
     2011
     2012    m = len(depth)
     2013
     2014    # Make everything into arrays
     2015    depth = np.asarray(depth,'d')
     2016    rho = np.asarray(rho,'d')
     2017    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
     2018    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
     2019
     2020    # Repeat rho,irho columns as needed
     2021    if len(rho.shape) == 1:
     2022        rho = rho[None,:]
     2023        irho = irho[None,:]
     2024
     2025    return calc(kz, depth, rho, irho, sigma)
     2026   
     2027################################################################################
    19112028# Stacking fault simulation codes
    19122029################################################################################
Note: See TracChangeset for help on using the changeset viewer.