source: trunk/GSASIIsasd.py @ 1237

Last change on this file since 1237 was 1237, checked in by vondreele, 9 years ago

Allow changes to limits when in SASD Models
"fix" matplotlib warning message when log plot is hit with a mouse button event
There is still an underlying problem.
Mods to SASD Model to have scaled errors on data for fitting purposes - uses wtFactor
Start on size distribution fitting
File size: 10.8 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 & volumes as class definitions
60###############################################################################
61
62class SASDParticles(object):
63    def __init__(self,name=None,parNames=None):
64        self.name = name
65        self.ParmNames = parmNames
66       
67    def FormFactor():
68        FF = 0.
69        return FF
70       
71    def Volume():
72        Vol = 0.
73        return Vol
74       
75class Sphere(SASDParticles):
76    def __init__(self,name=None,parmNames=None):
77        self.Name = name
78        if self.Name == None:
79            self.Name = 'Sphere'
80        self.ParmNames = parmNames
81        if self.ParmNames == None:
82            self.ParmNames = ['Radius',]
83       
84    def FormFactor(self,Q,Parms):
85        ''' Compute hard sphere form factor - can use numpy arrays
86        param float:Q Q value array (usually in A-1)
87        param float:R sphere radius (Usually in A - must match Q-1 units)
88        returns float: form factors as array as needed
89        '''
90        QR = Q*Parms[0]
91        return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
92       
93    def Volume(self,Parms):
94        ''' Compute volume of sphere
95        - numpy array friendly
96        param float:R sphere radius
97        returns float: volume
98        '''
99        return (4./3.)*np.pi*Parms[0]**3
100       
101class Spheroid(SASDParticles):
102    def __init__(self,name=None,parmNames=None):
103        self.Name = name
104        if self.Name == None:
105            self.Name = 'Spheroid'
106        self.ParmNames = parmNames
107        if self.ParmNames == None:
108            self.ParmNames = ['Radius','Aspect ratio']
109   
110    def FormFactor(self,Q,Parms):
111        ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
112        - can use numpy arrays for R & AR; will return corresponding numpy array
113        param float:Q Q value array (usually in A-1)
114        param float R: radius along 2 axes of spheroid
115        param float AR: aspect ratio so 3rd axis = R*AR
116        returns float: form factors as array as needed
117        '''
118        R,AR = Parms
119        NP = 50
120        if 0.99 < AR < 1.01:
121            return SphereFF(Q,R,0)
122        else:
123            cth = np.linspace(0,1.,NP)
124            Rct = R*np.sqrt(1.+(AR**2-1.)*cth**2)
125            return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=1)/NP)
126       
127    def Volume(self,Parms):
128        ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
129        - numpy array friendly
130        param float R: radius along 2 axes of spheroid
131        param float AR: aspect ratio so radius of 3rd axis = R*AR
132        returns float: volume
133        '''
134        R,AR = Parms
135        return AR*(4./3.)*np.pi*R**3
136       
137###############################################################################
138#### Particle form factors
139###############################################################################
140
141def SphereFF(Q,R,dummy=0):
142    ''' Compute hard sphere form factor - can use numpy arrays
143    param float:Q Q value array (usually in A-1)
144    param float:R sphere radius (Usually in A - must match Q-1 units)
145    returns float: form factors as array as needed
146    '''
147    QR = Q*R
148    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
149   
150def SpheroidFF(Q,R,AR):
151    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
152    - can use numpy arrays for R & AR; will return corresponding numpy array
153    param float:Q Q value array (usually in A-1)
154    param float R: radius along 2 axes of spheroid
155    param float AR: aspect ratio so 3rd axis = R*AR
156    returns float: form factors as array as needed
157    '''
158    NP = 50
159    if 0.99 < AR < 1.01:
160        return SphereFF(Q,R,0)
161    else:
162        cth = np.linspace(0,1.,NP)
163        Rct = R*np.sqrt(1.+(AR**2-1.)*cth**2)
164        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=1)/NP)
165           
166def CylinderFF(Q,R,L):
167    ''' Compute form factor for cylinders - can use numpy arrays
168    param float: Q Q value array (A-1)
169    param float: R cylinder radius (A)
170    param float: L cylinder length (A)
171    returns float: form factor
172    '''
173    NP = 200
174    alp = np.linspace(0,np.pi/2.,NP)
175    LBessArg = 0.5*Q[:,np.newaxis]*L*np.cos(alp)
176    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
177    SBessArg = Q[:,np.newaxis]*R*np.sin(alp)
178    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
179    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*(LBess*SBess)**2,axis=1)/NP)
180   
181def CylinderDFF(Q,L,D):
182    ''' Compute form factor for cylinders - can use numpy arrays
183    param float: Q Q value array (A-1)
184    param float: L cylinder length (A)
185    param float: D cylinder diameter (A)
186    returns float: form factor
187    '''
188    return CylinderFF(Q,D/2.,L)   
189   
190def CylinderARFF(Q,R,AR): 
191    ''' Compute form factor for cylinders - can use numpy arrays
192    param float: Q Q value array (A-1)
193    param float: R cylinder radius (A)
194    param float: AR cylinder aspect ratio = L/D = L/2R
195    returns float: form factor
196    '''
197    return CylinderFF(Q,R,2.*R*AR)   
198   
199def UniSphereFF(Q,R,dummy=0):
200    Rg = np.sqrt(3./5.)*R
201    B = 1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
202    QstV = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
203    return np.sqrt(np.exp((-Q**2*Rg**2)/3.)+(B/QstV**4))
204   
205def UniRodFF(Q,R,L):
206    Rg2 = np.sqrt(R**2/2+L**2/12)
207    B2 = np.pi/L
208    Rg1 = np.sqrt(3.)*R/2.
209    G1 = (2./3.)*R/L
210    B1 = 4.*(L+R)/(R**3*L**2)
211    QstV = Q/(scsp.erf(Q*Rg2/np.sqrt(6)))**3
212    FF = np.exp(-Q**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q**2/3.)
213    QstV = Q/(scsp.erf(Q*Rg1/np.sqrt(6)))**3
214    FF += G1*np.exp(-Q**2*Rg1**2/3.)+(B1/QstV**4)
215    return np.sqrt(FF)
216   
217def UniRodARFF(Q,R,AR):
218    return UniRodFF(Q,R,AR*R)
219   
220def UniDiskFF(Q,R,T):
221    Rg2 = np.sqrt(R**2/2.+T**2/12.)
222    B2 = 2./R**2
223    Rg1 = np.sqrt(3.)*T/2.
224    RgC2 = 1.1*Rg1
225    G1 = (2./3.)*(T/R)**2
226    B1 = 4.*(T+R)/(R**3*T**2)
227    QstV = Q/(scsp.erf(Q*Rg2/np.sqrt(6)))**3
228    FF = np.exp(-Q**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q**2/3.)
229    QstV = Q/(scsp.erf(Q*Rg1/np.sqrt(6)))**3
230    FF += G1*np.exp(-Q**2*Rg1**2/3.)+(B1/QstV**4)
231    return np.sqrt(FF)
232   
233def UniTubeFF(Q,R,L,T):
234    Ri = R-T
235    DR2 = R**2-Ri**2
236    Vt = np.pi*DR2*L
237    Rg3 = np.sqrt(DR2/2.+L**2/12.)
238    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
239    B2 = np.pi**2*T/Vt
240    B3 = np.pi/L
241    QstV = Q/(scsp.erf(Q*Rg3/np.sqrt(6)))**3
242    FF = np.exp(-Q**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q**2*R**2/3.)
243    QstV = Q/(scsp.erf(Q*R/np.sqrt(6)))**3
244    FF += (B2/QstV**2)*np.exp(-Q**2*T**2/3.)
245    QstV = Q/(scsp.erf(Q*T/np.sqrt(6)))**3
246    FF += B1/QstV**4
247    return np.sqrt(FF)
248   
249   
250   
251   
252
253###############################################################################
254#### Particle volumes
255###############################################################################
256
257def SphereVol(R):
258    ''' Compute volume of sphere
259    - numpy array friendly
260    param float:R sphere radius
261    returns float: volume
262    '''
263    return (4./3.)*np.pi*R**3
264
265def SpheroidVol(R,AR):
266    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
267    - numpy array friendly
268    param float R: radius along 2 axes of spheroid
269    param float AR: aspect ratio so radius of 3rd axis = R*AR
270    returns float: volume
271    '''
272    return AR*SphereVol(R)
273   
274def CylinderVol(R,L):
275    ''' Compute cylinder volume for radius & length
276    - numpy array friendly
277    param float: R diameter (A)
278    param float: L length (A)
279    returns float:volume (A^3)
280    '''
281    return np.pi*L*R**2
282   
283def CylinderDVol(L,D):
284    ''' Compute cylinder volume for length & diameter
285    - numpy array friendly
286    param float: L length (A)
287    param float: D diameter (A)
288    returns float:volume (A^3)
289    '''
290    return CylinderVol(D/2.,L)
291   
292def CylinderARVol(R,AR):
293    ''' Compute cylinder volume for radius & aspect ratio = L/D
294    - numpy array friendly
295    param float: R radius (A)
296    param float: AR=L/D=L/2R aspect ratio
297    returns float:volume
298    '''
299    return CylinderVol(R,2.*R*AR)
300       
301   
302   
303###############################################################################
304#### Size distribution
305###############################################################################
306
307def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
308    if data['Size']['logBins']:
309        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
310            data['Size']['Nbins']+1,True)
311    else:
312        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
313            data['Size']['Nbins']+1,True)
314    Dbins = np.diff(Bins)
315    BinMag = Dbins*np.ones_like(Dbins)
316    print np.sum(BinMag)
317       
318    print data['Size']
319#    print Limits
320#    print Substances
321#    print Sample
322#    print Profile
323   
Note: See TracBrowser for help on using the repository browser.