Changeset 1234
- Timestamp:
- Mar 4, 2014 3:24:30 PM (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIsasd.py
r1233 r1234 21 21 import numpy.linalg as nl 22 22 from numpy.fft import ifft, fft, fftshift 23 import scipy.special as scsp 23 24 import scipy.interpolate as si 24 25 import scipy.stats as st … … 56 57 57 58 ############################################################################### 58 #### Particle form factors & volumes 59 ############################################################################### 60 59 #### Particle form factors 60 ############################################################################### 61 62 class 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 77 class 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 110 def 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 61 119 def 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 135 def 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 150 def 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 159 def 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 168 def 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 174 def 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 186 def UniRodARFF(Q,R,AR): 187 return UniRodFF(Q,R,AR*R) 188 189 def 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 202 def 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 226 def 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 69 234 def SpheroidVol(R,AR): 70 235 ''' Compute volume of cylindrically symmetric ellipsoid (spheroid) 71 - can use numpy arrays for R & AR; will return corresponding numpy array236 - numpy array friendly 72 237 param float R: radius along 2 axes of spheroid 73 param float AR: aspect ratio so 3rd axis = R*AR238 param float AR: aspect ratio so radius of 3rd axis = R*AR 74 239 returns float: volume 75 240 ''' 76 return AR*(4./3.)*np.pi*R**3 77 241 return AR*SphereVol(R) 242 243 def 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 252 def 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 261 def 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 78 276 def SizeDistribution(Profile,Limits,Substances,Sample,data): 79 277 print data['Size']
Note: See TracChangeset
for help on using the changeset viewer.