source: trunk/GSASIIsasd.py @ 3298

Last change on this file since 3298 was 3136, checked in by vondreele, 8 years ago

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