source: trunk/GSASIIcomp.py @ 47

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

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