source: trunk/GSASIIsasd.py @ 1234

Last change on this file since 1234 was 1234, checked in by vondreele, 8 years ago

create particle formfactors for SASD calculations
tested by plotting them - they match what Irena uses
start alternative using a base class & defined Sphere

File size: 9.0 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII small angle calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2014-01-09 11:09:53 -0600 (Thu, 09 Jan 2014) $
10# $Author: vondreele $
11# $Revision: 1186 $
12# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIIsasd.py $
13# $Id: GSASIIsasd.py 1186 2014-01-09 17:09:53Z vondreele $
14########### SVN repository information ###################
15import sys
16import math
17import time
18
19import numpy as np
20import scipy as sp
21import numpy.linalg as nl
22from numpy.fft import ifft, fft, fftshift
23import scipy.special as scsp
24import scipy.interpolate as si
25import scipy.stats as st
26import scipy.optimize as so
27
28import GSASIIpath
29GSASIIpath.SetVersionNumber("$Revision: 1186 $")
30import GSASIIlattice as G2lat
31import GSASIIspc as G2spc
32import GSASIIElem as G2elem
33import GSASIIgrid as G2gd
34import GSASIIIO as G2IO
35import GSASIImath as G2mth
36import pypowder as pyd
37
38# trig functions in degrees
39sind = lambda x: math.sin(x*math.pi/180.)
40asind = lambda x: 180.*math.asin(x)/math.pi
41tand = lambda x: math.tan(x*math.pi/180.)
42atand = lambda x: 180.*math.atan(x)/math.pi
43atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
44cosd = lambda x: math.cos(x*math.pi/180.)
45acosd = lambda x: 180.*math.acos(x)/math.pi
46rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
47#numpy versions
48npsind = lambda x: np.sin(x*np.pi/180.)
49npasind = lambda x: 180.*np.arcsin(x)/math.pi
50npcosd = lambda x: np.cos(x*math.pi/180.)
51npacosd = lambda x: 180.*np.arccos(x)/math.pi
52nptand = lambda x: np.tan(x*math.pi/180.)
53npatand = lambda x: 180.*np.arctan(x)/np.pi
54npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
55npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
56npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
57   
58###############################################################################
59#### Particle form factors
60###############################################################################
61
62class SASDParticles(object):
63    def __init__(self,Q,name=None,parms=None,parNames=None):
64        self.Q = Q
65        self.name = name
66        self.Parms = parms
67        self.ParmNames = parmNames
68       
69    def FormFactor():
70        FF = 0.
71        return FF
72       
73    def Volume():
74        Vol = 0.
75        return Vol
76       
77class Sphere(SASDParticles):
78    def __init__(self,Q,name=None,parms=None,parmNames=None):
79        self.Q = Q
80        self.Name = name
81        if self.Name == None:
82            self.Name = 'Sphere'
83        self.Parms = parms
84        if self.Parms == None:
85            self.Parms = [50.0,]
86        self.ParmNames = parmNames
87        if self.ParmNames == None:
88            self.ParmNames = ['Radius',]
89       
90    def FormFactor(self):
91        ''' Compute hard sphere form factor - can use numpy arrays
92        param float:Q Q value array (usually in A-1)
93        param float:R sphere radius (Usually in A - must match Q-1 units)
94        returns float: form factors as array as needed
95        '''
96        QR = self.Q*self.Parms[0]
97        return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
98       
99    def Volume(self):
100        ''' Compute volume of sphere
101        - numpy array friendly
102        param float:R sphere radius
103        returns float: volume
104        '''
105        return (4./3.)*np.pi*self.Parms[0]**3
106       
107       
108   
109
110def SphereFF(Q,R,dummy=0):
111    ''' Compute hard sphere form factor - can use numpy arrays
112    param float:Q Q value array (usually in A-1)
113    param float:R sphere radius (Usually in A - must match Q-1 units)
114    returns float: form factors as array as needed
115    '''
116    QR = Q*R
117    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
118   
119def SpheroidFF(Q,R,AR):
120    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
121    - can use numpy arrays for R & AR; will return corresponding numpy array
122    param float:Q Q value array (usually in A-1)
123    param float R: radius along 2 axes of spheroid
124    param float AR: aspect ratio so 3rd axis = R*AR
125    returns float: form factors as array as needed
126    '''
127    NP = 50
128    if 0.99 < AR < 1.01:
129        return SphereFF(Q,R,0)
130    else:
131        cth = np.linspace(0,1.,NP)
132        Rct = R*np.sqrt(1.+(AR**2-1.)*cth**2)
133        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=1)/NP)
134           
135def CylinderFF(Q,R,L):
136    ''' Compute form factor for cylinders - can use numpy arrays
137    param float: Q Q value array (A-1)
138    param float: R cylinder radius (A)
139    param float: L cylinder length (A)
140    returns float: form factor
141    '''
142    NP = 200
143    alp = np.linspace(0,np.pi/2.,NP)
144    LBessArg = 0.5*Q[:,np.newaxis]*L*np.cos(alp)
145    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
146    SBessArg = Q[:,np.newaxis]*R*np.sin(alp)
147    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
148    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*(LBess*SBess)**2,axis=1)/NP)
149   
150def CylinderDFF(Q,L,D):
151    ''' Compute form factor for cylinders - can use numpy arrays
152    param float: Q Q value array (A-1)
153    param float: L cylinder length (A)
154    param float: D cylinder diameter (A)
155    returns float: form factor
156    '''
157    return CylinderFF(Q,D/2.,L)   
158   
159def CylinderARFF(Q,R,AR): 
160    ''' Compute form factor for cylinders - can use numpy arrays
161    param float: Q Q value array (A-1)
162    param float: R cylinder radius (A)
163    param float: AR cylinder aspect ratio = L/D = L/2R
164    returns float: form factor
165    '''
166    return CylinderFF(Q,R,2.*R*AR)   
167   
168def UniSphereFF(Q,R,dummy=0):
169    Rg = np.sqrt(3./5.)*R
170    B = 1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
171    QstV = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
172    return np.sqrt(np.exp((-Q**2*Rg**2)/3.)+(B/QstV**4))
173   
174def UniRodFF(Q,R,L):
175    Rg2 = np.sqrt(R**2/2+L**2/12)
176    B2 = np.pi/L
177    Rg1 = np.sqrt(3.)*R/2.
178    G1 = (2./3.)*R/L
179    B1 = 4.*(L+R)/(R**3*L**2)
180    QstV = Q/(scsp.erf(Q*Rg2/np.sqrt(6)))**3
181    FF = np.exp(-Q**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q**2/3.)
182    QstV = Q/(scsp.erf(Q*Rg1/np.sqrt(6)))**3
183    FF += G1*np.exp(-Q**2*Rg1**2/3.)+(B1/QstV**4)
184    return np.sqrt(FF)
185   
186def UniRodARFF(Q,R,AR):
187    return UniRodFF(Q,R,AR*R)
188   
189def UniDiskFF(Q,R,T):
190    Rg2 = np.sqrt(R**2/2.+T**2/12.)
191    B2 = 2./R**2
192    Rg1 = np.sqrt(3.)*T/2.
193    RgC2 = 1.1*Rg1
194    G1 = (2./3.)*(T/R)**2
195    B1 = 4.*(T+R)/(R**3*T**2)
196    QstV = Q/(scsp.erf(Q*Rg2/np.sqrt(6)))**3
197    FF = np.exp(-Q**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q**2/3.)
198    QstV = Q/(scsp.erf(Q*Rg1/np.sqrt(6)))**3
199    FF += G1*np.exp(-Q**2*Rg1**2/3.)+(B1/QstV**4)
200    return np.sqrt(FF)
201   
202def UniTubeFF(Q,R,L,T):
203    Ri = R-T
204    DR2 = R**2-Ri**2
205    Vt = np.pi*DR2*L
206    Rg3 = np.sqrt(DR2/2.+L**2/12.)
207    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
208    B2 = np.pi**2*T/Vt
209    B3 = np.pi/L
210    QstV = Q/(scsp.erf(Q*Rg3/np.sqrt(6)))**3
211    FF = np.exp(-Q**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q**2*R**2/3.)
212    QstV = Q/(scsp.erf(Q*R/np.sqrt(6)))**3
213    FF += (B2/QstV**2)*np.exp(-Q**2*T**2/3.)
214    QstV = Q/(scsp.erf(Q*T/np.sqrt(6)))**3
215    FF += B1/QstV**4
216    return np.sqrt(FF)
217   
218   
219   
220   
221
222###############################################################################
223#### Particle volumes
224###############################################################################
225
226def SphereVol(R):
227    ''' Compute volume of sphere
228    - numpy array friendly
229    param float:R sphere radius
230    returns float: volume
231    '''
232    return (4./3.)*np.pi*R**3
233
234def SpheroidVol(R,AR):
235    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
236    - numpy array friendly
237    param float R: radius along 2 axes of spheroid
238    param float AR: aspect ratio so radius of 3rd axis = R*AR
239    returns float: volume
240    '''
241    return AR*SphereVol(R)
242   
243def CylinderVol(R,L):
244    ''' Compute cylinder volume for radius & length
245    - numpy array friendly
246    param float: R diameter (A)
247    param float: L length (A)
248    returns float:volume (A^3)
249    '''
250    return np.pi*L*R**2
251   
252def CylinderDVol(L,D):
253    ''' Compute cylinder volume for length & diameter
254    - numpy array friendly
255    param float: L length (A)
256    param float: D diameter (A)
257    returns float:volume (A^3)
258    '''
259    return CylinderVol(D/2.,L)
260   
261def CylinderARVol(R,AR):
262    ''' Compute cylinder volume for radius & aspect ratio = L/D
263    - numpy array friendly
264    param float: R radius (A)
265    param float: AR=L/D=L/2R aspect ratio
266    returns float:volume
267    '''
268    return CylinderVol(R,2.*R*AR)
269       
270   
271   
272###############################################################################
273#### Size distribution
274###############################################################################
275
276def SizeDistribution(Profile,Limits,Substances,Sample,data):
277    print data['Size']
278#    print Limits
279#    print Substances
280#    print Sample
281#    print Profile
282   
Note: See TracBrowser for help on using the repository browser.