source: trunk/GSASIIsasd.py @ 2787

Last change on this file since 2787 was 2787, checked in by vondreele, 5 years ago

recover code removed in error from MakeDiamDist? in GSASIIsasd.py v#2546 Nov 22, 2016!

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 58.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: 2017-04-17 14:29:27 +0000 (Mon, 17 Apr 2017) $
10# $Author: vondreele $
11# $Revision: 2787 $
12# $URL: trunk/GSASIIsasd.py $
13# $Id: GSASIIsasd.py 2787 2017-04-17 14:29:27Z vondreele $
14########### SVN repository information ###################
15import os
16import math
17
18import numpy as np
19import scipy.special as scsp
20import scipy.optimize as so
21#import pdb
22
23import GSASIIpath
24GSASIIpath.SetVersionNumber("$Revision: 2787 $")
25import GSASIIpwd as G2pwd
26
27# trig functions in degrees
28sind = lambda x: math.sin(x*math.pi/180.)
29asind = lambda x: 180.*math.asin(x)/math.pi
30tand = lambda x: math.tan(x*math.pi/180.)
31atand = lambda x: 180.*math.atan(x)/math.pi
32atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
33cosd = lambda x: math.cos(x*math.pi/180.)
34acosd = lambda x: 180.*math.acos(x)/math.pi
35rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
36#numpy versions
37npsind = lambda x: np.sin(x*np.pi/180.)
38npasind = lambda x: 180.*np.arcsin(x)/math.pi
39npcosd = lambda x: np.cos(x*math.pi/180.)
40npacosd = lambda x: 180.*np.arccos(x)/math.pi
41nptand = lambda x: np.tan(x*math.pi/180.)
42npatand = lambda x: 180.*np.arctan(x)/np.pi
43npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
44npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
45npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
46   
47###############################################################################
48#### Particle form factors
49###############################################################################
50
51def SphereFF(Q,R,args=()):
52    ''' Compute hard sphere form factor - can use numpy arrays
53    param float Q: Q value array (usually in A-1)
54    param float R: sphere radius (Usually in A - must match Q-1 units)
55    param array args: ignored
56    returns float: form factors as array as needed
57    '''
58    QR = Q[:,np.newaxis]*R
59    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
60   
61def SphericalShellFF(Q,R,args=()):
62    ''' Compute spherical shell form factor - can use numpy arrays
63    param float Q: Q value array (usually in A-1)
64    param float R: sphere radius (Usually in A - must match Q-1 units)
65    param array args: [float r]: controls the shell thickness: R_inner = min(r*R,R), R_outer = max(r*R,R)
66    returns float: form factors as array as needed
67        Contributed by: L.A. Avakyan, Southern Federal University, Russia
68    '''
69    r = args[0]
70    if r < 0: # truncate to positive value
71        r = 0
72    if r < 1:  # r controls inner sphere radius
73        return SphereFF(Q,R) - SphereFF(Q,R*r)
74    else:      # r controls outer sphere radius
75        return SphereFF(Q,R*r) - SphereFF(Q,R)
76
77def SpheroidFF(Q,R,args):
78    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
79    - can use numpy arrays for R & AR; will return corresponding numpy array
80    param float Q : Q value array (usually in A-1)
81    param float R: radius along 2 axes of spheroid
82    param array args: [float AR]: aspect ratio so 3rd axis = R*AR
83    returns float: form factors as array as needed
84    '''
85    NP = 50
86    AR = args[0]
87    if 0.99 < AR < 1.01:
88        return SphereFF(Q,R,0)
89    else:
90        cth = np.linspace(0,1.,NP)
91        Rct = R[:,np.newaxis]*np.sqrt(1.+(AR**2-1.)*cth**2)
92        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=2)/NP)
93           
94def CylinderFF(Q,R,args):
95    ''' Compute form factor for cylinders - can use numpy arrays
96    param float Q: Q value array (A-1)
97    param float R: cylinder radius (A)
98    param array args: [float L]: cylinder length (A)
99    returns float: form factor
100    '''
101    L = args[0]
102    NP = 200
103    alp = np.linspace(0,np.pi/2.,NP)
104    QL = Q[:,np.newaxis]*L
105    LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)
106    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
107    QR = Q[:,np.newaxis]*R
108    SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
109    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
110    LSBess = LBess*SBess
111    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
112   
113def CylinderDFF(Q,L,args):
114    ''' Compute form factor for cylinders - can use numpy arrays
115    param float Q: Q value array (A-1)
116    param float L: cylinder half length (A)
117    param array args: [float R]: cylinder radius (A)
118    returns float: form factor
119    '''
120    R = args[0]
121    return CylinderFF(Q,R,args=[2.*L])   
122   
123def CylinderARFF(Q,R,args): 
124    ''' Compute form factor for cylinders - can use numpy arrays
125    param float Q: Q value array (A-1)
126    param float R: cylinder radius (A)
127    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
128    returns float: form factor
129    '''
130    AR = args[0]
131    return CylinderFF(Q,R,args=[2.*R*AR])   
132   
133def UniSphereFF(Q,R,args=0):
134    ''' Compute form factor for unified sphere - can use numpy arrays
135    param float Q: Q value array (A-1)
136    param float R: cylinder radius (A)
137    param array args: ignored
138    returns float: form factor
139    '''
140    Rg = np.sqrt(3./5.)*R
141    B = np.pi*1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
142    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
143    return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
144   
145def UniRodFF(Q,R,args):
146    ''' Compute form factor for unified rod - can use numpy arrays
147    param float Q: Q value array (A-1)
148    param float R: cylinder radius (A)
149    param array args: [float R]: cylinder radius (A)
150    returns float: form factor
151    '''
152    L = args[0]
153    Rg2 = np.sqrt(R**2/2+L**2/12)
154    B2 = np.pi/L
155    Rg1 = np.sqrt(3.)*R/2.
156    G1 = (2./3.)*R/L
157    B1 = 4.*(L+R)/(R**3*L**2)
158    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
159    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q[:,np.newaxis]**2/3.)
160    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
161    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
162    return np.sqrt(FF)
163   
164def UniRodARFF(Q,R,args):
165    ''' Compute form factor for unified rod of fixed aspect ratio - can use numpy arrays
166    param float Q: Q value array (A-1)
167    param float R: cylinder radius (A)
168    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
169    returns float: form factor
170    '''
171    AR = args[0]
172    return UniRodFF(Q,R,args=[2.*AR*R,])
173   
174def UniDiskFF(Q,R,args):
175    ''' Compute form factor for unified disk - can use numpy arrays
176    param float Q: Q value array (A-1)
177    param float R: cylinder radius (A)
178    param array args: [float T]: disk thickness (A)
179    returns float: form factor
180    '''
181    T = args[0]
182    Rg2 = np.sqrt(R**2/2.+T**2/12.)
183    B2 = 2./R**2
184    Rg1 = np.sqrt(3.)*T/2.
185    RgC2 = 1.1*Rg1
186    G1 = (2./3.)*(T/R)**2
187    B1 = 4.*(T+R)/(R**3*T**2)
188    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
189    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q[:,np.newaxis]**2/3.)
190    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
191    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
192    return np.sqrt(FF)
193   
194def UniTubeFF(Q,R,args):
195    ''' Compute form factor for unified tube - can use numpy arrays
196    assumes that core of tube is same as the matrix/solvent so contrast
197    is from tube wall vs matrix
198    param float Q: Q value array (A-1)
199    param float R: cylinder radius (A)
200    param array args: [float L,T]: tube length & wall thickness(A)
201    returns float: form factor
202    '''
203    L,T = args[:2]
204    Ri = R-T
205    DR2 = R**2-Ri**2
206    Vt = np.pi*DR2*L
207    Rg3 = np.sqrt(DR2/2.+L**2/12.)
208    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
209    B2 = np.pi**2*T/Vt
210    B3 = np.pi/L
211    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
212    FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
213    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
214    FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
215    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
216    FF += B1/QstV**4
217    return np.sqrt(FF)
218
219###############################################################################
220#### Particle volumes
221###############################################################################
222
223def SphereVol(R,args=()):
224    ''' Compute volume of sphere
225    - numpy array friendly
226    param float R: sphere radius
227    param array args: ignored
228    returns float: volume
229    '''
230    return (4./3.)*np.pi*R**3
231
232def SphericalShellVol(R,args):
233    ''' Compute volume of spherical shell
234    - numpy array friendly
235    param float R: sphere radius
236    param array args: [float r]: controls shell thickness, see SphericalShellFF description
237    returns float: volume
238    '''
239    r = args[0]
240    if r < 0:
241        r = 0
242    if r < 1:
243        return SphereVol(R) - SphereVol(R*r)
244    else:
245        return SphereVol(R*r) - SphereVol(R)
246
247def SpheroidVol(R,args):
248    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
249    - numpy array friendly
250    param float R: radius along 2 axes of spheroid
251    param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR
252    returns float: volume
253    '''
254    AR = args[0]
255    return AR*SphereVol(R)
256   
257def CylinderVol(R,args):
258    ''' Compute cylinder volume for radius & length
259    - numpy array friendly
260    param float R: diameter (A)
261    param array args: [float L]: length (A)
262    returns float:volume (A^3)
263    '''
264    L = args[0]
265    return np.pi*L*R**2
266   
267def CylinderDVol(L,args):
268    ''' Compute cylinder volume for length & diameter
269    - numpy array friendly
270    param float: L half length (A)
271    param array args: [float D]: diameter (A)
272    returns float:volume (A^3)
273    '''
274    D = args[0]
275    return CylinderVol(D/2.,[2.*L,])
276   
277def CylinderARVol(R,args):
278    ''' Compute cylinder volume for radius & aspect ratio = L/D
279    - numpy array friendly
280    param float: R radius (A)
281    param array args: [float AR]: =L/D=L/2R aspect ratio
282    returns float:volume
283    '''
284    AR = args[0]
285    return CylinderVol(R,[2.*R*AR,])
286   
287def UniSphereVol(R,args=()):
288    ''' Compute volume of sphere
289    - numpy array friendly
290    param float R: sphere radius
291    param array args: ignored
292    returns float: volume
293    '''
294    return SphereVol(R)
295   
296def UniRodVol(R,args):
297    ''' Compute cylinder volume for radius & length
298    - numpy array friendly
299    param float R: diameter (A)
300    param array args: [float L]: length (A)
301    returns float:volume (A^3)
302    '''
303    L = args[0]
304    return CylinderVol(R,[L,])
305   
306def UniRodARVol(R,args):
307    ''' Compute rod volume for radius & aspect ratio
308    - numpy array friendly
309    param float R: diameter (A)
310    param array args: [float AR]: =L/D=L/2R aspect ratio
311    returns float:volume (A^3)
312    '''
313    AR = args[0]
314    return CylinderARVol(R,[AR,])
315   
316def UniDiskVol(R,args):
317    ''' Compute disk volume for radius & thickness
318    - numpy array friendly
319    param float R: diameter (A)
320    param array args: [float T]: thickness
321    returns float:volume (A^3)
322    '''
323    T = args[0]
324    return CylinderVol(R,[T,])
325   
326def UniTubeVol(R,args):
327    ''' Compute tube volume for radius, length & wall thickness
328    - numpy array friendly
329    param float R: diameter (A)
330    param array args: [float L,T]: tube length & wall thickness(A)
331    returns float: volume (A^3) of tube wall
332    '''
333    L,T = args[:2]
334    return CylinderVol(R,[L,])-CylinderVol(R-T,[L,])
335   
336################################################################################
337#### Distribution functions & their cumulative fxns
338################################################################################
339
340def LogNormalDist(x,pos,args):
341    ''' Standard LogNormal distribution - numpy friendly on x axis
342    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
343    param float x: independent axis (can be numpy array)
344    param float pos: location of distribution
345    param float scale: width of distribution (m)
346    param float shape: shape - (sigma of log(LogNormal))
347    returns float: LogNormal distribution
348    '''
349    scale,shape = args
350    return np.exp(-np.log((x-pos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(x-pos)*shape)
351   
352def GaussDist(x,pos,args):
353    ''' Standard Normal distribution - numpy friendly on x axis
354    param float x: independent axis (can be numpy array)
355    param float pos: location of distribution
356    param float scale: width of distribution (sigma)
357    param float shape: not used
358    returns float: Normal distribution
359    '''
360    scale = args[0]
361    return (1./(scale*np.sqrt(2.*np.pi)))*np.exp(-(x-pos)**2/(2.*scale**2))
362   
363def LSWDist(x,pos,args=[]):
364    ''' Lifshitz-Slyozov-Wagner Ostwald ripening distribution - numpy friendly on x axis
365    ref:
366    param float x: independent axis (can be numpy array)
367    param float pos: location of distribution
368    param float scale: not used
369    param float shape: not used
370    returns float: LSW distribution   
371    '''
372    redX = x/pos
373    result = (81.*2**(-5/3.))*(redX**2*np.exp(-redX/(1.5-redX)))/((1.5-redX)**(11/3.)*(3.-redX)**(7/3.))
374   
375    return np.nan_to_num(result/pos)
376   
377def SchulzZimmDist(x,pos,args):
378    ''' Schulz-Zimm macromolecule distribution - numpy friendly on x axis
379    ref: http://goldbook.iupac.org/S05502.html
380    param float x: independent axis (can be numpy array)
381    param float pos: location of distribution
382    param float scale: width of distribution (sigma)
383    param float shape: not used
384    returns float: Schulz-Zimm distribution
385    '''
386    scale = args[0]
387    b = (2.*pos/scale)**2
388    a = b/pos
389    if b<70:    #why bother?
390        return (a**(b+1.))*x**b*np.exp(-a*x)/scsp.gamma(b+1.)
391    else:
392        return np.exp((b+1.)*np.log(a)-scsp.gammaln(b+1.)+b*np.log(x)-(a*x))
393           
394def LogNormalCume(x,pos,args):
395    ''' Standard LogNormal cumulative distribution - numpy friendly on x axis
396    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
397    param float x: independent axis (can be numpy array)
398    param float pos: location of distribution
399    param float scale: width of distribution (sigma)
400    param float shape: shape parameter
401    returns float: LogNormal cumulative distribution
402    '''
403    scale,shape = args
404    lredX = np.log((x-pos)/scale)
405    return (scsp.erf((lredX/shape)/np.sqrt(2.))+1.)/2.
406   
407def GaussCume(x,pos,args):
408    ''' Standard Normal cumulative distribution - numpy friendly on x axis
409    param float x: independent axis (can be numpy array)
410    param float pos: location of distribution
411    param float scale: width of distribution (sigma)
412    param float shape: not used
413    returns float: Normal cumulative distribution
414    '''
415    scale = args[0]
416    redX = (x-pos)/scale
417    return (scsp.erf(redX/np.sqrt(2.))+1.)/2.
418   
419def LSWCume(x,pos,args=[]):
420    ''' Lifshitz-Slyozov-Wagner Ostwald ripening cumulative distribution - numpy friendly on x axis
421    param float x: independent axis (can be numpy array)
422    param float pos: location of distribution
423    param float scale: not used
424    param float shape: not used
425    returns float: LSW cumulative distribution
426    '''
427    nP = 500
428    xMin,xMax = [0.,2.*pos]
429    X = np.linspace(xMin,xMax,nP,True)
430    fxn = LSWDist(X,pos)
431    mat = np.outer(np.ones(nP),fxn)
432    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
433    return np.interp(x,X,cume,0,1)
434   
435def SchulzZimmCume(x,pos,args):
436    ''' Schulz-Zimm cumulative distribution - numpy friendly on x axis
437    param float x: independent axis (can be numpy array)
438    param float pos: location of distribution
439    param float scale: width of distribution (sigma)
440    param float shape: not used
441    returns float: Normal distribution
442    '''
443    scale = args[0]
444    nP = 500
445    xMin = np.max([0.,pos-4.*scale])
446    xMax = np.min([pos+4.*scale,1.e5])
447    X = np.linspace(xMin,xMax,nP,True)
448    fxn = LSWDist(X,pos)
449    mat = np.outer(np.ones(nP),fxn)
450    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
451    return np.interp(x,X,cume,0,1)
452   
453    return []
454   
455################################################################################
456#### Structure factors for condensed systems
457################################################################################
458
459def DiluteSF(Q,args=[]):
460    ''' Default: no structure factor correction for dilute system
461    '''
462    return np.ones_like(Q)  #or return 1.0
463
464def HardSpheresSF(Q,args):
465    ''' Computes structure factor for not dilute monodisperse hard spheres
466    Refs.: PERCUS,YEVICK PHYS. REV. 110 1 (1958),THIELE J. CHEM PHYS. 39 474 (1968),
467    WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
468   
469    param float Q: Q value array (A-1)
470    param array args: [float R, float VolFrac]: interparticle distance & volume fraction
471    returns numpy array S(Q)
472    '''
473   
474    R,VolFr = args
475    denom = (1.-VolFr)**4
476    num = (1.+2.*VolFr)**2
477    alp = num/denom
478    bet = -6.*VolFr*(1.+VolFr/2.)**2/denom
479    gamm = 0.5*VolFr*num/denom       
480    A = 2.*Q*R
481    A2 = A**2
482    A3 = A**3
483    A4 = A**4
484    Rca = np.cos(A)
485    Rsa = np.sin(A)
486    calp = alp*(Rsa/A2-Rca/A)
487    cbet = bet*(2.*Rsa/A2-(A2-2.)*Rca/A3-2./A3)
488    cgam = gamm*(-Rca/A+(4./A)*((3.*A2-6.)*Rca/A4+(A2-6.)*Rsa/A3+6./A4))
489    pfac = -24.*VolFr/A
490    C = pfac*(calp+cbet+cgam)
491    return 1./(1.-C)
492       
493def SquareWellSF(Q,args):
494    '''Computes structure factor for not dilute monodisperse hard sphere with a
495    square well potential interaction.
496    Refs.: SHARMA,SHARMA, PHYSICA 89A,(1977),213-
497   
498    :param float Q: Q value array (A-1)
499    :param array args: [float R, float VolFrac, float depth, float width]:
500        interparticle distance, volume fraction (<0.08), well depth (e/kT<1.5kT),
501        well width
502    :returns: numpy array S(Q)
503      well depth > 0 attractive & values outside above limits nonphysical cf.
504      Monte Carlo simulations
505    '''
506    R,VolFr,Depth,Width = args
507    eta = VolFr
508    eta2 = eta*eta
509    eta3 = eta*eta2
510    eta4 = eta*eta3       
511    etam1 = 1. - eta
512    etam14 = etam1**4
513    alp = (  ( (1. + 2.*eta)**2 ) + eta3*( eta-4.0 )  )/etam14
514    bet = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14
515    gam = 0.5*eta*( (1. + 2.*eta)**2 + eta3*(eta-4.) )/etam14
516
517    SK = 2.*Q*R
518    SK2 = SK*SK
519    SK3 = SK*SK2
520    SK4 = SK3*SK
521    T1 = alp * SK3 * ( np.sin(SK) - SK * np.cos(SK) )
522    T2 = bet * SK2 * ( 2.*SK*np.sin(SK) - (SK2-2.)*np.cos(SK) - 2.0 )
523    T3 =   ( 4.0*SK3 - 24.*SK ) * np.sin(SK) 
524    T3 = T3 - ( SK4 - 12.0*SK2 + 24.0 )*np.cos(SK) + 24.0   
525    T3 = gam*T3
526    T4 = -Depth*SK3*(np.sin(Width*SK) - Width*SK*np.cos(Width*SK)+ SK*np.cos(SK) - np.sin(SK) )
527    CK =  -24.0*eta*( T1 + T2 + T3 + T4 )/SK3/SK3
528    return 1./(1.-CK)
529   
530def StickyHardSpheresSF(Q,args):
531    ''' Computes structure factor for not dilute monodisperse hard spheres
532    Refs.: PERCUS,YEVICK PHYS. REV. 110 1 (1958),THIELE J. CHEM PHYS. 39 474 (1968),
533    WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
534   
535    param float Q: Q value array (A-1)
536    param array args: [float R, float VolFrac]: sphere radius & volume fraction
537    returns numpy array S(Q)
538    '''
539    R,VolFr,epis,sticky = args
540    eta = VolFr/(1.0-epis)/(1.0-epis)/(1.0-epis)       
541    sig = 2.0 * R
542    aa = sig/(1.0 - epis)
543    etam1 = 1.0 - eta
544#  SOLVE QUADRATIC FOR LAMBDA
545    qa = eta/12.0
546    qb = -1.0*(sticky + eta/(etam1))
547    qc = (1.0 + eta/2.0)/(etam1*etam1)
548    radic = qb*qb - 4.0*qa*qc
549#   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
550    lam1 = (-1.0*qb-np.sqrt(radic))/(2.0*qa)
551    lam2 = (-1.0*qb+np.sqrt(radic))/(2.0*qa)
552    lam = min(lam1,lam2)
553    mu = lam*eta*etam1
554    alp = (1.0 + 2.0*eta - mu)/(etam1*etam1)
555    bet = (mu - 3.0*eta)/(2.0*etam1*etam1)
556#   CALCULATE THE STRUCTURE FACTOR<P></P>
557    kk = Q*aa
558    k2 = kk*kk
559    k3 = kk*k2
560    ds = np.sin(kk)
561    dc = np.cos(kk)
562    aq1 = ((ds - kk*dc)*alp)/k3
563    aq2 = (bet*(1.0-dc))/k2
564    aq3 = (lam*ds)/(12.0*kk)
565    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
566
567    bq1 = alp*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
568    bq2 = bet*(1.0/kk - ds/k2)
569    bq3 = (lam/12.0)*((1.0 - dc)/kk)
570    bq = 12.0*eta*(bq1+bq2-bq3)
571    sq = 1.0/(aq*aq +bq*bq)
572
573    return sq
574
575#def HayterPenfoldSF(Q,args): #big & ugly function - do later (if ever)
576#    pass
577   
578def InterPrecipitateSF(Q,args):
579    ''' Computes structure factor for precipitates in a matrix
580    Refs.: E-Wen Huang, Peter K. Liaw, Lionel Porcar, Yun Liu, Yee-Lang Liu,
581    Ji-Jung Kai, and Wei-Ren Chen,APPLIED PHYSICS LETTERS 93, 161904 (2008)
582    R. Giordano, A. Grasso, and J. Teixeira, Phys. Rev. A 43, 6894    1991   
583    param float Q: Q value array (A-1)
584    param array args: [float R, float VolFr]: "radius" & volume fraction
585    returns numpy array S(Q)
586    '''
587    R,VolFr = args
588    QV2 = Q**2*VolFr**2
589    top = 1 - np.exp(-QV2/4)*np.cos(2.*Q*R)
590    bot = 1-2*np.exp(-QV2/4)*np.cos(2.*Q*R) + np.exp(-QV2/2)
591    return 2*(top/bot) - 1
592         
593################################################################################
594##### SB-MaxEnt
595################################################################################
596
597def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
598    '''Calculates the response matrix :math:`G(Q,r)`
599   
600    :param float q: :math:`Q`
601    :param float r: :math:`r`
602    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
603    :param function FFfxn: form factor function FF(q,r,args)
604    :param function Volfxn: volume function Vol(r,args)
605    :returns float: G(Q,r)
606    '''
607    FF = FFfxn(q,r,args)
608    Vol = Volfxn(r,args)
609    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
610   
611'''
612sbmaxent
613
614Entropy maximization routine as described in the article
615J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
616("MNRAS": "Monthly Notices of the Royal Astronomical Society")
617
618:license: Copyright (c) 2013, UChicago Argonne, LLC
619:license: This file is distributed subject to a Software License Agreement found
620     in the file LICENSE that is included with this distribution.
621
622References:
623
6241. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
6252. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
626   Neutron Scattering Data Analysis, Rutherford
627   Appleton Laboratory, UK, 1986; ed. MW Johnson,
628   IOP Conference Series 81 (1986) 81 - 86, Institute
629   of Physics, Bristol, UK.
6303. ID Culverwell and GP Clarke; Ibid. 87 - 96.
6314. JA Potton, GK Daniell, & BD Rainford,
632   J APPL CRYST 21 (1988) 663 - 668.
6335. JA Potton, GJ Daniell, & BD Rainford,
634   J APPL CRYST 21 (1988) 891 - 897.
635
636'''
637
638class MaxEntException(Exception): 
639    '''Any exception from this module'''
640    pass
641
642def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
643    '''
644    do the complete Maximum Entropy algorithm of Skilling and Bryan
645   
646    :param float datum[]:
647    :param float sigma[]:
648    :param float[][] G: transformation matrix
649    :param float base[]:
650    :param int IterMax:
651    :param obj image_to_data: opus function (defaults to opus)
652    :param obj data_to_image: tropus function (defaults to tropus)
653   
654    :returns float[]: :math:`f(r) dr`
655    '''
656       
657    TEST_LIMIT        = 0.05                    # for convergence
658    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
659    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
660    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
661    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
662   
663    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
664    MOVE_PASSES       = 0.001                   # convergence test in routine: move
665
666    def tropus (data, G):
667        '''
668        tropus: transform data-space -> solution-space:  [G] * data
669       
670        default definition, caller can use this definition or provide an alternative
671       
672        :param float[M] data: observations, ndarray of shape (M)
673        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
674        :returns float[N]: calculated image, ndarray of shape (N)
675        '''
676        return G.dot(data)
677
678    def opus (image, G):
679        '''
680        opus: transform solution-space -> data-space:  [G]^tr * image
681       
682        default definition, caller can use this definition or provide an alternative
683       
684        :param float[N] image: solution, ndarray of shape (N)
685        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
686        :returns float[M]: calculated data, ndarray of shape (M)
687        '''
688        return np.dot(G.T,image)    #G.transpose().dot(image)
689
690    def Dist(s2, beta):
691        '''measure the distance of this possible solution'''
692        w = 0
693        n = beta.shape[0]
694        for k in range(n):
695            z = -sum(s2[k] * beta)
696            w += beta[k] * z
697        return w
698   
699    def ChiNow(ax, c1, c2, s1, s2):
700        '''
701        ChiNow
702       
703        :returns tuple: (ChiNow computation of ``w``, beta)
704        '''
705       
706        bx = 1 - ax
707        a =   bx * c2  -  ax * s2
708        b = -(bx * c1  -  ax * s1)
709   
710        beta = ChoSol(a, b)
711        w = 1.0
712        for k in range(SEARCH_DIRECTIONS):
713            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
714        return w, beta
715   
716    def ChoSol(a, b):
717        '''
718        ChoSol: ? chop the solution vectors ?
719       
720        :returns: new vector beta
721        '''
722        n = b.shape[0]
723        fl = np.zeros((n,n))
724        bl = np.zeros_like(b)
725       
726        #print_arr("ChoSol: a", a)
727        #print_vec("ChoSol: b", b)
728   
729        if (a[0][0] <= 0):
730            msg = "ChoSol: a[0][0] = " 
731            msg += str(a[0][0])
732            msg += '  Value must be positive'
733            raise MaxEntException(msg)
734   
735        # first, compute fl from a
736        # note fl is a lower triangular matrix
737        fl[0][0] = math.sqrt (a[0][0])
738        for i in (1, 2):
739            fl[i][0] = a[i][0] / fl[0][0]
740            for j in range(1, i+1):
741                z = 0.0
742                for k in range(j):
743                    z += fl[i][k] * fl[j][k]
744                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
745                z = a[i][j] - z
746                if j == i:
747                    y = math.sqrt(z)
748                else:
749                    y = z / fl[j][j]
750                fl[i][j] = y
751        #print_arr("ChoSol: fl", fl)
752   
753        # next, compute bl from fl and b
754        bl[0] = b[0] / fl[0][0]
755        for i in (1, 2):
756            z = 0.0
757            for k in range(i):
758                z += fl[i][k] * bl[k]
759                #print "\t", i, k, z
760            bl[i] = (b[i] - z) / fl[i][i]
761        #print_vec("ChoSol: bl", bl)
762   
763        # last, compute beta from bl and fl
764        beta = np.empty((n))
765        beta[-1] = bl[-1] / fl[-1][-1]
766        for i in (1, 0):
767            z = 0.0
768            for k in range(i+1, n):
769                z += fl[k][i] * beta[k]
770                #print "\t\t", i, k, 'z=', z
771            beta[i] = (bl[i] - z) / fl[i][i]
772        #print_vec("ChoSol: beta", beta)
773   
774        return beta
775
776    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
777        '''
778        move beta one step closer towards the solution
779        '''
780        a_lower, a_upper = 0., 1.          # bracket  "a"
781        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
782        #print "MaxEntMove: cmin = %g" % cmin
783        if cmin*chisq > chizer:
784            ctarg = (1.0 + cmin)/2
785        else:
786            ctarg = chizer/chisq
787        f_lower = cmin - ctarg
788        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
789        f_upper = c_upper - ctarg
790   
791        fx = 2*MOVE_PASSES      # just to start off
792        loop = 1
793        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
794            a_new = (a_lower + a_upper) * 0.5           # search by bisection
795            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
796            fx = c_new - ctarg
797            # tighten the search range for the next pass
798            if f_lower*fx > 0:
799                a_lower, f_lower = a_new, fx
800            if f_upper*fx > 0:
801                a_upper, f_upper = a_new, fx
802            loop += 1
803   
804        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
805            msg = "MaxEntMove: Loop counter = " 
806            msg += str(MAX_MOVE_LOOPS)
807            msg += '  No convergence in alpha chop'
808            raise MaxEntException(msg)
809   
810        w = Dist (s2, beta);
811        m = SEARCH_DIRECTIONS
812        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
813            for k in range(m):
814                beta[k] *= math.sqrt (fSum/(blank*w))
815        chtarg = ctarg * chisq
816        return w, chtarg, loop, a_new, fx, beta
817       
818#MaxEnt_SB starts here   
819
820    if image_to_data == None:
821        image_to_data = opus
822    if data_to_image == None:
823        data_to_image = tropus
824    n   = len(base)
825    npt = len(datum)
826
827    # Note that the order of subscripts for
828    # "xi" and "eta" has been reversed from
829    # the convention used in the FORTRAN version
830    # to enable parts of them to be passed as
831    # as vectors to "image_to_data" and "data_to_image".
832    xi      = np.zeros((SEARCH_DIRECTIONS, n))
833    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
834    beta    = np.zeros((SEARCH_DIRECTIONS))
835    # s1      = np.zeros((SEARCH_DIRECTIONS))
836    # c1      = np.zeros((SEARCH_DIRECTIONS))
837    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
838    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
839
840    # TODO: replace blank (scalar) with base (vector)
841    blank = sum(base) / len(base)   # use the average value of base
842
843    chizer, chtarg = npt*1.0, npt*1.0
844    f = base * 1.0                              # starting distribution is base
845
846    fSum  = sum(f)                              # find the sum of the f-vector
847    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
848    chisq = sum(z*z)                            # Chi^2, SB eq. 4
849
850    for iter in range(IterMax):
851        ox = -2 * z / sigma                        # gradient of Chi^2
852
853        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
854        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
855        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
856        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
857        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
858
859        a = 1.0
860        b = 1.0 / cnorm
861        if iter == 0:
862            test = 0.0     # mismatch between entropy and ChiSquared gradients
863        else:
864            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
865            a = 0.5 / (snorm * test)
866            b *= 0.5 / test
867        xi[0] = f * cgrad / cnorm
868        xi[1] = f * (a * sgrad - b * cgrad)
869
870        eta[0] = image_to_data (xi[0], G);          # image --> data
871        eta[1] = image_to_data (xi[1], G);          # image --> data
872        ox = eta[1] / (sigma * sigma)
873        xi[2] = data_to_image (ox, G);              # data --> image
874        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
875        xi[2] = f * xi[2] * a
876        eta[2] = image_to_data (xi[2], G)           # image --> data
877       
878#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
879#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
880
881        # prepare the search directions for the conjugate gradient technique
882        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
883        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
884#         print_vec("MaxEnt: c1", c1)
885#         print_vec("MaxEnt: s1", s1)
886
887        for k in range(SEARCH_DIRECTIONS):
888            for l in range(k+1):
889                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
890                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
891#         print_arr("MaxEnt: c2", c2)
892#         print_arr("MaxEnt: s2", s2)
893
894        # reflect across the body diagonal
895        for k, l in ((0,1), (0,2), (1,2)):
896            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
897            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
898 
899        beta[0] = -0.5 * c1[0] / c2[0][0]
900        beta[1] = 0.0
901        beta[2] = 0.0
902        if (iter > 0):
903            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
904
905        f_old = f.copy()    # preserve the last image
906        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
907       
908        # As mentioned at the top of p.119,
909        # need to protect against stray negative values.
910        # In this case, set them to RESET_STRAYS * base[i]
911        #f = f.clip(RESET_STRAYS * blank, f.max())
912        for i in range(n):
913            if f[i] <= 0.0:
914                f[i] = RESET_STRAYS * base[i]
915        df = f - f_old
916        fSum = sum(f)
917        fChange = sum(df)
918
919        # calculate the normalized entropy
920        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
921        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
922        chisq = sum(z*z)                            # report this ChiSq
923
924        if report:
925            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
926            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
927            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
928
929        # See if we have finished our task.
930        # do the hardest test first
931        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
932            print ' Convergence achieved.'
933            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
934    print ' No convergence! Try increasing Error multiplier.'
935    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
936
937   
938###############################################################################
939#### IPG/TNNLS Routines
940###############################################################################
941
942def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
943    ''' An implementation of the Interior-Point Gradient method of
944    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
945    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
946    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
947    Problem addressed: Total Non-Negative Least Squares (TNNLS)
948    :param float datum[]:
949    :param float sigma[]:
950    :param float[][] G: transformation matrix
951    :param int IterMax:
952    :param float Qvec: data positions for Power = 0-4
953    :param float approach: 0.8 default fitting parameter
954    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
955   
956    '''
957    if Power < 0: 
958        GmatE = G/sigma[:np.newaxis]
959        IntE = datum/sigma
960        pwr = 0
961        QvecP = np.ones_like(datum)
962    else:
963        GmatE = G[:]
964        IntE = datum[:]
965        pwr = Power
966        QvecP = Qvec**pwr
967    Amat = GmatE*QvecP[:np.newaxis]
968    AAmat = np.inner(Amat,Amat)
969    Bvec = IntE*QvecP
970    Xw = np.ones_like(Bins)*1.e-32
971    calc = np.dot(G.T,Xw)
972    nIter = 0
973    err = 10.
974    while (nIter<IterMax) and (err > 1.):
975        #Step 1 in M&Z paper:
976        Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
977        Dk = Xw/np.inner(AAmat,Xw)
978        Pk = -Dk*Qk
979        #Step 2 in M&Z paper:
980        alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
981        alpWv = -Xw/Pk
982        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
983        #Step 3 in M&Z paper:
984        shift = AkSt*Pk
985        Xw += shift
986        #done IPG; now results
987        nIter += 1
988        calc = np.dot(G.T,Xw)
989        chisq = np.sum(((datum-calc)/sigma)**2)
990        err = chisq/len(datum)
991        if report:
992            print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
993    return chisq,Xw,calc
994
995###############################################################################
996#### SASD Utilities
997###############################################################################
998
999def SetScale(Data,refData):
1000    Profile,Limits,Sample = Data
1001    refProfile,refLimits,refSample = refData
1002    x,y = Profile[:2]
1003    rx,ry = refProfile[:2]
1004    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
1005    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
1006    iBeg = np.searchsorted(x,Beg)
1007    iFin = np.searchsorted(x,Fin)+1        #include last point
1008    sum = np.sum(y[iBeg:iFin])
1009    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
1010    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
1011   
1012def Bestimate(G,Rg,P):
1013    return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))
1014   
1015def SmearData(Ic,Q,slitLen,Back):
1016    Np = Q.shape[0]
1017    Qtemp = np.concatenate([Q,Q[-1]+20*Q])
1018    Ictemp = np.concatenate([Ic,Ic[-1]*(1-(Qtemp[Np:]-Qtemp[Np])/(20*Qtemp[Np-1]))])
1019    Icsm = np.zeros_like(Q)
1020    Qsm = 2*slitLen*(np.interp(np.arange(2*Np)/2.,np.arange(Np),Q)-Q[0])/(Q[-1]-Q[0])
1021    Sp = np.searchsorted(Qsm,slitLen)
1022    DQsm = np.diff(Qsm)[:Sp]
1023    Ism = np.interp(np.sqrt(Q[:,np.newaxis]**2+Qsm**2),Qtemp,Ictemp)
1024    Icsm = np.sum((Ism[:,:Sp]*DQsm),axis=1)
1025    Icsm /= slitLen
1026    return Icsm
1027   
1028###############################################################################
1029#### Size distribution
1030###############################################################################
1031
1032def SizeDistribution(Profile,ProfDict,Limits,Sample,data):
1033    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1034        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1035        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1036        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1037        'Cylinder diam':[CylinderDFF,CylinderDVol],
1038        'Spherical shell': [SphericalShellFF, SphericalShellVol]}
1039    Shape = data['Size']['Shape'][0]
1040    Parms = data['Size']['Shape'][1:]
1041    if data['Size']['logBins']:
1042        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
1043            data['Size']['Nbins']+1,True)/2.        #make radii
1044    else:
1045        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
1046            data['Size']['Nbins']+1,True)/2.        #make radii
1047    Dbins = np.diff(Bins)
1048    Bins = Bins[:-1]+Dbins/2.
1049    Contrast = Sample['Contrast'][1]
1050    Scale = Sample['Scale'][0]
1051    Sky = 10**data['Size']['MaxEnt']['Sky']
1052    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
1053    Back = data['Back']
1054    Q,Io,wt,Ic,Ib = Profile[:5]
1055    Qmin = Limits[1][0]
1056    Qmax = Limits[1][1]
1057    wtFactor = ProfDict['wtFactor']
1058    Ibeg = np.searchsorted(Q,Qmin)
1059    Ifin = np.searchsorted(Q,Qmax)+1        #include last point
1060    BinMag = np.zeros_like(Bins)
1061    Ic[:] = 0.
1062    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
1063    if 'MaxEnt' == data['Size']['Method']:
1064        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
1065            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
1066            data['Size']['MaxEnt']['Niter'],report=True)
1067    elif 'IPG' == data['Size']['Method']:
1068        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
1069            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
1070            Power=data['Size']['IPG']['Power'],report=True)
1071    Ib[:] = Back[0]
1072    Ic[Ibeg:Ifin] += Back[0]
1073    print ' Final chi^2: %.3f'%(chisq)
1074    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
1075       
1076################################################################################
1077#### Modelling
1078################################################################################
1079
1080def ModelFit(Profile,ProfDict,Limits,Sample,Model):
1081    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1082        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1083        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1084        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1085        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol],
1086        'Spherical shell':[SphericalShellFF,SphericalShellVol]}
1087           
1088    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
1089            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
1090               
1091    parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
1092        'PkInt','PkPos','PkSig','PkGam',]
1093       
1094    FFparmOrder = ['Aspect ratio','Length','Diameter','Thickness','Shell thickness']
1095   
1096    SFparmOrder = ['Dist','VolFr','epis','Sticky','Depth','Width']
1097
1098    def GetModelParms():
1099        parmDict = {'Scale':Sample['Scale'][0],'SlitLen':Sample.get('SlitLen',0.0),}
1100        varyList = []
1101        values = []
1102        levelTypes = []
1103        Back = Model['Back']
1104        if Back[1]:
1105            varyList += ['Back',]
1106            values.append(Back[0])
1107        parmDict['Back'] = Back[0]
1108        partData = Model['Particle']
1109        for i,level in enumerate(partData['Levels']):
1110            cid = str(i)+';'
1111            controls = level['Controls']
1112            Type = controls['DistType']
1113            levelTypes.append(Type)
1114            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1115                if Type not in ['Monodosperse',]:
1116                    parmDict[cid+'NumPoints'] = controls['NumPoints']
1117                    parmDict[cid+'Cutoff'] = controls['Cutoff']
1118                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
1119                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
1120                parmDict[cid+'StrFact'] = sfxns[controls['StrFact']]
1121                parmDict[cid+'Contrast'] = controls['Contrast']
1122                for item in FFparmOrder:
1123                    if item in controls['FFargs']:
1124                        parmDict[cid+item] = controls['FFargs'][item][0]
1125                        if controls['FFargs'][item][1]:
1126                            varyList.append(cid+item)
1127                            values.append(controls['FFargs'][item][0])
1128                for item in SFparmOrder:
1129                    if item in controls.get('SFargs',{}):
1130                        parmDict[cid+item] = controls['SFargs'][item][0]
1131                        if controls['SFargs'][item][1]:
1132                            varyList.append(cid+item)
1133                            values.append(controls['SFargs'][item][0])
1134            distDict = controls['DistType']
1135            for item in parmOrder:
1136                if item in level[distDict]:
1137                    parmDict[cid+item] = level[distDict][item][0]
1138                    if level[distDict][item][1]:
1139                        values.append(level[distDict][item][0])
1140                        varyList.append(cid+item)
1141        return levelTypes,parmDict,varyList,values
1142       
1143    def SetModelParms():
1144        print ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
1145        if 'Back' in varyList:
1146            Model['Back'][0] = parmDict['Back']
1147            print %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back'])
1148        partData = Model['Particle']
1149        for i,level in enumerate(partData['Levels']):
1150            controls = level['Controls']
1151            Type = controls['DistType']
1152            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1153                print ' Component %d: Type: %s: Structure Factor: %s Contrast: %12.3f'  \
1154                    %(i,Type,controls['StrFact'],controls['Contrast'])               
1155            else:
1156                print ' Component %d: Type: %s: '%(i,Type,)
1157            cid = str(i)+';'
1158            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1159                for item in FFparmOrder:
1160                    if cid+item in varyList:
1161                        controls['FFargs'][item][0] = parmDict[cid+item]
1162                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1163                for item in SFparmOrder:
1164                    if cid+item in varyList:
1165                        controls['SFargs'][item][0] = parmDict[cid+item]
1166                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1167            distDict = controls['DistType']
1168            for item in level[distDict]:
1169                if cid+item in varyList:
1170                    level[distDict][item][0] = parmDict[cid+item]
1171                    print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1172                   
1173    def calcSASD(values,Q,Io,wt,Ifb,levelTypes,parmDict,varyList):
1174        parmDict.update(zip(varyList,values))
1175        M = np.sqrt(wt)*(getSASD(Q,levelTypes,parmDict)+Ifb-parmDict['Scale']*Io)
1176        return M
1177       
1178    def getSASD(Q,levelTypes,parmDict):
1179        Ic = np.zeros_like(Q)
1180        for i,Type in enumerate(levelTypes):
1181            cid = str(i)+';'
1182            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1183                FFfxn = parmDict[cid+'FormFact']
1184                Volfxn = parmDict[cid+'FFVolume']
1185                SFfxn = parmDict[cid+'StrFact']
1186                FFargs = []
1187                SFargs = []
1188                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',cid+'Shell thickness']:
1189                    if item in parmDict: 
1190                        FFargs.append(parmDict[item])
1191                for item in [cid+'Dist',cid+'VolFr',cid+'epis',cid+'Sticky',cid+'Depth',cid+'Width']:
1192                    if item in parmDict: 
1193                        SFargs.append(parmDict[item])
1194                distDict = {}
1195                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
1196                    if item in parmDict:
1197                        distDict[item.split(';')[1]] = parmDict[item]
1198                contrast = parmDict[cid+'Contrast']
1199                rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)
1200                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
1201                dist *= parmDict[cid+'Volume']
1202                Ic += np.dot(Gmat,dist)*SFfxn(Q,args=SFargs)
1203            elif 'Unified' in Type:
1204                Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
1205                    parmDict[cid+'P'],parmDict[cid+'Cutoff']
1206                Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
1207                Guin = G*np.exp(-(Q*Rg)**2/3)
1208                Porod = (B/Qstar**P)*np.exp(-(Q*Rgco)**2/3)
1209                Ic += Guin+Porod
1210            elif 'Porod' in Type:
1211                B,P,Rgco = parmDict[cid+'B'],parmDict[cid+'P'],parmDict[cid+'Cutoff']
1212                Porod = (B/Q**P)*np.exp(-(Q*Rgco)**2/3)
1213                Ic += Porod
1214            elif 'Mono' in Type:
1215                FFfxn = parmDict[cid+'FormFact']
1216                Volfxn = parmDict[cid+'FFVolume']
1217                SFfxn = parmDict[cid+'StrFact']
1218                FFargs = []
1219                SFargs = []
1220                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',cid+'Shell thickness']:
1221                    if item in parmDict: 
1222                        FFargs.append(parmDict[item])
1223                for item in [cid+'Dist',cid+'VolFr',cid+'epis',cid+'Sticky',cid+'Depth',cid+'Width']:
1224                    if item in parmDict: 
1225                        SFargs.append(parmDict[item])
1226                contrast = parmDict[cid+'Contrast']
1227                R = parmDict[cid+'Radius']
1228                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
1229                Ic += Gmat[0]*parmDict[cid+'Volume']*SFfxn(Q,args=SFargs)
1230            elif 'Bragg' in Type:
1231                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
1232                    parmDict[cid+'PkSig'],parmDict[cid+'PkGam'],Q)
1233        Ic += parmDict['Back']  #/parmDict['Scale']
1234        slitLen = Sample['SlitLen']
1235        if slitLen:
1236            Ic = SmearData(Ic,Q,slitLen,parmDict['Back'])
1237        return Ic
1238       
1239    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1240    Qmin = Limits[1][0]
1241    Qmax = Limits[1][1]
1242    wtFactor = ProfDict['wtFactor']
1243    Ibeg = np.searchsorted(Q,Qmin)
1244    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
1245    Ic[:] = 0
1246    levelTypes,parmDict,varyList,values = GetModelParms()
1247    if varyList:
1248        result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
1249            args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList))
1250        parmDict.update(zip(varyList,result[0]))
1251        chisq = np.sum(result[2]['fvec']**2)
1252        ncalc = result[2]['nfev']
1253        covM = result[1]
1254    else:   #nothing varied
1255        M = calcSASD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList)
1256        chisq = np.sum(M**2)
1257        ncalc = 0
1258        covM = []
1259        sig = []
1260        sigDict = {}
1261        result = []
1262    Rvals = {}
1263    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
1264    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
1265    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],levelTypes,parmDict)
1266    Msg = 'Failed to converge'
1267    try:
1268        Nans = np.isnan(result[0])
1269        if np.any(Nans):
1270            name = varyList[Nans.nonzero(True)[0]]
1271            Msg = 'Nan result for '+name+'!'
1272            raise ValueError
1273        Negs = np.less_equal(result[0],0.)
1274        if np.any(Negs):
1275            name = varyList[Negs.nonzero(True)[0]]
1276            Msg = 'negative coefficient for '+name+'!'
1277            raise ValueError
1278        if len(covM):
1279            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
1280            sigDict = dict(zip(varyList,sig))
1281        print ' Results of small angle data modelling fit:'
1282        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
1283        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1284        SetModelParms()
1285        covMatrix = covM*Rvals['GOF']
1286        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
1287    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
1288        return False,0,0,0,0,0,0,Msg
1289   
1290def ModelFxn(Profile,ProfDict,Limits,Sample,sasdData):
1291   
1292    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1293        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1294        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1295        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1296        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol],
1297        'Spherical shell':[SphericalShellFF,SphericalShellVol]}
1298    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
1299            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
1300
1301#    pdb.set_trace()
1302    partData = sasdData['Particle']
1303    Back = sasdData['Back']
1304    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1305    Qmin = Limits[1][0]
1306    Qmax = Limits[1][1]
1307    Ibeg = np.searchsorted(Q,Qmin)
1308    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
1309    Ib[:] = Back[0]
1310    Ic[:] = 0
1311    Rbins = []
1312    Dist = []
1313    for level in partData['Levels']:
1314        controls = level['Controls']
1315        distFxn = controls['DistType']
1316        if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1317            parmDict = level[controls['DistType']]
1318            FFfxn = shapes[controls['FormFact']][0]
1319            Volfxn = shapes[controls['FormFact']][1]
1320            SFfxn = sfxns[controls['StrFact']]
1321            FFargs = []
1322            SFargs = []
1323            for item in ['Dist','VolFr','epis','Sticky','Depth','Width',]:
1324                if item in controls.get('SFargs',{}):
1325                    SFargs.append(controls['SFargs'][item][0])
1326            for item in ['Aspect ratio','Length','Thickness','Diameter','Shell thickness']:
1327                if item in controls['FFargs']: 
1328                    FFargs.append(controls['FFargs'][item][0])
1329            contrast = controls['Contrast']
1330            distDict = {}
1331            for item in parmDict:
1332                distDict[item] = parmDict[item][0]
1333            rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],distDict)
1334            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
1335            dist *= level[distFxn]['Volume'][0]
1336            Ic[Ibeg:Ifin] += np.dot(Gmat,dist)*SFfxn(Q[Ibeg:Ifin],args=SFargs)
1337            Rbins.append(rBins)
1338            Dist.append(dist/(4.*dBins))
1339        elif 'Unified' in distFxn:
1340            parmDict = level[controls['DistType']]
1341            Rg,G,B,P,Rgco = parmDict['Rg'][0],parmDict['G'][0],parmDict['B'][0],    \
1342                parmDict['P'][0],parmDict['Cutoff'][0]
1343            Qstar = Q[Ibeg:Ifin]/(scsp.erf(Q[Ibeg:Ifin]*Rg/np.sqrt(6)))**3
1344            Guin = G*np.exp(-(Q[Ibeg:Ifin]*Rg)**2/3)
1345            Porod = (B/Qstar**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1346            Ic[Ibeg:Ifin] += Guin+Porod
1347            Rbins.append([])
1348            Dist.append([])
1349        elif 'Porod' in distFxn:
1350            parmDict = level[controls['DistType']]
1351            B,P,Rgco = parmDict['B'][0],parmDict['P'][0],parmDict['Cutoff'][0]
1352            Porod = (B/Q[Ibeg:Ifin]**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1353            Ic[Ibeg:Ifin] += Porod
1354            Rbins.append([])
1355            Dist.append([])
1356        elif 'Mono' in distFxn:
1357            parmDict = level[controls['DistType']]
1358            R = level[controls['DistType']]['Radius'][0]
1359            FFfxn = shapes[controls['FormFact']][0]
1360            Volfxn = shapes[controls['FormFact']][1]
1361            SFfxn = sfxns[controls['StrFact']]
1362            FFargs = []
1363            SFargs = []
1364            for item in ['Dist','VolFr','epis','Sticky','Depth','Width',]:
1365                if item in controls.get('SFargs',{}):
1366                    SFargs.append(controls['SFargs'][item][0])
1367            for item in ['Aspect ratio','Length','Thickness','Diameter','Shell thickness']:
1368                if item in controls['FFargs']: 
1369                    FFargs.append(controls['FFargs'][item][0])
1370            contrast = controls['Contrast']
1371            Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
1372            Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]*SFfxn(Q[Ibeg:Ifin],args=SFargs)
1373            Rbins.append([])
1374            Dist.append([])
1375        elif 'Bragg' in distFxn:
1376            parmDict = level[controls['DistType']]
1377            Ic[Ibeg:Ifin] += parmDict['PkInt'][0]*G2pwd.getPsVoigt(parmDict['PkPos'][0],
1378                parmDict['PkSig'][0],parmDict['PkGam'][0],Q[Ibeg:Ifin])
1379            Rbins.append([])
1380            Dist.append([])
1381    Ic[Ibeg:Ifin] += Back[0]
1382    slitLen = Sample.get('SlitLen',0.)
1383    if slitLen:
1384        Ic[Ibeg:Ifin] = SmearData(Ic,Q,slitLen,Back[0])[Ibeg:Ifin]
1385    sasdData['Size Calc'] = [Rbins,Dist]
1386   
1387def MakeDiamDist(DistName,nPoints,cutoff,distDict):
1388   
1389    if 'LogNormal' in DistName:
1390        distFxn = 'LogNormalDist'
1391        cumeFxn = 'LogNormalCume'
1392        pos = distDict['MinSize']
1393        args = [distDict['Mean'],distDict['StdDev']]
1394        step = 4.*np.sqrt(np.exp(distDict['StdDev']**2)*(np.exp(distDict['StdDev']**2)-1.))
1395        mode = distDict['MinSize']+distDict['Mean']/np.exp(distDict['StdDev']**2)
1396        minX = 1. #pos
1397        maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
1398    elif 'Gauss' in DistName:
1399        distFxn = 'GaussDist'
1400        cumeFxn = 'GaussCume'
1401        pos = distDict['Mean']
1402        args = [distDict['StdDev']]
1403        step = 0.02*distDict['StdDev']
1404        mode = distDict['Mean']
1405        minX = np.max([mode-4.*distDict['StdDev'],1.])
1406        maxX = np.min([mode+4.*distDict['StdDev'],1.e5])
1407    elif 'LSW' in DistName:
1408        distFxn = 'LSWDist'
1409        cumeFxn = 'LSWCume'
1410        pos = distDict['Mean']
1411        args = []
1412        minX,maxX = [0.,2.*pos]
1413    elif 'Schulz' in DistName:
1414        distFxn = 'SchulzZimmDist'
1415        cumeFxn = 'SchulzZimmCume'
1416        pos = distDict['Mean']
1417        args = [distDict['StdDev']]
1418        minX = np.max([1.,pos-4.*distDict['StdDev']])
1419        maxX = np.min([pos+4.*distDict['StdDev'],1.e5])
1420    nP = 500
1421    Diam = np.logspace(0.,5.,nP,True)
1422    TCW = eval(cumeFxn+'(Diam,pos,args)')
1423    CumeTgts = np.linspace(cutoff,(1.-cutoff),nPoints+1,True)
1424    rBins = np.interp(CumeTgts,TCW,Diam,0,0)
1425    dBins = np.diff(rBins)
1426    rBins = rBins[:-1]+dBins/2.
1427    return rBins,dBins,eval(distFxn+'(rBins,pos,args)')
1428
1429################################################################################
1430#### MaxEnt testing stuff
1431################################################################################
1432
1433def print_vec(text, a):
1434    '''print the contents of a vector to the console'''
1435    n = a.shape[0]
1436    print "%s[ = (" % text,
1437    for i in range(n):
1438        s = " %g, " % a[i]
1439        print s,
1440    print ")"
1441
1442def print_arr(text, a):
1443    '''print the contents of an array to the console'''
1444    n, m = a.shape
1445    print "%s[][] = (" % text
1446    for i in range(n):
1447        print " (",
1448        for j in range(m):
1449            print " %g, " % a[i][j],
1450        print "),"
1451    print ")"
1452
1453def test_MaxEnt_SB(report=True):
1454    def readTextData(filename):
1455        '''return q, I, dI from a 3-column text file'''
1456        if not os.path.exists(filename):
1457            raise Exception("file not found: " + filename)
1458        buf = [line.split() for line in open(filename, 'r').readlines()]
1459        buf = zip(*buf)         # transpose rows and columns
1460        q  = np.array(buf[0], dtype=np.float64)
1461        I  = np.array(buf[1], dtype=np.float64)
1462        dI = np.array(buf[2], dtype=np.float64)
1463        return q, I, dI
1464    print "MaxEnt_SB: "
1465    test_data_file = os.path.join( 'testinp', 'test.sas')
1466    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
1467    bkg   = 0.1     #   I = I - bkg
1468    dMin, dMax, nRadii = 25, 9000, 40
1469    defaultDistLevel = 1.0e-6
1470    IterMax = 40
1471    errFac = 1.05
1472   
1473    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
1474    dr   = r * (r[1]/r[0] - 1)          # step size
1475    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
1476    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
1477   
1478    qVec, I, dI = readTextData(test_data_file)
1479    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
1480   
1481    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
1482    if f_dr is None:
1483        print "no solution"
1484        return
1485   
1486    print "solution reached"
1487    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
1488        print '%10.4f %10.4f %12.4g'%(a,b,c)
1489
1490def tests():
1491    test_MaxEnt_SB(report=True)
1492
1493if __name__ == '__main__':
1494    tests()
1495
Note: See TracBrowser for help on using the repository browser.