Changeset 1234 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Mar 4, 2014 3:24:30 PM (8 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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIsasd.py

    r1233 r1234  
    2121import numpy.linalg as nl
    2222from numpy.fft import ifft, fft, fftshift
     23import scipy.special as scsp
    2324import scipy.interpolate as si
    2425import scipy.stats as st
     
    5657   
    5758###############################################################################
    58 #### Particle form factors & volumes
    59 ###############################################################################
    60 
     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   
    61119def SpheroidFF(Q,R,AR):
    62     '''
    63     '''
    64     QR = Q*R
    65     if 0.98 < AR < 1.02:
    66         return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
    67        
    68    
     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
    69234def SpheroidVol(R,AR):
    70235    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
    71     - can use numpy arrays for R & AR; will return corresponding numpy array
     236    - numpy array friendly
    72237    param float R: radius along 2 axes of spheroid
    73     param float AR: aspect ratio so 3rd axis = R*AR
     238    param float AR: aspect ratio so radius of 3rd axis = R*AR
    74239    returns float: volume
    75240    '''
    76     return AR*(4./3.)*np.pi*R**3
    77    
     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
    78276def SizeDistribution(Profile,Limits,Substances,Sample,data):
    79277    print data['Size']
Note: See TracChangeset for help on using the changeset viewer.