Changeset 1234

Timestamp:

Mar 4, 2014 3:24:30 PM (9 years ago)

Author:

vondreele

Message:

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

File:

• trunk/GSASIIsasd.py (modified) (2 diffs)

trunk/GSASIIsasd.py

@@ -21 +21 @@
 import numpy.linalg as nl
 from numpy.fft import ifft, fft, fftshift
+import scipy.special as scsp
 import scipy.interpolate as si
 import scipy.stats as st
@@ -56 +57 @@
     
 ###############################################################################
-#### Particle form factors & volumes
-###############################################################################
-
+#### Particle form factors
+###############################################################################
+
+class SASDParticles(object):
+    def __init__(self,Q,name=None,parms=None,parNames=None):
+        self.Q = Q
+        self.name = name
+        self.Parms = parms
+        self.ParmNames = parmNames
+        
+    def FormFactor():
+        FF = 0.
+        return FF
+        
+    def Volume():
+        Vol = 0.
+        return Vol
+        
+class Sphere(SASDParticles):
+    def __init__(self,Q,name=None,parms=None,parmNames=None):
+        self.Q = Q
+        self.Name = name
+        if self.Name == None:
+            self.Name = 'Sphere'
+        self.Parms = parms
+        if self.Parms == None:
+            self.Parms = [50.0,]
+        self.ParmNames = parmNames
+        if self.ParmNames == None:
+            self.ParmNames = ['Radius',]
+        
+    def FormFactor(self):
+        ''' Compute hard sphere form factor - can use numpy arrays
+        param float:Q Q value array (usually in A-1)
+        param float:R sphere radius (Usually in A - must match Q-1 units)
+        returns float: form factors as array as needed
+        '''
+        QR = self.Q*self.Parms[0]
+        return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
+        
+    def Volume(self):
+        ''' Compute volume of sphere
+        - numpy array friendly
+        param float:R sphere radius
+        returns float: volume
+        '''
+        return (4./3.)*np.pi*self.Parms[0]**3
+        
+        
+    
+
+def SphereFF(Q,R,dummy=0):
+    ''' Compute hard sphere form factor - can use numpy arrays
+    param float:Q Q value array (usually in A-1)
+    param float:R sphere radius (Usually in A - must match Q-1 units)
+    returns float: form factors as array as needed
+    '''
+    QR = Q*R
+    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
+    
 def SpheroidFF(Q,R,AR):
-    '''
-    '''
-    QR = Q*R
-    if 0.98 < AR < 1.02:
-        return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
-        
-    
+    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid) 
+    - can use numpy arrays for R & AR; will return corresponding numpy array
+    param float:Q Q value array (usually in A-1)
+    param float R: radius along 2 axes of spheroid
+    param float AR: aspect ratio so 3rd axis = R*AR
+    returns float: form factors as array as needed
+    '''
+    NP = 50
+    if 0.99 < AR < 1.01:
+        return SphereFF(Q,R,0)
+    else:
+        cth = np.linspace(0,1.,NP)
+        Rct = R*np.sqrt(1.+(AR**2-1.)*cth**2)
+        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=1)/NP)
+            
+def CylinderFF(Q,R,L):
+    ''' Compute form factor for cylinders - can use numpy arrays
+    param float: Q Q value array (A-1)
+    param float: R cylinder radius (A)
+    param float: L cylinder length (A)
+    returns float: form factor
+    '''
+    NP = 200
+    alp = np.linspace(0,np.pi/2.,NP)
+    LBessArg = 0.5*Q[:,np.newaxis]*L*np.cos(alp)
+    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
+    SBessArg = Q[:,np.newaxis]*R*np.sin(alp)
+    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
+    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*(LBess*SBess)**2,axis=1)/NP)
+    
+def CylinderDFF(Q,L,D):
+    ''' Compute form factor for cylinders - can use numpy arrays
+    param float: Q Q value array (A-1)
+    param float: L cylinder length (A)
+    param float: D cylinder diameter (A)
+    returns float: form factor
+    '''
+    return CylinderFF(Q,D/2.,L)    
+    
+def CylinderARFF(Q,R,AR): 
+    ''' Compute form factor for cylinders - can use numpy arrays
+    param float: Q Q value array (A-1)
+    param float: R cylinder radius (A)
+    param float: AR cylinder aspect ratio = L/D = L/2R
+    returns float: form factor
+    '''
+    return CylinderFF(Q,R,2.*R*AR)    
+    
+def UniSphereFF(Q,R,dummy=0):
+    Rg = np.sqrt(3./5.)*R
+    B = 1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
+    QstV = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
+    return np.sqrt(np.exp((-Q**2*Rg**2)/3.)+(B/QstV**4))
+    
+def UniRodFF(Q,R,L):
+    Rg2 = np.sqrt(R**2/2+L**2/12)
+    B2 = np.pi/L
+    Rg1 = np.sqrt(3.)*R/2.
+    G1 = (2./3.)*R/L
+    B1 = 4.*(L+R)/(R**3*L**2)
+    QstV = Q/(scsp.erf(Q*Rg2/np.sqrt(6)))**3
+    FF = np.exp(-Q**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q**2/3.)
+    QstV = Q/(scsp.erf(Q*Rg1/np.sqrt(6)))**3
+    FF += G1*np.exp(-Q**2*Rg1**2/3.)+(B1/QstV**4)
+    return np.sqrt(FF)
+    
+def UniRodARFF(Q,R,AR):
+    return UniRodFF(Q,R,AR*R)
+    
+def UniDiskFF(Q,R,T):
+    Rg2 = np.sqrt(R**2/2.+T**2/12.)
+    B2 = 2./R**2
+    Rg1 = np.sqrt(3.)*T/2.
+    RgC2 = 1.1*Rg1
+    G1 = (2./3.)*(T/R)**2
+    B1 = 4.*(T+R)/(R**3*T**2)
+    QstV = Q/(scsp.erf(Q*Rg2/np.sqrt(6)))**3
+    FF = np.exp(-Q**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q**2/3.)
+    QstV = Q/(scsp.erf(Q*Rg1/np.sqrt(6)))**3
+    FF += G1*np.exp(-Q**2*Rg1**2/3.)+(B1/QstV**4)
+    return np.sqrt(FF)
+    
+def UniTubeFF(Q,R,L,T):
+    Ri = R-T
+    DR2 = R**2-Ri**2
+    Vt = np.pi*DR2*L
+    Rg3 = np.sqrt(DR2/2.+L**2/12.)
+    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
+    B2 = np.pi**2*T/Vt
+    B3 = np.pi/L
+    QstV = Q/(scsp.erf(Q*Rg3/np.sqrt(6)))**3
+    FF = np.exp(-Q**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q**2*R**2/3.)
+    QstV = Q/(scsp.erf(Q*R/np.sqrt(6)))**3
+    FF += (B2/QstV**2)*np.exp(-Q**2*T**2/3.)
+    QstV = Q/(scsp.erf(Q*T/np.sqrt(6)))**3
+    FF += B1/QstV**4
+    return np.sqrt(FF)
+    
+    
+    
+    
+
+###############################################################################
+#### Particle volumes
+###############################################################################
+
+def SphereVol(R):
+    ''' Compute volume of sphere
+    - numpy array friendly
+    param float:R sphere radius
+    returns float: volume
+    '''
+    return (4./3.)*np.pi*R**3
+
 def SpheroidVol(R,AR):
     ''' Compute volume of cylindrically symmetric ellipsoid (spheroid) 
-    - can use numpy arrays for R & AR; will return corresponding numpy array
+    - numpy array friendly
     param float R: radius along 2 axes of spheroid
-    param float AR: aspect ratio so 3rd axis = R*AR
+    param float AR: aspect ratio so radius of 3rd axis = R*AR
     returns float: volume
     '''
-    return AR*(4./3.)*np.pi*R**3
-    
+    return AR*SphereVol(R)
+    
+def CylinderVol(R,L):
+    ''' Compute cylinder volume for radius & length
+    - numpy array friendly
+    param float: R diameter (A)
+    param float: L length (A)
+    returns float:volume (A^3)
+    '''
+    return np.pi*L*R**2
+    
+def CylinderDVol(L,D):
+    ''' Compute cylinder volume for length & diameter
+    - numpy array friendly
+    param float: L length (A)
+    param float: D diameter (A)
+    returns float:volume (A^3)
+    '''
+    return CylinderVol(D/2.,L)
+    
+def CylinderARVol(R,AR):
+    ''' Compute cylinder volume for radius & aspect ratio = L/D
+    - numpy array friendly
+    param float: R radius (A)
+    param float: AR=L/D=L/2R aspect ratio
+    returns float:volume
+    '''
+    return CylinderVol(R,2.*R*AR)
+        
+    
+    
+###############################################################################
+#### Size distribution
+###############################################################################
+
 def SizeDistribution(Profile,Limits,Substances,Sample,data):
     print data['Size']