source: trunk/GSASIIcomp.py @ 49

Last change on this file since 49 was 49, checked in by vondreel, 13 years ago

resolved zoom issues with peak picking, etc.
further steps for image integration
File size: 56.5 KB
Line 
1#GSASII computational module
2import sys
3import math
4import wx
5import time
6import numpy as np
7import numpy.linalg as nl
8import os.path as ospath
9# determine a binary path pased on the host OS and the python version, path is relative to
10# location of this file
11if sys.platform == "win32":
12    bindir = 'binwin%d.%d' % sys.version_info[0:2]
13elif sys.platform == "darwin":
14    bindir = 'binmac%d.%d' % sys.version_info[0:2]
15else:
16    bindir = 'bin'
17
18if ospath.exists(ospath.join(sys.path[0],bindir)) and ospath.join(sys.path[0],bindir) not in sys.path: 
19    sys.path.insert(0,ospath.join(sys.path[0],bindir))
20
21try: 
22    import pypowder as pyp
23except:
24    # create an app to display the error, since we are still loading routines at this point
25    app = wx.App()
26    app.MainLoop()
27    msg = wx.MessageDialog(None, message="Unable to load the GSAS powder computation module, pypowder",
28                     caption="Import Error",
29                     style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
30    msg.ShowModal()
31    # this error is non-recoverable, so just quit
32    exit()
33import GSASIIplot as G2plt
34
35# trig functions in degrees
36sind = lambda x: math.sin(x*math.pi/180.)
37asind = lambda x: 180.*math.asin(x)/math.pi
38tand = lambda x: math.tan(x*math.pi/180.)
39atand = lambda x: 180.*math.atan(x)/math.pi
40atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
41cosd = lambda x: math.cos(x*math.pi/180.)
42acosd = lambda x: 180.*math.acos(x)/math.pi
43rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
44#numpy versions
45npsind = lambda x: np.sin(x*np.pi/180.)
46npasind = lambda x: 180.*np.arcsin(x)/math.pi
47npcosd = lambda x: np.cos(x*math.pi/180.)
48nptand = lambda x: np.tan(x*math.pi/180.)
49npatand = lambda x: 180.*np.arctan(x)/np.pi
50npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
51
52def sec2HMS(sec):
53    H = int(sec/3600)
54    M = int(sec/60-H*60)
55    S = sec-3600*H-60*M
56    return '%d:%2d:%.2f'%(H,M,S)
57
58#def valueEsd(value,esd,precision):
59#    nDecimals = lambda esd: max(0.,1.545-math.log10(abs(esd)))
60#    if esd:
61#    else:
62#       
63#       
64
65def DoPeakFit(peaks,background,limits,inst,data):
66   
67    def backgroundPrint(background,sigback):
68        if background[1]:
69            print 'Background coefficients for',background[0],'function'
70            ptfmt = "%12.5f"
71            ptstr =  'values:'
72            sigstr = 'esds  :'
73            for i,back in enumerate(background[3:]):
74                ptstr += ptfmt % (back)
75                sigstr += ptfmt % (sigback[i+3])
76            print ptstr
77            print sigstr
78        else:
79            print 'Background not refined'
80   
81    def instPrint(instVal,siginst,insLabels):
82        print 'Instrument Parameters:'
83        ptfmt = "%12.6f"
84        ptlbls = 'names :'
85        ptstr =  'values:'
86        sigstr = 'esds  :'
87        for i,value in enumerate(instVal):
88            ptlbls += "%s" % (insLabels[i].center(12))
89            ptstr += ptfmt % (value)
90            if siginst[i]:
91                sigstr += ptfmt % (siginst[i])
92            else:
93                sigstr += 12*' '
94        print ptlbls
95        print ptstr
96        print sigstr
97   
98    def peaksPrint(peaks,sigpeaks):
99        print 'Peak list='
100
101    begin = time.time()
102    Np = len(peaks[0])
103    DataType = inst[1][0]
104    instVal = inst[1][1:]
105    insref = inst[2][1:]
106    insLabels = inst[3][1:]
107    Ka2 = False
108    Ioff = 3
109    if len(instVal) == 11:
110        lamratio = instVal[1]/instVal[0]
111        Ka2 = True
112        Ioff = 5
113    insref = insref[len(insref)-6:]
114    for peak in peaks:
115        dip = []
116        dip.append(tand(peak[0]/2.0)**2)
117        dip.append(tand(peak[0]/2.0))
118        dip.append(1.0/cosd(peak[0]/2.0))
119        dip.append(tand(peak[0]/2.0))
120        peak.append(dip)
121    B = background[2]
122    bcof = background[3:3+B]
123    Bv = 0
124    if background[1]:
125        Bv = B
126    x,y,w,yc,yb,yd = data               #these are numpy arrays!
127    V = []
128    A = []
129    swobs = 0.0
130    smin = 0.0
131    first = True
132    GoOn = True
133    Go = True
134    dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \
135        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
136    screenSize = wx.DisplaySize()
137    Size = dlg.GetSize()
138    dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
139    try:
140        i = 0
141        for xi in x:
142            Go = dlg.Update(i)[0]
143            if GoOn:
144                GoOn = Go
145            if limits[0] <= xi <= limits[1]:
146                yb[i] = 0.0
147                dp = []
148                for j in range(B):
149                    t = (xi-limits[0])**j
150                    yb[i] += t*bcof[j]
151                    if background[1]:
152                        dp.append(t)
153                yc[i] = yb[i]
154                Iv = 0
155                for j in range(6):
156                    if insref[j]:
157                        dp.append(0.0)
158                        Iv += 1
159                for peak in peaks:
160                    dip = peak[-1]
161                    f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],peak[8],0.0)
162                    yc[i] += f[0]*peak[2]
163                    if f[0] > 0.0:
164                        j = 0
165                        if insref[0]:              #U
166                            dp[Bv+j] += f[3]*dip[0]
167                            j += 1
168                        if insref[1]:              #V
169                            dp[Bv+j] += f[3]*dip[1]
170                            j += 1
171                        if insref[2]:              #W
172                            dp[Bv+j] += f[3]
173                            j += 1
174                        if insref[3]:              #X
175                            dp[Bv+j] += f[4]*dip[2]
176                            j += 1
177                        if insref[4]:              #Y
178                            dp[Bv+j] += f[4]*dip[3]
179                            j += 1
180                        if insref[5]:              #SH/L
181                            dp[Bv+j] += f[5]
182                    if Ka2:
183                       pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0))
184                       f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],peak[8],0.0)
185                       yc[i] += f2[0]*peak[2]*instVal[3]
186                       if f[0] > 0.0:
187                           j = 0
188                           if insref[0]:              #U
189                               dp[Bv+j] += f2[3]*dip[0]*instVal[3]
190                               j += 1
191                           if insref[1]:              #V
192                               dp[Bv+j] += f2[3]*dip[1]*instVal[3]
193                               j += 1
194                           if insref[2]:              #W
195                               dp[Bv+j] += f2[3]*instVal[3]
196                               j += 1
197                           if insref[3]:              #X
198                               dp[Bv+j] += f2[4]*dip[2]*instVal[3]
199                               j += 1
200                           if insref[4]:              #Y
201                               dp[Bv+j] += f2[4]*dip[3]*instVal[3]
202                               j += 1
203                           if insref[5]:              #SH/L
204                               dp[Bv+j] += f2[5]*instVal[3]                       
205                    for j in range(0,Np,2):          #14s
206                        if peak[j+1]: dp.append(f[j/2+1])
207                yd[i] = y[i]-yc[i]
208                swobs += w[i]*y[i]**2
209                t2 = w[i]*yd[i]
210                smin += t2*yd[i]                 
211                if first:
212                    first = False
213                    M = len(dp)
214                    A = np.zeros(shape=(M,M))
215                    V = np.zeros(shape=(M))
216                A,V = pyp.buildmv(t2,w[i],M,dp,A,V)
217            i += 1
218    finally:
219        dlg.Destroy()
220    print GoOn
221    Rwp = smin/swobs
222    Rwp = math.sqrt(Rwp)*100.0
223    norm = np.diag(A)
224    for i,elm in enumerate(norm):
225        if elm <= 0.0:
226            print norm
227            return False,0,0,0
228    for i in xrange(len(V)):
229        norm[i] = 1.0/math.sqrt(norm[i])
230        V[i] *= norm[i]
231        a = A[i]
232        for j in xrange(len(V)):
233            a[j] *= norm[i]
234    A = np.transpose(A)
235    for i in xrange(len(V)):
236        a = A[i]
237        for j in xrange(len(V)):
238            a[j] *= norm[i]
239    b = nl.solve(A,V)
240    A = nl.inv(A)
241    sig = np.diag(A)
242    for i in xrange(len(V)):
243        b[i] *= norm[i]
244        sig[i] *= norm[i]*norm[i]
245        sig[i] = math.sqrt(abs(sig[i]))
246    sigback = [0,0,0]
247    for j in range(Bv):
248        background[j+3] += b[j]
249        sigback.append(sig[j])
250    backgroundPrint(background,sigback)
251    k = 0
252    delt = []
253    siginst = []
254    for i in range(len(instVal)-6):
255        siginst.append(0.0)
256    for j in range(6):
257        if insref[j]:
258            instVal[j+Ioff] += b[Bv+k]*0.5
259            siginst.append(sig[Bv+k])
260            delt.append(b[Bv+k]*0.5)
261            k += 1
262        else:
263            delt.append(0.0)
264            siginst.append(0.0)
265    instPrint(instVal,siginst,insLabels)
266    inst[1] = [DataType,]
267    for val in instVal:
268        inst[1].append(val)
269    B = Bv+Iv
270    for peak in peaks:
271        del peak[-1]                        # remove dip from end
272        delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2]
273        delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0)
274        delshl = delt[5]   
275        for j in range(0,len(peak[:-1]),2):
276            if peak[j+1]: 
277                peak[j] += b[B]*0.5
278                B += 1
279        peak[4] += delsig
280        if peak[4] < 0.0:
281            print 'ERROR - negative sigma'
282            return False,0,0,0,False           
283        peak[6] += delgam
284        if peak[6] < 0.0:
285            print 'ERROR - negative gamma'
286            return False,0,0,0,False
287        peak[8] += delshl           
288        if peak[8] < 0.0:
289            print 'ERROR - negative SH/L'
290            return False,0,0,0,False           
291    runtime = time.time()-begin   
292    data = [x,y,w,yc,yb,yd]
293    return True,smin,Rwp,runtime,GoOn
294   
295#some cell utilities
296#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
297#           cell - a,b,c,alp,bet,gam in A & deg
298   
299def calc_rDsq(H,A):
300    rdsq = H[0]*H[0]*A[0]+H[1]*H[1]*A[1]+H[2]*H[2]*A[2]+H[0]*H[1]*A[3]+H[0]*H[2]*A[4]+H[1]*H[2]*A[5]
301    return rdsq
302   
303def calc_rDsqZ(H,A,Z,tth,lam):
304    rpd = math.pi/180.
305    rdsq = calc_rDsq(H,A)+Z*math.sin(tth*rpd)*2.0*rpd/(lam*lam)
306    return rdsq
307   
308def calc_rVsq(A):
309    rVsq = A[0]*A[1]*A[2]+0.25*(A[3]*A[4]*A[5]-A[0]*A[5]**2-A[1]*A[4]**2-A[2]*A[3]**2)
310    if rVsq < 0:
311        return 1
312    return rVsq
313   
314def calc_rV(A):
315    return math.sqrt(calc_rVsq(A))
316   
317def calc_V(A):
318    return 1./calc_rV(A)
319   
320def scaleAbyV(A,V):
321    v = calc_V(A)
322    scale = math.exp(math.log(v/V)/3.)**2
323    for i in range(6):
324        A[i] *= scale
325   
326def ranaxis(dmin,dmax):
327    import random as rand
328    return rand.random()*(dmax-dmin)+dmin
329   
330def ran2axis(k,N):
331    import random as rand
332    T = 1.5+0.49*k/N
333#    B = 0.99-0.49*k/N
334#    B = 0.99-0.049*k/N
335    B = 0.99-0.149*k/N
336    R = (T-B)*rand.random()+B
337    return R
338   
339def ranNaxis(k,N):
340    import random as rand
341    T = 1.0+1.0*k/N
342    B = 1.0-1.0*k/N
343    R = (T-B)*rand.random()+B
344    return R
345   
346def ranAbyV(Bravais,dmin,dmax,V):
347    cell = [0,0,0,0,0,0]
348    bad = True
349    while bad:
350        bad = False
351        cell = rancell(Bravais,dmin,dmax)
352        G,g = cell2Gmat(cell)
353        A = Gmat2A(G)
354        if calc_rVsq(A) < 1:
355            scaleAbyV(A,V)
356            cell = A2cell(A)
357            for i in range(3):
358                bad |= cell[i] < dmin
359    return A
360   
361def ranAbyR(Bravais,A,k,N,ranFunc):
362    R = ranFunc(k,N)
363    if Bravais in [0,1,2]:          #cubic - not used
364        A[0] = A[1] = A[2] = A[0]*R
365        A[3] = A[4] = A[5] = 0.
366    elif Bravais in [3,4]:          #hexagonal/trigonal
367        A[0] = A[1] = A[3] = A[0]*R
368        A[2] *= R
369        A[4] = A[5] = 0.       
370    elif Bravais in [5,6]:          #tetragonal
371        A[0] = A[1] = A[0]*R
372        A[2] *= R
373        A[3] = A[4] = A[5] = 0.       
374    elif Bravais in [7,8,9,10]:     #orthorhombic
375        A[0] *= R
376        A[1] *= R
377        A[2] *= R
378        A[3] = A[4] = A[5] = 0.       
379    elif Bravais in [11,12]:        #monoclinic
380        A[0] *= R
381        A[1] *= R
382        A[2] *= R
383        A[4] *= R
384        A[3] = A[5] = 0.       
385    else:                           #triclinic
386        A[0] *= R
387        A[1] *= R
388        A[2] *= R
389        A[3] *= R
390        A[4] *= R
391        A[5] *= R
392    return A
393   
394def rancell(Bravais,dmin,dmax):
395    if Bravais in [0,1,2]:          #cubic
396        a = b = c = ranaxis(dmin,dmax)
397        alp = bet = gam = 90
398    elif Bravais in [3,4]:          #hexagonal/trigonal
399        a = b = ranaxis(dmin,dmax)
400        c = ranaxis(dmin,dmax)
401        alp = bet =  90
402        gam = 120
403    elif Bravais in [5,6]:          #tetragonal
404        a = b = ranaxis(dmin,dmax)
405        c = ranaxis(dmin,dmax)
406        alp = bet = gam = 90
407    elif Bravais in [7,8,9,10]:       #orthorhombic - F,I,C,P - a<b<c convention
408        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
409        abc.sort()
410        a = abc[0]
411        b = abc[1]
412        c = abc[2]
413        alp = bet = gam = 90
414    elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
415        ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
416        ac.sort()
417        a = ac[0]
418        b = ranaxis(dmin,dmax)
419        c = ac[1]
420        alp = gam = 90
421        bet = ranaxis(90.,130.)
422    else:                           #triclinic - a<b<c convention
423        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
424        abc.sort()
425        a = abc[0]
426        b = abc[1]
427        c = abc[2]
428        r = 0.5*b/c
429        alp = ranaxis(acosd(r),acosd(-r))
430        r = 0.5*a/c
431        bet = ranaxis(acosd(r),acosd(-r))
432        r = 0.5*a/b
433        gam = ranaxis(acosd(r),acosd(-r)) 
434    return [a,b,c,alp,bet,gam]
435   
436def A2Gmat(A):
437    G = np.zeros(shape=(3,3))
438    G = [[A[0],A[3]/2.,A[4]/2.], [A[3]/2.,A[1],A[5]/2.], [A[4]/2.,A[5]/2.,A[2]]]
439    g = nl.inv(G)
440    return G,g
441   
442def fillgmat(cell):
443    a,b,c,alp,bet,gam = cell
444    g = np.array([[a*a,a*b*cosd(gam),a*c*cosd(bet)],[a*b*cosd(gam),b*b,b*c*cosd(alp)],
445        [a*c*cosd(bet),b*c*cosd(alp),c*c]])
446    return g
447           
448def cell2Gmat(cell):
449    #returns reciprocal (G) & real (g) metric tensors
450    g = fillgmat(cell)
451    G = nl.inv(g)       
452    return G,g
453   
454def invcell2Gmat(invcell):
455    G = fillgmat(invcell)
456    g = nl.inv(G)
457    return G,g
458   
459def Gmat2cell(g):
460    #returns lattice parameters from real metric tensor (g)
461    a = math.sqrt(max(0,g[0][0]))
462    b = math.sqrt(max(0,g[1][1]))
463    c = math.sqrt(max(0,g[2][2]))
464    alp = acosd(g[2][1]/(b*c))
465    bet = acosd(g[2][0]/(a*c))
466    gam = acosd(g[0][1]/(a*b))
467    return a,b,c,alp,bet,gam
468   
469def Gmat2A(G):
470    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
471   
472def cell2A(cell):
473    G,g = cell2Gmat(cell)
474    return Gmat2A(G)
475   
476def A2cell(A):
477    G,g = A2Gmat(A)
478    return Gmat2cell(g)
479   
480def A2invcell(A):
481    ainv = math.sqrt(max(0.,A[0]))
482    binv = math.sqrt(max(0.,A[1]))
483    cinv = math.sqrt(max(0.,A[2]))
484    gaminv = acosd(max(-0.5,min(0.5,0.5*A[3]/(ainv*binv))))
485    betinv = acosd(max(-0.5,min(0.5,0.5*A[4]/(ainv*cinv))))
486    alpinv = acosd(max(-0.5,min(0.5,0.5*A[5]/(binv*cinv))))
487    return ainv,binv,cinv,alpinv,betinv,gaminv
488   
489def cell2AB(cell):
490    #from real lattice parameters - cell
491    # returns A for Cartesian to crystal transformations A*X = x
492    # and inverse B for crystal to Cartesian transformation B*x = X
493    G,g = cell2Gmat(cell)       #reciprocal & real metric tensors
494    cosAlpStar = G[2][1]/math.sqrt(G[1][1]*G[2][2])
495    sinAlpStar = math.sqrt(1.0-cosAlpStar**2)
496    B = np.eye(3)
497    B *= cell[:3]
498    A = np.zeros(shape=(3,3))
499    A[0][0] = 1.0
500    A[0][1] = cosd(cell[5])
501    A[1][1] = sinAlpStar*sind(cell[5])
502    A[1][2] = -cosAlpStar*sind(cell[5])
503    A[0][2] = cosd(cell[4])
504    A[2][2] = sind(cell[4])
505    B = np.dot(A,B)
506    A = nl.inv(B)
507    return A,B
508   
509def makeMat(Angle,Axis):
510    #Make rotation matrix from Angle in degrees,Axis =0 for rotation about x, =1 for about y, etc.
511    cs = cosd(Angle)
512    ss = sind(Angle)
513    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
514    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
515                   
516def MaxIndex(dmin,A):
517    #finds maximum allowed hkl for given A within dmin
518    Hmax = [0,0,0]
519    try:
520        cell = A2cell(A)
521    except:
522        cell = [1,1,1,90,90,90]
523    for i in range(3):
524        Hmax[i] = int(round(cell[i]/dmin))
525    return Hmax
526   
527def sortHKLd(HKLd,ifreverse,ifdup):
528    #HKLd is a list of [h,k,l,d,...]; ifreverse=True for largest d first
529    #ifdup = True if duplicate d-spacings allowed
530    T = []
531    for i,H in enumerate(HKLd):
532        if ifdup:
533            T.append((H[3],i))
534        else:
535            T.append(H[3])           
536    D = dict(zip(T,HKLd))
537    T.sort()
538    if ifreverse:
539        T.reverse()
540    X = []
541    okey = ''
542    for key in T: 
543        if key != okey: X.append(D[key])    #remove duplicate d-spacings
544        okey = key
545    return X
546   
547def sortM20(cells):
548    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
549    #sort highest M20 1st
550    T = []
551    for i,M in enumerate(cells):
552        T.append((M[0],i))
553    D = dict(zip(T,cells))
554    T.sort()
555    T.reverse()
556    X = []
557    for key in T:
558        X.append(D[key])
559    return X
560               
561 
562def GenHBravais(dmin,Bravais,A):
563# dmin - minimum d-spacing
564# Bravais in range(14) to indicate Bravais lattice; 0-2 cubic, 3,4 - hexagonal/trigonal,
565# 5,6 - tetragonal, 7-10 - orthorhombic, 11,12 - monoclinic, 13 - triclinic
566# A - as defined in calc_rDsq
567# returns HKL = [h,k,l,d,0] sorted so d largest first
568    import math
569    if Bravais in [9,11]:
570        Cent = 'C'
571    elif Bravais in [1,5,8]:
572        Cent = 'I'
573    elif Bravais in [0,7]:
574        Cent = 'F'
575    elif Bravais in [3]:
576        Cent = 'R'
577    else:
578        Cent = 'P'
579    Hmax = MaxIndex(dmin,A)
580    dminsq = 1./(dmin**2)
581    HKL = []
582    if Bravais == 13:                       #triclinic
583        for l in range(-Hmax[2],Hmax[2]+1):
584            for k in range(-Hmax[1],Hmax[1]+1):
585                hmin = 0
586                if (k < 0): hmin = 1
587                if (k ==0 and l < 0): hmin = 1
588                for h in range(hmin,Hmax[0]+1):
589                    H=[h,k,l]
590                    rdsq = calc_rDsq(H,A)
591                    if 0 < rdsq <= dminsq:
592                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
593    elif Bravais in [11,12]:                #monoclinic - b unique
594        Hmax = SwapIndx(2,Hmax)
595        for h in range(Hmax[0]+1):
596            for k in range(-Hmax[1],Hmax[1]+1):
597                lmin = 0
598                if k < 0:lmin = 1
599                for l in range(lmin,Hmax[2]+1):
600                    [h,k,l] = SwapIndx(-2,[h,k,l])
601                    H = []
602                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
603                    if H:
604                        rdsq = calc_rDsq(H,A)
605                        if 0 < rdsq <= dminsq:
606                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
607                    [h,k,l] = SwapIndx(2,[h,k,l])
608    elif Bravais in [7,8,9,10]:            #orthorhombic
609        for h in range(Hmax[0]+1):
610            for k in range(Hmax[1]+1):
611                for l in range(Hmax[2]+1):
612                    H = []
613                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
614                    if H:
615                        rdsq = calc_rDsq(H,A)
616                        if 0 < rdsq <= dminsq:
617                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
618    elif Bravais in [5,6]:                  #tetragonal
619        for l in range(Hmax[2]+1):
620            for k in range(Hmax[1]+1):
621                for h in range(k,Hmax[0]+1):
622                    H = []
623                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
624                    if H:
625                        rdsq = calc_rDsq(H,A)
626                        if 0 < rdsq <= dminsq:
627                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
628    elif Bravais in [3,4]:
629        lmin = 0
630        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
631        for l in range(lmin,Hmax[2]+1):
632            for k in range(Hmax[1]+1):
633                hmin = k
634                if l < 0: hmin += 1
635                for h in range(hmin,Hmax[0]+1):
636                    H = []
637                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
638                    if H:
639                        rdsq = calc_rDsq(H,A)
640                        if 0 < rdsq <= dminsq:
641                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
642
643    else:                                   #cubic
644        for l in range(Hmax[2]+1):
645            for k in range(l,Hmax[1]+1):
646                for h in range(k,Hmax[0]+1):
647                    H = []
648                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
649                    if H:
650                        rdsq = calc_rDsq(H,A)
651                        if 0 < rdsq <= dminsq:
652                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
653    return sortHKLd(HKL,True,False)
654   
655def SwapXY(x,y):
656    return [y,x]
657   
658def SwapIndx(Axis,H):
659    if Axis in [1,-1]:
660        return H
661    elif Axis in [2,-3]:
662        return [H[1],H[2],H[0]]
663    else:
664        return [H[2],H[0],H[1]]
665       
666def CentCheck(Cent,H):
667    h,k,l = H
668    if Cent == 'A' and (k+l)%2:
669        return False
670    elif Cent == 'B' and (h+l)%2:
671        return False
672    elif Cent == 'C' and (h+k)%2:
673        return False
674    elif Cent == 'I' and (h+k+l)%2:
675        return False
676    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
677        return False
678    elif Cent == 'R' and (-h+k+l)%3:
679        return False
680    else:
681        return True
682                                   
683def GenHLaue(dmin,Laue,Cent,Axis,A):
684# dmin - minimum d-spacing
685# Laue - Laue group symbol = '-1','2/m','mmmm','4/m','6/m','4/mmm','6/mmm',
686#                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
687# Cent - lattice centering = 'P','A','B','C','I','F'
688# Axis - code for unique monoclinic axis = 'a','b','c'
689# A - 6 terms as defined in calc_rDsq
690# returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
691# part of reciprocal space ignoring anomalous dispersion
692    import math
693    Hmax = MaxIndex(dmin,A)
694    dminsq = 1./(dmin**2)
695    HKL = []
696    if Laue == '-1':                       #triclinic
697        for l in range(-Hmax[2],Hmax[2]+1):
698            for k in range(-Hmax[1],Hmax[1]+1):
699                hmin = 0
700                if (k < 0) or (k ==0 and l < 0): hmin = 1
701                for h in range(hmin,Hmax[0]+1):
702                    H = []
703                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
704                    rdsq = calc_rDsq(H,A)
705                    if 0 < rdsq <= dminsq:
706                        HKL.append([h,k,l,1/math.sqrt(rdsq)])
707    elif Laue == '2/m':                #monoclinic
708        Hmax = SwapIndx(Axis,Hmax)
709        for h in range(Hmax[0]+1):
710            for k in range(-Hmax[1],Hmax[1]+1):
711                lmin = 0
712                if k < 0:lmin = 1
713                for l in range(lmin,Hmax[2]+1):
714                    [h,k,l] = SwapIndx(-Axis,[h,k,l])
715                    H = []
716                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
717                    if H:
718                        rdsq = calc_rDsq(H,A)
719                        if 0 < rdsq <= dminsq:
720                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
721                    [h,k,l] = SwapIndx(Axis,[h,k,l])
722    elif Laue in ['mmm','4/m','6/m']:            #orthorhombic
723        for l in range(Hmax[2]+1):
724            for h in range(Hmax[0]+1):
725                kmin = 1
726                if Laue == 'mmm' or h ==0: kmin = 0
727                for k in range(kmin,Hmax[2]+1):
728                    H = []
729                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
730                    if H:
731                        rdsq = calc_rDsq(H,A)
732                        if 0 < rdsq <= dminsq:
733                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
734    elif Laue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
735        for l in range(Hmax[2]+1):
736            for h in range(Hmax[0]+1):
737                for k in range(h+1):
738                    H = []
739                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
740                    if H:
741                        rdsq = calc_rDsq(H,A)
742                        if 0 < rdsq <= dminsq:
743                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
744    elif Laue in ['3m1','31m','3','3R','3mR']:                  #trigonals
745        for l in range(-Hmax[2],Hmax[2]+1):
746            hmin = 0
747            if l < 0: hmin = 1
748            for h in range(hmin,Hmax[0]+1):
749                if Laue in ['3R','3']:
750                    kmax = h
751                    kmin = -int((h-1)/2.)
752                else:
753                    kmin = 0
754                    kmax = h
755                    if Laue in ['3m1','3mR'] and l < 0: kmax = h-1
756                    if Laue == '31m' and l < 0: kmin = 1
757                for k in range(kmin,kmax+1):
758                    H = []
759                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
760                    if H:
761                        rdsq = calc_rDsq(H,A)
762                        if 0 < rdsq <= dminsq:
763                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
764    else:                                   #cubic
765        for h in range(Hmax[0]+1):
766            for k in range(h+1):
767                lmin = 0
768                lmax = k
769                if Laue =='m3':
770                    lmax = h-1
771                    if h == k: lmax += 1
772                for l in range(lmin,lmax+1):
773                    H = []
774                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
775                    if H:
776                        rdsq = calc_rDsq(H,A)
777                        if 0 < rdsq <= dminsq:
778                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
779    return sortHKLd(HKL,True,True)
780   
781def calc_M20(peaks,HKL):
782    diff = 0
783    X20 = 0
784    for Nobs20,peak in enumerate(peaks):
785        if peak[3]:
786            Qobs = 1.0/peak[7]**2
787            Qcalc = 1.0/peak[8]**2
788            diff += abs(Qobs-Qcalc)
789        elif peak[2]:
790            X20 += 1
791        if Nobs20 == 19: 
792            d20 = peak[7]
793            break
794    else:
795        d20 = peak[7]
796        Nobs20 = len(peaks)
797    for N20,hkl in enumerate(HKL):
798        if hkl[3] < d20:
799            break               
800    eta = diff/Nobs20
801    Q20 = 1.0/d20**2
802    if diff:
803        M20 = Q20/(2.0*diff)
804    else:
805        M20 = 0
806    M20 /= (1.+X20)
807    return M20,X20
808   
809def IndexPeaks(peaks,HKL):
810    import bisect
811    hklds = [1000.0]                                    #make bounded list of available d-spacings
812    N = len(HKL)
813    if N == 0: return False
814    for hkl in HKL:
815        hklds.append(hkl[3])
816    hklds.append(0.0)
817    hklds.sort()                                        # ascending sort - upper bound at end
818    hklmax = [0,0,0]
819    for ipk,peak in enumerate(peaks):
820        if peak[2]:
821            i = bisect.bisect_right(hklds,peak[7])          # find peak position in hkl list
822            dm = peak[7]-hklds[i-1]                         # peak to neighbor hkls in list
823            dp = hklds[i]-peak[7]
824            pos = N-i                                       # reverse the order
825            if dp > dm: pos += 1                            # closer to upper than lower
826            hkl = HKL[pos]                                 # put in hkl
827            if hkl[4] >= 0:                                 # peak already assigned - test if this one better
828                opeak = peaks[hkl[4]]
829                dold = abs(opeak[7]-hkl[3])
830                dnew = min(dm,dp)
831                if dold > dnew:                             # new better - zero out old
832                    for j in range(3):
833                        opeak[j+4] = 0
834                    opeak[8] = 0.
835                else:                                       # old better - do nothing
836                    continue               
837            hkl[4] = ipk
838            for j in range(3):
839                peak[j+4] = hkl[j]
840            peak[8] = hkl[3]                                # fill in d-calc
841    for peak in peaks:
842        peak[3] = False
843        if peak[2]:
844            if peak[8] > 0.:
845                for j in range(3):
846                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
847                peak[3] = True
848    if hklmax[0]*hklmax[1]*hklmax[2] > 0:
849        return True
850    else:
851        return False
852   
853def FitHKL(ibrav,peaks,A,wtP):
854    def ShiftTest(a,b):
855        if b < -0.1*a: 
856            b = -0.0001*a
857        return b
858    smin = 0.
859    first = True
860    for peak in peaks:
861        if peak[2] and peak[3]:
862            h,k,l = H = peak[4:7]
863            Qo = 1./peak[7]**2
864            Qc = calc_rDsq(H,A)
865            try:
866                peak[8] = 1./math.sqrt(Qc)
867            except:
868                print A2invcell(A)
869            delt = Qo-Qc
870            smin += delt**2
871            dp = []
872            if ibrav in [0,1,2]:            #m3m
873                dp.append(h*h+k*k+l*l)
874            elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
875                dp.append(h*h+k*k+h*k)
876                dp.append(l*l)
877            elif ibrav in [5,6]:            #4/mmm
878                dp.append(h*h+k*k)
879                dp.append(l*l)
880            elif ibrav in [7,8,9,10]:       #mmm
881                dp.append(h*h)
882                dp.append(k*k)
883                dp.append(l*l)
884            elif ibrav in [11,12]:          #2/m
885                dp.append(h*h)
886                dp.append(k*k)
887                dp.append(l*l)
888                dp.append(h*l)
889            else:                           #1
890#    derivatives for a0*h^2+a1*k^2+a2*l^2+a3*h*k+a4*h*l+a5*k*l
891                dp.append(h*h)
892                dp.append(k*k)
893                dp.append(l*l)
894                dp.append(h*k)
895                dp.append(h*l)
896                dp.append(k*l)
897            if first:
898                first = False
899                M = len(dp)
900                B = np.zeros(shape=(M,M))
901                V = np.zeros(shape=(M))
902            dwt = peak[7]**wtP
903            B,V = pyp.buildmv(delt*dwt,dwt,M,dp,B,V)
904    if nl.det(B) > 0.0:
905        try:
906            b = nl.solve(B,V)
907            B = nl.inv(B)
908            sig = np.diag(B)
909        except SingularMatrix:
910            return False,0
911        if ibrav in [0,1,2]:            #m3m
912            A[0] += ShiftTest(A[0],b[0])
913            A[1] = A[2] = A[0]
914        elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
915            A[0] += ShiftTest(A[0],b[0])
916            A[1] = A[3] = A[0]
917            A[2] += ShiftTest(A[2],b[1])
918        elif ibrav in [5,6]:            #4/mmm
919            A[0] += ShiftTest(A[0],b[0])
920            A[1] = A[0]
921            A[2] += ShiftTest(A[2],b[1])
922        elif ibrav in [7,8,9,10]:       #mmm
923            for i in range(3):
924                A[i] += ShiftTest(A[i],b[i])
925        elif ibrav in [11,12]:          #2/m
926            for i in range(3):
927                A[i] += ShiftTest(A[i],b[i])
928            A[4] += ShiftTest(A[4],b[3])
929            A[4] = min(1.4*math.sqrt(A[0]*A[2]),A[4])   #min beta star = 45
930        else:                           #1
931            oldV = math.sqrt(1./calc_rVsq(A))
932            oldA = A[:]
933            for i in range(6):
934                A[i] += b[i]*0.2
935            A[3] = min(1.1*math.sqrt(max(0,A[1]*A[2])),A[3])
936            A[4] = min(1.1*math.sqrt(max(0,A[0]*A[2])),A[4])
937            A[5] = min(1.1*math.sqrt(max(0,A[0]*A[1])),A[5])
938            ratio = math.sqrt(1./calc_rVsq(A))/oldV
939            if 0.9 > ratio or ratio > 1.1:
940                A = oldA
941#    else:
942#        print B
943#        print V
944#        for peak in peaks:
945#            print peak
946    return True,smin
947       
948def FitHKLZ(ibrav,peaks,Z,A):
949    return A,Z
950   
951def rotOrthoA(A):
952    return [A[1],A[2],A[0],0,0,0]
953   
954def swapMonoA(A):
955    return [A[2],A[1],A[0],0,A[4],0]
956   
957def oddPeak(indx,peaks):
958    noOdd = True
959    for peak in peaks:
960        H = peak[4:7]
961        if H[indx] % 2:
962            noOdd = False
963    return noOdd
964   
965def halfCell(ibrav,A,peaks):
966    if ibrav in [0,1,2]:
967        if oddPeak(0,peaks):
968            A[0] *= 2
969            A[1] = A[2] = A[0]
970    elif ibrav in [3,4,5,6]:
971        if oddPeak(0,peaks):
972            A[0] *= 2
973            A[1] = A[0]
974        if oddPeak(2,peaks):
975            A[2] *=2
976    else:
977        if oddPeak(0,peaks):
978            A[0] *=2
979        if oddPeak(1,peaks):
980            A[1] *=2
981        if oddPeak(2,peaks):
982            A[2] *=2
983    return A
984   
985def getDmin(peaks):
986    return peaks[-1][7]
987   
988def getDmax(peaks):
989    return peaks[0][7]
990   
991def refinePeaks(peaks,ibrav,A):
992    dmin = getDmin(peaks)
993    smin = 1.0e10
994    pwr = 6
995    maxTries = 10
996    if ibrav == 13:
997        pwr = 4
998        maxTries = 10
999    OK = False
1000    tries = 0
1001    HKL = GenHBravais(dmin,ibrav,A)
1002    while IndexPeaks(peaks,HKL):
1003        Pwr = pwr - 2*(tries % 2)
1004        HKL = []
1005        tries += 1
1006        osmin = smin
1007        oldA = A
1008        OK,smin = FitHKL(ibrav,peaks,A,Pwr)
1009        for a in A[:3]:
1010            if a < 0:
1011                A = oldA
1012                OK = False
1013                break
1014        if OK:
1015            HKL = GenHBravais(dmin,ibrav,A)
1016        if len(HKL) == 0: break                         #absurd cell obtained!
1017        rat = (osmin-smin)/smin
1018        if abs(rat) < 1.0e-5 or not OK: break
1019        if tries > maxTries: break
1020    if OK:
1021        OK,smin = FitHKL(ibrav,peaks,A,4)
1022    M20,X20 = calc_M20(peaks,HKL)
1023    return len(HKL),M20,X20
1024       
1025def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
1026# dlg & ncMax are used for wx progress bar
1027# A != 0 find the best A near input A,
1028# A = 0 for random cell, volume normalized to V1;
1029# returns number of generated hkls, M20, X20 & A for best found
1030    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
1031    dmin = getDmin(peaks)-0.05
1032    amin = 2.5
1033    amax = 5.*getDmax(peaks)
1034    Asave = []
1035    GoOn = True
1036    if A:
1037        HKL = GenHBravais(dmin,ibrav,A[:])
1038        if len(HKL) > mHKL[ibrav]:
1039            IndexPeaks(peaks,HKL)
1040            Asave.append([calc_M20(peaks,HKL),A[:]])
1041    tries = 0
1042    while tries < Ntries:
1043        if A:
1044            Anew = ranAbyR(ibrav,A[:],tries+1,Ntries,ran2axis)
1045            if ibrav in [11,12,13]:
1046                Anew = ranAbyR(ibrav,A[:],tries/10+1,Ntries,ran2axis)
1047        else:
1048            Anew = ranAbyV(ibrav,amin,amax,V1)
1049        HKL = GenHBravais(dmin,ibrav,Anew)
1050        if len(HKL) > mHKL[ibrav] and IndexPeaks(peaks,HKL):
1051            Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
1052            Asave.append([calc_M20(peaks,HKL),Anew[:]])
1053            if ibrav == 9:                          #C-centered orthorhombic
1054                for i in range(2):
1055                    Anew = rotOrthoA(Anew[:])
1056                    Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
1057                    HKL = GenHBravais(dmin,ibrav,Anew)
1058                    IndexPeaks(peaks,HKL)
1059                    Asave.append([calc_M20(peaks,HKL),Anew[:]])
1060            elif ibrav == 11:                      #C-centered monoclinic
1061                Anew = swapMonoA(Anew[:])
1062                Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
1063                HKL = GenHBravais(dmin,ibrav,Anew)
1064                IndexPeaks(peaks,HKL)
1065                Asave.append([calc_M20(peaks,HKL),Anew[:]])
1066        else:
1067            break
1068        Nc = len(HKL)
1069        if Nc >= ncMax:
1070            GoOn = False
1071        elif dlg:
1072            GoOn = dlg.Update(Nc)[0]
1073            if not GoOn:
1074                break
1075        tries += 1   
1076    X = sortM20(Asave)
1077    if X:
1078        Lhkl,M20,X20 = refinePeaks(peaks,ibrav,X[0][1])
1079        return GoOn,Lhkl,M20,X20,X[0][1]
1080    else:
1081        return GoOn,0,0,0,0
1082       
1083def monoCellReduce(ibrav,A):
1084    a,b,c,alp,bet,gam = A2cell(A)
1085    G,g = A2Gmat(A)
1086    if ibrav in [11]:
1087        u = [0,0,-1]
1088        v = [1,0,2]
1089        anew = math.sqrt(np.dot(np.dot(v,g),v))
1090        if anew < a:
1091            cang = np.dot(np.dot(u,g),v)/(anew*c)
1092            beta = acosd(-abs(cang))
1093            A = cell2A([anew,b,c,90,beta,90])
1094    else:
1095        u = [-1,0,0]
1096        v = [1,0,1]
1097        cnew = math.sqrt(np.dot(np.dot(v,g),v))
1098        if cnew < c:
1099            cang = np.dot(np.dot(u,g),v)/(a*cnew)
1100            beta = acosd(-abs(cang))
1101            A = cell2A([a,b,cnew,90,beta,90])
1102    return A
1103
1104def DoIndexPeaks(peaks,inst,controls,bravais):
1105   
1106    def peakDspace(peaks,A):
1107        for peak in peaks:
1108            if peak[3]:
1109                dsq = calc_rDsq(peak[4:7],A)
1110                if dsq > 0:
1111                    peak[8] = 1./math.sqrt(dsq)
1112        return
1113    delt = 0.05                                     #lowest d-spacing cushion - can be fixed?
1114    amin = 2.5
1115    amax = 5.0*getDmax(peaks)
1116    dmin = getDmin(peaks)-delt
1117    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
1118        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
1119        'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
1120    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
1121    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
1122    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
1123    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
1124    Nobs = len(peaks)
1125    wave = inst[1]
1126    if len(inst) > 10:
1127        zero = inst[3]
1128    else:
1129        zero = inst[2]
1130    print "%s %8.5f %6.3f" % ('wavelength, zero =',wave,zero)
1131    print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax)
1132    ifzero,maxzero,ncno = controls[:3]
1133    ncMax = Nobs*ncno
1134    print "%s %d %s %d %s %d" % ('change zero =',ifzero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs)
1135    cells = []
1136    for ibrav in range(14):
1137        begin = time.time()
1138        if bravais[ibrav]:
1139            print 'cell search for ',bravaisNames[ibrav]
1140            print '      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test'
1141            V1 = controls[3]
1142            bestM20 = 0
1143            topM20 = 0
1144            cycle = 0
1145            while cycle < 5:
1146                dlg = wx.ProgressDialog("Generated reflections",tries[cycle]+" cell search for "+bravaisNames[ibrav],ncMax, \
1147                    style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
1148                screenSize = wx.DisplaySize()
1149                Size = dlg.GetSize()
1150                dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
1151                try:
1152                    GoOn = True
1153                    while GoOn:                                                 #Loop over increment of volume
1154                        N2 = 0
1155                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
1156                            if ibrav > 2:
1157                                if not N2:
1158                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
1159                                if A:
1160                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
1161                            else:
1162                                GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
1163                            if Nc >= ncMax:
1164                                GoOn = False
1165                                break
1166                            elif 3*Nc < Nobs:
1167                                N2 = 10
1168                                break
1169                            else:
1170                                if not GoOn:
1171                                    break
1172                                if M20 > 1.0:
1173                                    bestM20 = max(bestM20,M20)
1174                                    A = halfCell(ibrav,A[:],peaks)
1175                                    if ibrav in [12]:
1176                                        A = monoCellReduce(ibrav,A[:])
1177                                    HKL = GenHBravais(dmin,ibrav,A)
1178                                    IndexPeaks(peaks,HKL)
1179                                    a,b,c,alp,bet,gam = A2cell(A)
1180                                    V = calc_V(A)
1181                                    print "%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f" % (M20,X20,Nc,a,b,c,alp,bet,gam,V,V1)
1182                                    if M20 >= 2.0:
1183                                        cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False])
1184                            if not GoOn:
1185                                break
1186                            N2 += 1
1187                        if ibrav < 11:
1188                            V1 *= 1.1
1189                        elif ibrav in range(11,14):
1190                            V1 *= 1.05
1191                        if not GoOn:
1192                            if bestM20 > topM20:
1193                                topM20 = bestM20
1194                                if cells:
1195                                    V1 = cells[0][9]
1196                                else:
1197                                    V1 = 25
1198                                ncMax += Nobs
1199                                cycle += 1
1200                                print 'Restart search, new Max Nc = ',ncMax
1201                            else:
1202                                cycle = 10
1203                finally:
1204                    dlg.Destroy()
1205            print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
1206                ', elapsed time = ',sec2HMS(time.time()-begin))
1207           
1208    if cells:
1209        cells = sortM20(cells)
1210        cells[0][-1] = True
1211        return True,dmin,cells
1212    else:
1213        return False,0,0
1214       
1215def FitRing(ring):
1216    Err,parms = FitCircle(ring)
1217    Err /= len(ring)
1218#    print 'circle error:','%8f'%(Err)
1219    if Err > 20000.:
1220        eparms = FitEllipse(ring)
1221        if eparms:
1222            parms = eparms
1223    return parms
1224       
1225def FitCircle(ring):
1226    import numpy.linalg as nl
1227   
1228    def makeParmsCircle(B):
1229        cent = [-B[0]/2,-B[1]/2]
1230        phi = 0.
1231        sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2])
1232        return cent,phi,[sr1,sr2]
1233       
1234    ring = np.array(ring)
1235    x = np.asarray(ring.T[0])
1236    y = np.asarray(ring.T[1])
1237   
1238    M = np.array((x,y,np.ones_like(x)))
1239    B = np.array(-(x**2+y**2))
1240    result = nl.lstsq(M.T,B)
1241    return result[1],makeParmsCircle(result[0])
1242       
1243def FitEllipse(ring):
1244    import numpy.linalg as nl
1245           
1246    def makeParmsEllipse(B):
1247        det = 4.*(1.-B[0]**2)-B[1]**2
1248        if det < 0.:
1249            print 'hyperbola!'
1250            return 0
1251        elif det == 0.:
1252            print 'parabola!'
1253            return 0
1254        cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
1255            (B[1]*B[2]-2.*(1.+B[0])*B[3])/det]
1256        phi = 0.5*atand(0.5*B[1]/B[0])
1257       
1258        a = (1.+B[0])/cosd(2*phi)
1259        b = 2.-a
1260        f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4]
1261        if f/a < 0 or f/b < 0:
1262            return 0
1263        sr1 = math.sqrt(f/a)
1264        sr2 = math.sqrt(f/b)
1265        if sr1 > sr2:
1266            sr1,sr2 = SwapXY(sr1,sr2)
1267            phi -= 90.
1268            if phi < -90.:
1269                phi += 180.
1270        return cent,phi,[sr1,sr2]
1271               
1272    ring = np.array(ring)
1273    x = np.asarray(ring.T[0])
1274    y = np.asarray(ring.T[1])
1275    M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x)))
1276    B = np.array(-(x**2+y**2))
1277    result = nl.lstsq(M.T,B)
1278    return makeParmsEllipse(result[0])
1279   
1280def FitDetector(rings,p0,wave):
1281    from scipy.optimize import leastsq
1282    def ellipseCalc(B,xyd,wave):
1283        x = xyd[0]
1284        y = xyd[1]
1285        dsp = xyd[2]
1286        dist,x0,y0,phi,tilt = B
1287        tth = 2.0*npasind(wave/(2.*dsp))
1288        ttth = nptand(tth)
1289        radius = dist*ttth
1290        stth = npsind(tth)
1291        cosb = npcosd(tilt)
1292        R1 = dist*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
1293        R0 = np.sqrt(R1*radius*cosb)
1294        zdis = R1*ttth*nptand(tilt)
1295        X = x-x0+zdis*npsind(phi)
1296        Y = y-y0-zdis*npcosd(phi)
1297        XR = X*npcosd(phi)-Y*npsind(phi)
1298        YR = X*npsind(phi)+Y*npcosd(phi)
1299        return (XR/R0)**2+(YR/R1)**2-1
1300    result = leastsq(ellipseCalc,p0,args=(rings.T,wave))
1301    return result[0]
1302           
1303def ImageLocalMax(image,w,Xpix,Ypix):
1304    w2 = w*2
1305    size = len(image)
1306    if (w < Xpix < size-w) and (w < Ypix < size-w) and image[Ypix,Xpix]:
1307        Z = image[Ypix-w:Ypix+w,Xpix-w:Xpix+w]
1308        Zmax = np.argmax(Z)
1309        Zmin = np.argmin(Z)
1310        Xpix += Zmax%w2-w
1311        Ypix += Zmax/w2-w
1312        return Xpix,Ypix,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
1313    else:
1314        return 0,0,0,0
1315   
1316def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
1317    cent,phi,radii = ellipse
1318    cphi = cosd(phi)
1319    sphi = sind(phi)
1320    ring = []
1321    for a in range(0,360,2):
1322        x = radii[0]*cosd(a)
1323        y = radii[1]*sind(a)
1324        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
1325        Y = (sphi*x+cphi*y+cent[1])*scaley
1326        X,Y,I,J = ImageLocalMax(image,pix,X,Y)     
1327        if I and J and I/J > reject:
1328            X /= scalex                         #convert to mm
1329            Y /= scaley
1330            ring.append([X,Y,dsp])
1331    if len(ring) < 45:             #want more than 1/4 of a circle
1332        return []
1333    return ring
1334   
1335def calcDist(radii,tth):
1336    stth = sind(tth)
1337    ctth = cosd(tth)
1338    ttth = tand(tth)
1339    return math.sqrt(radii[0]**4/(ttth**2*((radii[0]*ctth)**2+(radii[1]*stth)**2)))
1340   
1341def calcZdisCosB(radius,tth,radii):
1342    cosB = sinb = radii[0]**2/(radius*radii[1])
1343    if cosB > 1.:
1344        return 0.,1.
1345    else:
1346        cosb = math.sqrt(1.-sinb**2)
1347        ttth = tand(tth)
1348        zdis = radii[1]*ttth*cosb/sinb
1349        return zdis,cosB
1350   
1351def GetEllipse(dsp,data):
1352    dist = data['distance']
1353    cent = data['center']
1354    tilt = data['tilt']
1355    phi = data['rotation']
1356    radii = [0,0]
1357    tth = 2.0*asind(data['wavelength']/(2.*dsp))
1358    ttth = tand(tth)
1359    stth = sind(tth)
1360    ctth = cosd(tth)
1361    cosb = cosd(tilt)
1362    radius = dist*ttth
1363    radii[1] = dist*stth*ctth*cosb/(cosb**2-stth**2)
1364    if radii[1] > 0:
1365        radii[0] = math.sqrt(radii[1]*radius*cosb)
1366        zdis = radii[1]*ttth*tand(tilt)
1367        elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
1368        return elcent,phi,radii
1369    else:
1370        return False
1371       
1372def GetDetectorXY(dsp,azm,data):
1373    from scipy.optimize import fsolve
1374    def func(xy,*args):
1375       azm,phi,R0,R1,A,B = args
1376       cp = cosd(phi)
1377       sp = sind(phi)
1378       x,y = xy
1379       out = []
1380       out.append(y-x*tand(azm))
1381       out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
1382       return out
1383    elcent,phi,radii = GetEllipse(dsp,data)
1384    cent = data['center']
1385    tilt = data['tilt']
1386    phi = data['rotation']
1387    wave = data['wavelength']
1388    dist = data['distance']
1389    tth = 2.0*asind(wave/(2.*dsp))
1390    ttth = tand(tth)
1391    radius = dist*ttth
1392    stth = sind(tth)
1393    cosb = cosd(tilt)
1394    R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
1395    R0 = math.sqrt(R1*radius*cosb)
1396    zdis = R1*ttth*tand(tilt)
1397    A = zdis*sind(phi)
1398    B = -zdis*cosd(phi)
1399    xy0 = [radius*cosd(azm),radius*sind(azm)]
1400    xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent
1401    return xy
1402                   
1403def GetTthDspAzm(x,y,data):
1404    wave = data['wavelength']
1405    dist = data['distance']
1406    cent = data['center']
1407    tilt = data['tilt']
1408    phi = data['rotation']
1409    dx = np.array(x-cent[0],dtype=np.float32)
1410    dy = np.array(y-cent[1],dtype=np.float32)
1411    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
1412    X = np.dot(X,makeMat(phi,2))
1413    Z = np.dot(X,makeMat(tilt,0)).T[2]
1414    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
1415    dsp = wave/(2.*npsind(tth/2.))
1416    azm = npatan2d(dy,dx)
1417    return tth,azm,dsp
1418   
1419def GetTth(x,y,data):
1420    return GetTthDspAzm(x,y,data)[0]
1421   
1422def GetTthAzm(x,y,data):
1423    return GetTthDspAzm(x,y,data)[0:2]
1424   
1425def GetDsp(x,y,data):
1426    return GetTthDspAzm(x,y,data)[2]
1427       
1428def ImageCompress(image,scale):
1429    if scale == 1:
1430        return image
1431    else:
1432        return image[::scale,::scale]
1433       
1434def ImageCalibrate(self,data):
1435    import copy
1436    import ImageCalibrants as calFile
1437    print 'image calibrate'
1438    ring = data['ring']
1439    pixelSize = data['pixelSize']
1440    scalex = 1000./pixelSize[0]
1441    scaley = 1000./pixelSize[1]
1442    cutoff = data['cutoff']
1443    if len(ring) < 5:
1444        print 'not enough inner ring points for ellipse'
1445        return False
1446       
1447    #fit start points on inner ring
1448    data['ellipses'] = []
1449    outE = FitRing(ring)
1450    if outE:
1451        print 'start ellipse:',outE
1452        ellipse = outE
1453    else:
1454        return False
1455       
1456    #setup 180 points on that ring for "good" fit
1457    Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ)
1458    if Ring:
1459        ellipse = FitRing(Ring)
1460        Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ)    #do again
1461        ellipse = FitRing(Ring)
1462    else:
1463        print '1st ring not sufficiently complete to proceed'
1464        return False
1465    print 'inner ring:',ellipse
1466    data['center'] = copy.copy(ellipse[0])           #not right!! (but useful for now)
1467    data['ellipses'].append(ellipse[:]+('r',))
1468    G2plt.PlotImage(self)
1469   
1470    #setup for calibration
1471    data['rings'] = []
1472    data['ellipses'] = []
1473    if not data['calibrant']:
1474        print 'no calibration material selected'
1475        return True
1476       
1477    Bravais,cell = calFile.Calibrants[data['calibrant']]
1478    A = cell2A(cell)
1479    wave = data['wavelength']
1480    cent = data['center']
1481    pixLimit = data['pixLimit']
1482    elcent,phi,radii = ellipse
1483    HKL = GenHBravais(0.5,Bravais,A)
1484    dsp = HKL[0][3]
1485    tth = 2.0*asind(wave/(2.*dsp))
1486    ttth = tand(tth)
1487    data['distance'] = dist = calcDist(radii,tth)
1488    radius = dist*tand(tth)
1489    zdis,cosB = calcZdisCosB(radius,tth,radii)
1490    cent1 = []
1491    cent2 = []
1492    xSum = 0
1493    ySum = 0
1494    zxSum = 0
1495    zySum = 0
1496    phiSum = 0
1497    tiltSum = 0
1498    distSum = 0
1499    Zsum = 0
1500    for i,H in enumerate(HKL):
1501        dsp = H[3]
1502        tth = 2.0*asind(0.5*wave/dsp)
1503        stth = sind(tth)
1504        ctth = cosd(tth)
1505        ttth = tand(tth)
1506        radius = dist*ttth
1507        elcent,phi,radii = ellipse
1508        radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
1509        radii[0] = math.sqrt(radii[1]*radius*cosB)
1510        zdis,cosB = calcZdisCosB(radius,tth,radii)
1511        zsinp = zdis*sind(phi)
1512        zcosp = zdis*cosd(phi)
1513        cent = data['center']
1514        elcent = [cent[0]+zsinp,cent[1]-zcosp]
1515        ratio = radii[1]/radii[0]
1516        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
1517        if Ring:
1518            numZ = len(Ring)
1519            data['rings'].append(np.array(Ring))
1520            ellipse = FitRing(Ring)
1521            elcent,phi,radii = ellipse               
1522            if abs(phi) > 45. and phi < 0.:
1523                phi += 180.
1524            dist = calcDist(radii,tth)
1525            distR = 1.-dist/data['distance']
1526            if distR > 0.001:
1527                print 'Wavelength too large?'
1528            elif distR < -0.001:
1529                print 'Wavelength too small?'
1530            else:
1531                if abs((radii[1]/radii[0]-ratio)/ratio) > 0.01:
1532                    print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i)
1533                    return False
1534            zdis,cosB = calcZdisCosB(radius,tth,radii)
1535            Tilt = acosd(cosB)          # 0 <= tilt <= 90
1536            zsinp = zdis*sind(ellipse[1])
1537            zcosp = zdis*cosd(ellipse[1])
1538            cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp]))
1539            cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp]))
1540            if i:
1541                d1 = cent1[-1]-cent1[-2]        #get shift of 2 possible center solutions
1542                d2 = cent2[-1]-cent2[-2]
1543                if np.dot(d2,d2) > np.dot(d1,d1):  #right solution is the larger shift
1544                    data['center'] = cent1[-1]
1545                else:
1546                    data['center'] = cent2[-1]
1547                Zsum += numZ
1548                phiSum += numZ*phi
1549                distSum += numZ*dist
1550                xSum += numZ*data['center'][0]
1551                ySum += numZ*data['center'][1]
1552                tiltSum += numZ*abs(Tilt)
1553            cent = data['center']
1554            print 'for ring # %2i dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d' \
1555                %(i,dist,phi,Tilt,cent[0],cent[1],numZ)
1556            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
1557            G2plt.PlotImage(self)
1558        else:
1559            break
1560    fullSize = len(self.ImageZ)/scalex
1561    if 2*radii[1] < .9*fullSize:
1562        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
1563    if not Zsum:
1564        print 'Only one ring fitted. Check your wavelength.'
1565        return False
1566    cent = data['center'] = [xSum/Zsum,ySum/Zsum]
1567    wave = data['wavelength']
1568    dist = data['distance'] = distSum/Zsum
1569   
1570    #possible error if no. of rings < 3! Might need trap here
1571    d1 = cent1[-1]-cent1[1]             #compare last ring to 2nd ring
1572    d2 = cent2[-1]-cent2[1]
1573    Zsign = 1
1574    len1 = math.sqrt(np.dot(d1,d1))
1575    len2 = math.sqrt(np.dot(d2,d2))
1576    t1 = d1/len1
1577    t2 = d2/len2
1578    if len2 > len1:
1579        if -135. < atan2d(t2[1],t2[0]) < 45.:
1580            Zsign = -1
1581    else:
1582        if -135. < atan2d(t1[1],t1[0]) < 45.:
1583            Zsign = -1
1584   
1585    tilt = data['tilt'] = Zsign*tiltSum/Zsum
1586    phi = data['rotation'] = phiSum/Zsum
1587    rings = np.concatenate((data['rings']),axis=0)
1588    p0 = [dist,cent[0],cent[1],phi,tilt]
1589    result = FitDetector(rings,p0,wave)
1590    data['distance'] = result[0]
1591    data['center'] = result[1:3]
1592    data['rotation'] = np.mod(result[3],360.0)
1593    data['tilt'] = result[4]
1594    N = len(data['ellipses'])
1595    data['ellipses'] = []           #clear away individual ellipse fits
1596    for H in HKL[:N]:
1597        ellipse = GetEllipse(H[3],data)
1598        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
1599    G2plt.PlotImage(self)       
1600    return True
1601   
1602def ImageIntegrate(self,data):
1603    print 'image integrate'
1604    pixelSize = data['pixelSize']
1605    scalex = pixelSize[0]/1000.
1606    scaley = pixelSize[1]/1000.
1607    LUtth = data['IOtth']
1608    if data['fullIntegrate']:
1609        LRazm = [0,360]
1610    else:
1611        LRazm = data['LRazimuth']
1612    numAzms = data['outAzimuths']
1613    numChans = data['outChannels']
1614    outGrid = np.zeros(shape=(numAzms,numChans))
1615    outNum = np.zeros(shape=(numAzms,numChans))           
1616    imageN = len(self.ImageZ)
1617    t0 = time.time()
1618    print 'Create ',imageN,' X ',imageN,' 2-theta,azimuth map'
1619    tax,tay = np.mgrid[0:imageN,0:imageN]
1620    tax = np.asfarray(tax)
1621    tay = np.asfarray(tay)
1622    tax *= scalex
1623    tay *= scaley
1624    t1 = time.time()
1625    print "Elapsed time:","%8.3f"%(t1-t0), "s"
1626    print 'Fill map with 2-theta/azimuth values'
1627    self.TA = np.reshape(GetTthAzm(tay,tax,data),(2,imageN,imageN))
1628    t2 = time.time()
1629    print "Elapsed time:","%8.3f"%(t2-t1), "s"
1630    G2plt.PlotTRImage(self)
1631    print 'Form 1-D histograms for ',numAzms,' azimuthal angles'
1632    print 'Integration limits:',LUtth,LRazm
1633    NST = np.histogram2d(self.TA[1].flatten(),self.TA[0].flatten(), \
1634        bins=(numChans,numAzms),range=[LUtth,LRazm])
1635    HST = np.histogram2d(self.TA[1].flatten(),self.TA[0].flatten(),normed=True, \
1636        bins=(numChans,numAzms),weights=self.ImageZ.flatten(),range=[LUtth,LRazm])
1637    t3 = time.time()
1638    print "Elapsed time:","%8.3f"%(t3-t2), "s"
1639    print NST[0]
1640    print HST[0]
1641    print HST[1]
1642    print HST[2]
1643   
1644       
1645def test():
1646    cell = [5,5,5,90,90,90]
1647    A = cell2A(cell)
1648    assert ( calc_V(A) == 125. )
1649   
1650    print 'test passed'
1651
1652print __name__
1653if __name__ == '__main__':
1654    test()
1655   
1656       
Note: See TracBrowser for help on using the repository browser.