source: trunk/GSASIIsasd.py @ 2308

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

add sp.shell model or SASD - Thx to L.A, Avakyan
fix (again) page switching problem with Texture
Change import string for G2img files

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